t-bdiv.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:8k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
  2. This program is free software; you can redistribute it and/or modify it under
  3. the terms of the GNU General Public License as published by the Free Software
  4. Foundation; either version 3 of the License, or (at your option) any later
  5. version.
  6. This program is distributed in the hope that it will be useful, but WITHOUT ANY
  7. WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
  8. PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  9. You should have received a copy of the GNU General Public License along with
  10. this program.  If not, see http://www.gnu.org/licenses/.  */
  11. #include <stdlib.h> /* for strtol */
  12. #include <stdio.h> /* for printf */
  13. #include "gmp.h"
  14. #include "gmp-impl.h"
  15. #include "longlong.h"
  16. #include "tests/tests.h"
  17. static void
  18. dumpy (mp_srcptr p, mp_size_t n)
  19. {
  20.   mp_size_t i;
  21.   if (n > 20)
  22.     {
  23.       for (i = n - 1; i >= n - 4; i--)
  24. {
  25.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  26.   printf (" ");
  27. }
  28.       printf ("... ");
  29.       for (i = 3; i >= 0; i--)
  30. {
  31.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  32.   printf (" " + (i == 0));
  33. }
  34.     }
  35.   else
  36.     {
  37.       for (i = n - 1; i >= 0; i--)
  38. {
  39.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  40.   printf (" " + (i == 0));
  41. }
  42.     }
  43.   puts ("");
  44. }
  45. static unsigned long test;
  46. void
  47. check_one (mp_ptr qp, mp_srcptr rp, mp_limb_t rh,
  48.    mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, char *fname)
  49. {
  50.   mp_size_t qn;
  51.   int cmp;
  52.   mp_ptr tp;
  53.   mp_limb_t cy = 4711; /* silence warnings */
  54.   TMP_DECL;
  55.   qn = nn - dn;
  56.   if (qn == 0)
  57.     return;
  58.   TMP_MARK;
  59.   tp = TMP_ALLOC_LIMBS (nn + 1);
  60.   if (dn >= qn)
  61.     mpn_mul (tp, dp, dn, qp, qn);
  62.   else
  63.     mpn_mul (tp, qp, qn, dp, dn);
  64.   if (rp != NULL)
  65.     {
  66.       cy = mpn_add_n (tp + qn, tp + qn, rp, dn);
  67.       cmp = cy != rh || mpn_cmp (tp, np, nn) != 0;
  68.     }
  69.   else
  70.     cmp = mpn_cmp (tp, np, nn - dn) != 0;
  71.   if (cmp != 0)
  72.     {
  73.       printf ("r*******************************************************************************n");
  74.       printf ("%s inconsistent in test %lun", fname, test);
  75.       printf ("N=   "); dumpy (np, nn);
  76.       printf ("D=   "); dumpy (dp, dn);
  77.       printf ("Q=   "); dumpy (qp, qn);
  78.       if (rp != NULL)
  79. {
  80.   printf ("R=   "); dumpy (rp, dn);
  81.   printf ("Rb=  %d, Cy=%dn", (int) cy, (int) rh);
  82. }
  83.       printf ("T=   "); dumpy (tp, nn);
  84.       printf ("nn = %ld, dn = %ld, qn = %ld", nn, dn, qn);
  85.       printf ("n*******************************************************************************n");
  86.       abort ();
  87.     }
  88.   TMP_FREE;
  89. }
  90. /* These are *bit* sizes. */
  91. #define SIZE_LOG 16
  92. #define MAX_DN (1L << SIZE_LOG)
  93. #define MAX_NN (1L << (SIZE_LOG + 1))
  94. #define COUNT 500
  95. mp_limb_t
  96. random_word (gmp_randstate_ptr rs)
  97. {
  98.   mpz_t x;
  99.   mp_limb_t r;
  100.   TMP_DECL;
  101.   TMP_MARK;
  102.   MPZ_TMP_INIT (x, 2);
  103.   mpz_urandomb (x, rs, 32);
  104.   r = mpz_get_ui (x);
  105.   TMP_FREE;
  106.   return r;
  107. }
  108. int
  109. main (int argc, char **argv)
  110. {
  111.   gmp_randstate_ptr rands;
  112.   unsigned long maxnbits, maxdbits, nbits, dbits;
  113.   mpz_t n, d, tz;
  114.   mp_size_t maxnn, maxdn, nn, dn, clearn, i;
  115.   mp_ptr np, dp, qp, rp;
  116.   mp_limb_t rh;
  117.   mp_limb_t t;
  118.   mp_limb_t dinv;
  119.   int count = COUNT;
  120.   mp_ptr scratch;
  121.   mp_limb_t ran;
  122.   mp_size_t alloc, itch;
  123.   mp_limb_t rran0, rran1, qran0, qran1;
  124.   TMP_DECL;
  125.   if (argc > 1)
  126.     {
  127.       char *end;
  128.       count = strtol (argv[1], &end, 0);
  129.       if (*end || count <= 0)
  130. {
  131.   fprintf (stderr, "Invalid test count: %s.n", argv[1]);
  132.   return 1;
  133. }
  134.     }
  135.   maxdbits = MAX_DN;
  136.   maxnbits = MAX_NN;
  137.   tests_start ();
  138.   rands = RANDS;
  139.   mpz_init (n);
  140.   mpz_init (d);
  141.   mpz_init (tz);
  142.   maxnn = maxnbits / GMP_NUMB_BITS + 1;
  143.   maxdn = maxdbits / GMP_NUMB_BITS + 1;
  144.   TMP_MARK;
  145.   qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  146.   rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  147.   alloc = 1;
  148.   scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);
  149.   for (test = 0; test < count;)
  150.     {
  151.       nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;
  152.       if (maxdbits > nbits)
  153. dbits = random_word (rands) % nbits + 1;
  154.       else
  155. dbits = random_word (rands) % maxdbits + 1;
  156. #if RAND_UNIFORM
  157. #define RANDFUNC mpz_urandomb
  158. #else
  159. #define RANDFUNC mpz_rrandomb
  160. #endif
  161.       do
  162. {
  163.   RANDFUNC (n, rands, nbits);
  164.   do
  165.     {
  166.       RANDFUNC (d, rands, dbits);
  167.     }
  168.   while (mpz_sgn (d) == 0);
  169.   np = PTR (n);
  170.   dp = PTR (d);
  171.   nn = SIZ (n);
  172.   dn = SIZ (d);
  173. }
  174.       while (nn < dn);
  175.       dp[0] |= 1;
  176.       mpz_urandomb (tz, rands, 32);
  177.       t = mpz_get_ui (tz);
  178.       if (t % 17 == 0)
  179. dp[0] = GMP_NUMB_MAX;
  180.       switch (t % 16)
  181. {
  182. case 0:
  183.   clearn = random_word (rands) % nn;
  184.   for (i = 0; i <= clearn; i++)
  185.     np[i] = 0;
  186.   break;
  187. case 1:
  188.   mpn_sub_1 (np + nn - dn, dp, dn, random_word (rands));
  189.   break;
  190. case 2:
  191.   mpn_add_1 (np + nn - dn, dp, dn, random_word (rands));
  192.   break;
  193. }
  194.       test++;
  195.       binvert_limb (dinv, dp[0]);
  196.       rran0 = random_word (rands);
  197.       rran1 = random_word (rands);
  198.       qran0 = random_word (rands);
  199.       qran1 = random_word (rands);
  200.       qp[-1] = qran0;
  201.       qp[nn - dn + 1] = qran1;
  202.       rp[-1] = rran0;
  203.       ran = random_word (rands);
  204.       if ((double) (nn - dn) * dn < 1e5)
  205. {
  206.   if (nn > dn)
  207.     {
  208.       /* Test mpn_sbpi1_bdiv_qr */
  209.       MPN_ZERO (qp, nn - dn);
  210.       MPN_ZERO (rp, dn);
  211.       MPN_COPY (rp, np, nn);
  212.       rh = mpn_sbpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
  213.       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  214.       ASSERT_ALWAYS (rp[-1] == rran0);
  215.       check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_sbpi1_bdiv_qr");
  216.     }
  217.   if (nn > dn)
  218.     {
  219.       /* Test mpn_sbpi1_bdiv_q */
  220.       MPN_COPY (rp, np, nn);
  221.       MPN_ZERO (qp, nn - dn);
  222.       mpn_sbpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
  223.       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  224.       ASSERT_ALWAYS (rp[-1] == rran0);
  225.       check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_sbpi1_bdiv_q");
  226.     }
  227. }
  228.       if (dn >= 4 && nn - dn >= 2)
  229. {
  230.   /* Test mpn_dcpi1_bdiv_qr */
  231.   MPN_COPY (rp, np, nn);
  232.   MPN_ZERO (qp, nn - dn);
  233.   rh = mpn_dcpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
  234.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  235.   ASSERT_ALWAYS (rp[-1] == rran0);
  236.   check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_dcpi1_bdiv_qr");
  237. }
  238.       if (dn >= 4 && nn - dn >= 2)
  239. {
  240.   /* Test mpn_dcpi1_bdiv_q */
  241.   MPN_COPY (rp, np, nn);
  242.   MPN_ZERO (qp, nn - dn);
  243.   mpn_dcpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
  244.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  245.   ASSERT_ALWAYS (rp[-1] == rran0);
  246.   check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_dcpi1_bdiv_q");
  247. }
  248.       if (nn - dn < 2 || dn < 2)
  249. continue;
  250.       /* Test mpn_mu_bdiv_qr */
  251.       itch = mpn_mu_bdiv_qr_itch (nn, dn);
  252.       if (itch + 1 > alloc)
  253. {
  254.   scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  255.   alloc = itch + 1;
  256. }
  257.       scratch[itch] = ran;
  258.       MPN_ZERO (qp, nn - dn);
  259.       MPN_ZERO (rp, dn);
  260.       rp[dn] = rran1;
  261.       rh = mpn_mu_bdiv_qr (qp, rp, np, nn, dp, dn, scratch);
  262.       ASSERT_ALWAYS (ran == scratch[itch]);
  263.       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  264.       ASSERT_ALWAYS (rp[-1] == rran0);  ASSERT_ALWAYS (rp[dn] == rran1);
  265.       check_one (qp, rp, rh, np, nn, dp, dn, "mpn_mu_bdiv_qr");
  266.       /* Test mpn_mu_bdiv_q */
  267.       itch = mpn_mu_bdiv_q_itch (nn, dn);
  268.       if (itch + 1 > alloc)
  269. {
  270.   scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  271.   alloc = itch + 1;
  272. }
  273.       scratch[itch] = ran;
  274.       MPN_ZERO (qp, nn - dn + 1);
  275.       mpn_mu_bdiv_q (qp, np, nn - dn, dp, dn, scratch);
  276.       ASSERT_ALWAYS (ran == scratch[itch]);
  277.       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  278.       check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_mu_bdiv_q");
  279.     }
  280.   __GMP_FREE_FUNC_LIMBS (scratch, alloc);
  281.   TMP_FREE;
  282.   mpz_clear (n);
  283.   mpz_clear (d);
  284.   mpz_clear (tz);
  285.   tests_end ();
  286.   return 0;
  287. }