开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 数学计算 > 查看源码
t-div.c
资源名称：gmp-5.0.1.tar.gz　[点击查看]
上传用户：qaz666999
上传日期：2022-08-06
资源大小：2570k
文件大小：11k
源码类别：

数学计算

开发平台：

Unix_Linux

t-div.c：源码内容
  1. /* Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
  3. This program is free software; you can redistribute it and/or modify it under
  4. the terms of the GNU General Public License as published by the Free Software
  5. Foundation; either version 3 of the License, or (at your option) any later
  6. version.
  8. This program is distributed in the hope that it will be useful, but WITHOUT ANY
  9. WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
  10. PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  12. You should have received a copy of the GNU General Public License along with
  13. this program.  If not, see http://www.gnu.org/licenses/.  */
  16. #include <stdlib.h> /* for strtol */
  17. #include <stdio.h> /* for printf */
  19. #include "gmp.h"
  20. #include "gmp-impl.h"
  21. #include "longlong.h"
  22. #include "tests/tests.h"
  25. static void
  26. dumpy (mp_srcptr p, mp_size_t n)
  27. {
  28.   mp_size_t i;
  29.   if (n > 20)
  30.     {
  31.       for (i = n - 1; i >= n - 4; i--)
  32. {
  33.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  34.   printf (" ");
  35. }
  36.       printf ("... ");
  37.       for (i = 3; i >= 0; i--)
  38. {
  39.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  40.   printf (" " + (i == 0));
  41. }
  42.     }
  43.   else
  44.     {
  45.       for (i = n - 1; i >= 0; i--)
  46. {
  47.   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
  48.   printf (" " + (i == 0));
  49. }
  50.     }
  51.   puts ("");
  52. }
  54. static unsigned long test;
  56. static void
  57. check_one (mp_ptr qp, mp_srcptr rp,
  58.    mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn,
  59.    char *fname, mp_limb_t q_allowed_err)
  60. {
  61.   mp_size_t qn = nn - dn + 1;
  62.   mp_ptr tp;
  63.   const char *msg;
  64.   const char *tvalue;
  65.   mp_limb_t i;
  66.   TMP_DECL;
  67.   TMP_MARK;
  69.   tp = TMP_ALLOC_LIMBS (nn + 1);
  70.   if (dn >= qn)
  71.     refmpn_mul (tp, dp, dn, qp, qn);
  72.   else
  73.     refmpn_mul (tp, qp, qn, dp, dn);
  75.   for (i = 0; i < q_allowed_err && (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0); i++)
  76.     ASSERT_NOCARRY (refmpn_sub (tp, tp, nn+1, dp, dn));
  78.   if (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0)
  79.     {
  80.       msg = "q too large";
  81.       tvalue = "Q*D";
  82.     error:
  83.       printf ("r*******************************************************************************n");
  84.       printf ("%s failed test %lu: %sn", fname, test, msg);
  85.       printf ("N=    "); dumpy (np, nn);
  86.       printf ("D=    "); dumpy (dp, dn);
  87.       printf ("Q=    "); dumpy (qp, qn);
  88.       if (rp)
  89. { printf ("R=    "); dumpy (rp, dn); }
  90.       printf ("%5s=", tvalue); dumpy (tp, nn+1);
  91.       printf ("nn = %ld, dn = %ld, qn = %ldn", nn, dn, qn);
  92.       abort ();
  93.     }
  95.   ASSERT_NOCARRY (refmpn_sub_n (tp, np, tp, nn));
  96.   tvalue = "N-Q*D";
  97.   if (!mpn_zero_p (tp + dn, nn - dn) || mpn_cmp (tp, dp, dn) >= 0)
  98.     {
  99.       msg = "q too small";
  100.       goto error;
  101.     }
  103.   if (rp && mpn_cmp (rp, tp, dn) != 0)
  104.     {
  105.       msg = "r incorrect";
  106.       goto error;
  107.     }
  109.   TMP_FREE;
  110. }
  113. /* These are *bit* sizes. */
  114. #ifndef SIZE_LOG
  115. #define SIZE_LOG 17
  116. #endif
  117. #define MAX_DN (1L << SIZE_LOG)
  118. #define MAX_NN (1L << (SIZE_LOG + 1))
  120. #define COUNT 200
  122. mp_limb_t
  123. random_word (gmp_randstate_ptr rs)
  124. {
  125.   mpz_t x;
  126.   mp_limb_t r;
  127.   TMP_DECL;
  128.   TMP_MARK;
  130.   MPZ_TMP_INIT (x, 2);
  131.   mpz_urandomb (x, rs, 32);
  132.   r = mpz_get_ui (x);
  133.   TMP_FREE;
  134.   return r;
  135. }
  137. int
  138. main (int argc, char **argv)
  139. {
  140.   gmp_randstate_ptr rands;
  141.   unsigned long maxnbits, maxdbits, nbits, dbits;
  142.   mpz_t n, d, q, r, tz;
  143.   mp_size_t maxnn, maxdn, nn, dn, clearn, i;
  144.   mp_ptr np, dp, qp, rp;
  145.   mp_limb_t t;
  146.   gmp_pi1_t dinv;
  147.   int count = COUNT;
  148.   mp_ptr scratch;
  149.   mp_limb_t ran;
  150.   mp_size_t alloc, itch;
  151.   mp_limb_t rran0, rran1, qran0, qran1;
  152.   TMP_DECL;
  154.   if (argc > 1)
  155.     {
  156.       char *end;
  157.       count = strtol (argv[1], &end, 0);
  158.       if (*end || count <= 0)
  159. {
  160.   fprintf (stderr, "Invalid test count: %s.n", argv[1]);
  161.   return 1;
  162. }
  163.     }
  166.   maxdbits = MAX_DN;
  167.   maxnbits = MAX_NN;
  169.   tests_start ();
  170.   rands = RANDS;
  172.   mpz_init (n);
  173.   mpz_init (d);
  174.   mpz_init (q);
  175.   mpz_init (r);
  176.   mpz_init (tz);
  178.   maxnn = maxnbits / GMP_NUMB_BITS + 1;
  179.   maxdn = maxdbits / GMP_NUMB_BITS + 1;
  181.   TMP_MARK;
  183.   qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  184.   rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  186.   alloc = 1;
  187.   scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);
  189.   for (test = 0; test < count;)
  190.     {
  191.       do
  192. {
  193.   nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;
  194.   if (maxdbits > nbits)
  195.     dbits = random_word (rands) % nbits + 1;
  196.   else
  197.     dbits = random_word (rands) % maxdbits + 1;
  198. }
  199.       while (nbits < dbits);
  201. #if RAND_UNIFORM
  202. #define RANDFUNC mpz_urandomb
  203. #else
  204. #define RANDFUNC mpz_rrandomb
  205. #endif
  207.       do
  208. RANDFUNC (d, rands, dbits);
  209.       while (mpz_sgn (d) == 0);
  210.       dn = SIZ (d);
  211.       dp = PTR (d);
  212.       dp[dn - 1] |= GMP_NUMB_HIGHBIT;
  214.       if (test % 2 == 0)
  215. {
  216.   RANDFUNC (n, rands, nbits);
  217.   nn = SIZ (n);
  218.   ASSERT_ALWAYS (nn >= dn);
  219. }
  220.       else
  221. {
  222.   do
  223.     {
  224.       RANDFUNC (q, rands, random_word (rands) % (nbits - dbits + 1));
  225.       RANDFUNC (r, rands, random_word (rands) % mpz_sizeinbase (d, 2));
  226.       mpz_mul (n, q, d);
  227.       mpz_add (n, n, r);
  228.       nn = SIZ (n);
  229.     }
  230.   while (nn > maxnn || nn < dn);
  231. }
  233.       ASSERT_ALWAYS (nn <= maxnn);
  234.       ASSERT_ALWAYS (dn <= maxdn);
  236.       np = PTR (n);
  238.       mpz_urandomb (tz, rands, 32);
  239.       t = mpz_get_ui (tz);
  241.       if (t % 17 == 0)
  242. dp[dn - 1] = GMP_NUMB_MAX;
  244.       switch (t % 16)
  245. {
  246. case 0:
  247.   clearn = random_word (rands) % nn;
  248.   for (i = clearn; i < nn; i++)
  249.     np[i] = 0;
  250.   break;
  251. case 1:
  252.   mpn_sub_1 (np + nn - dn, dp, dn, random_word (rands));
  253.   break;
  254. case 2:
  255.   mpn_add_1 (np + nn - dn, dp, dn, random_word (rands));
  256.   break;
  257. }
  259.       test++;
  261.       invert_pi1 (dinv, dp[dn - 1], dp[dn - 2]);
  263.       rran0 = random_word (rands);
  264.       rran1 = random_word (rands);
  265.       qran0 = random_word (rands);
  266.       qran1 = random_word (rands);
  268.       qp[-1] = qran0;
  269.       qp[nn - dn + 1] = qran1;
  270.       rp[-1] = rran0;
  272.       ran = random_word (rands);
  274.       if ((double) (nn - dn) * dn < 1e5)
  275. {
  276.   /* Test mpn_sbpi1_div_qr */
  277.   if (dn > 2)
  278.     {
  279.       MPN_COPY (rp, np, nn);
  280.       if (nn > dn)
  281. MPN_ZERO (qp, nn - dn);
  282.       qp[nn - dn] = mpn_sbpi1_div_qr (qp, rp, nn, dp, dn, dinv.inv32);
  283.       check_one (qp, rp, np, nn, dp, dn, "mpn_sbpi1_div_qr", 0);
  284.     }
  286.   /* Test mpn_sbpi1_divappr_q */
  287.   if (dn > 2)
  288.     {
  289.       MPN_COPY (rp, np, nn);
  290.       if (nn > dn)
  291. MPN_ZERO (qp, nn - dn);
  292.       qp[nn - dn] = mpn_sbpi1_divappr_q (qp, rp, nn, dp, dn, dinv.inv32);
  293.       check_one (qp, NULL, np, nn, dp, dn, "mpn_sbpi1_divappr_q", 1);
  294.     }
  296.   /* Test mpn_sbpi1_div_q */
  297.   if (dn > 2)
  298.     {
  299.       MPN_COPY (rp, np, nn);
  300.       if (nn > dn)
  301. MPN_ZERO (qp, nn - dn);
  302.       qp[nn - dn] = mpn_sbpi1_div_q (qp, rp, nn, dp, dn, dinv.inv32);
  303.       check_one (qp, NULL, np, nn, dp, dn, "mpn_sbpi1_div_q", 0);
  304.     }
  305. }
  307.       /* Test mpn_dcpi1_div_qr */
  308.       if (dn >= 6 && nn - dn >= 3)
  309. {
  310.   MPN_COPY (rp, np, nn);
  311.   if (nn > dn)
  312.     MPN_ZERO (qp, nn - dn);
  313.   qp[nn - dn] = mpn_dcpi1_div_qr (qp, rp, nn, dp, dn, &dinv);
  314.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  315.   ASSERT_ALWAYS (rp[-1] == rran0);
  316.   check_one (qp, rp, np, nn, dp, dn, "mpn_dcpi1_div_qr", 0);
  317. }
  319.       /* Test mpn_dcpi1_divappr_q */
  320.       if (dn >= 6 && nn - dn >= 3)
  321. {
  322.   MPN_COPY (rp, np, nn);
  323.   if (nn > dn)
  324.     MPN_ZERO (qp, nn - dn);
  325.   qp[nn - dn] = mpn_dcpi1_divappr_q (qp, rp, nn, dp, dn, &dinv);
  326.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  327.   ASSERT_ALWAYS (rp[-1] == rran0);
  328.   check_one (qp, NULL, np, nn, dp, dn, "mpn_dcpi1_divappr_q", 1);
  329. }
  331.       /* Test mpn_dcpi1_div_q */
  332.       if (dn >= 6 && nn - dn >= 3)
  333. {
  334.   MPN_COPY (rp, np, nn);
  335.   if (nn > dn)
  336.     MPN_ZERO (qp, nn - dn);
  337.   qp[nn - dn] = mpn_dcpi1_div_q (qp, rp, nn, dp, dn, &dinv);
  338.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  339.   ASSERT_ALWAYS (rp[-1] == rran0);
  340.   check_one (qp, NULL, np, nn, dp, dn, "mpn_dcpi1_div_q", 0);
  341. }
  343.      /* Test mpn_mu_div_qr */
  344.       if (nn - dn > 2 && dn >= 2)
  345. {
  346.   itch = mpn_mu_div_qr_itch (nn, dn, 0);
  347.   if (itch + 1 > alloc)
  348.     {
  349.       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  350.       alloc = itch + 1;
  351.     }
  352.   scratch[itch] = ran;
  353.   MPN_ZERO (qp, nn - dn);
  354.   MPN_ZERO (rp, dn);
  355.   rp[dn] = rran1;
  356.   qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
  357.   ASSERT_ALWAYS (ran == scratch[itch]);
  358.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  359.   ASSERT_ALWAYS (rp[-1] == rran0);  ASSERT_ALWAYS (rp[dn] == rran1);
  360.   check_one (qp, rp, np, nn, dp, dn, "mpn_mu_div_qr", 0);
  361. }
  363.       /* Test mpn_mu_divappr_q */
  364.       if (nn - dn > 2 && dn >= 2)
  365. {
  366.   itch = mpn_mu_divappr_q_itch (nn, dn, 0);
  367.   if (itch + 1 > alloc)
  368.     {
  369.       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  370.       alloc = itch + 1;
  371.     }
  372.   scratch[itch] = ran;
  373.   MPN_ZERO (qp, nn - dn);
  374.   qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dp, dn, scratch);
  375.   ASSERT_ALWAYS (ran == scratch[itch]);
  376.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  377.   check_one (qp, NULL, np, nn, dp, dn, "mpn_mu_divappr_q", 4);
  378. }
  380.       /* Test mpn_mu_div_q */
  381.       if (nn - dn > 2 && dn >= 2)
  382. {
  383.   itch = mpn_mu_div_q_itch (nn, dn, 0);
  384.   if (itch + 1> alloc)
  385.     {
  386.       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  387.       alloc = itch + 1;
  388.     }
  389.   scratch[itch] = ran;
  390.   MPN_ZERO (qp, nn - dn);
  391.   qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dp, dn, scratch);
  392.   ASSERT_ALWAYS (ran == scratch[itch]);
  393.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  394.   check_one (qp, NULL, np, nn, dp, dn, "mpn_mu_div_q", 0);
  395. }
  398.       if (1)
  399. {
  400.   itch = nn + 1;
  401.   if (itch + 1> alloc)
  402.     {
  403.       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  404.       alloc = itch + 1;
  405.     }
  406.   scratch[itch] = ran;
  407.   mpn_div_q (qp, np, nn, dp, dn, scratch);
  408.   ASSERT_ALWAYS (ran == scratch[itch]);
  409.   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  410.   check_one (qp, NULL, np, nn, dp, dn, "mpn_div_q", 0);
  411. }
  413.       /* Finally, test mpn_div_q without msb set.  */
  414.       dp[dn - 1] &= ~GMP_NUMB_HIGHBIT;
  415.       if (dp[dn - 1] == 0)
  416. continue;
  418.       itch = nn + 1;
  419.       if (itch + 1> alloc)
  420. {
  421.   scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
  422.   alloc = itch + 1;
  423. }
  424.       scratch[itch] = ran;
  425.       mpn_div_q (qp, np, nn, dp, dn, scratch);
  426.       ASSERT_ALWAYS (ran == scratch[itch]);
  427.       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
  428.       check_one (qp, NULL, np, nn, dp, dn, "mpn_div_q", 0);
  429.     }
  431.   __GMP_FREE_FUNC_LIMBS (scratch, alloc);
  433.   TMP_FREE;
  435.   mpz_clear (n);
  436.   mpz_clear (d);
  437.   mpz_clear (q);
  438.   mpz_clear (r);
  439.   mpz_clear (tz);
  441.   tests_end ();
  442.   return 0;
  443. }