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

数学计算

开发平台:

Unix_Linux

  1. /* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui.
  2. Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004 Free
  3. Software Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. #include <stdio.h>
  16. #include <stdlib.h>
  17. #include "gmp.h"
  18. #include "gmp-impl.h"
  19. #include "tests.h"
  20. static mp_size_t one_test __GMP_PROTO ((mpz_t, mpz_t, int));
  21. static void debug_mp __GMP_PROTO ((mpz_t, int));
  22. #define MIN_OPERAND_SIZE 2
  23. /* Fixed values, for regression testing of mpn_hgcd. */
  24. struct value { int res; const char *a; const char *b; };
  25. static const struct value hgcd_values[] = {
  26. #if GMP_NUMB_BITS == 32
  27.   { 5,
  28.     "0x1bddff867272a9296ac493c251d7f46f09a5591fe",
  29.     "0xb55930a2a68a916450a7de006031068c5ddb0e5c" },
  30.   { 4,
  31.     "0x2f0ece5b1ee9c15e132a01d55768dc13",
  32.     "0x1c6f4fd9873cdb24466e6d03e1cc66e7" },
  33.   { 3, "0x7FFFFC003FFFFFFFFFC5", "0x3FFFFE001FFFFFFFFFE3"},
  34. #endif
  35.   { -1, NULL, NULL }
  36. };
  37. struct hgcd_ref
  38. {
  39.   mpz_t m[2][2];
  40. };
  41. static void hgcd_ref_init __GMP_PROTO ((struct hgcd_ref *hgcd));
  42. static void hgcd_ref_clear __GMP_PROTO ((struct hgcd_ref *hgcd));
  43. static int hgcd_ref __GMP_PROTO ((struct hgcd_ref *hgcd, mpz_t a, mpz_t b));
  44. static int hgcd_ref_equal __GMP_PROTO ((const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref));
  45. int
  46. main (int argc, char **argv)
  47. {
  48.   mpz_t op1, op2, temp1, temp2;
  49.   int i, j, chain_len;
  50.   gmp_randstate_ptr rands;
  51.   mpz_t bs;
  52.   unsigned long size_range;
  53.   tests_start ();
  54.   rands = RANDS;
  55.   mpz_init (bs);
  56.   mpz_init (op1);
  57.   mpz_init (op2);
  58.   mpz_init (temp1);
  59.   mpz_init (temp2);
  60.   for (i = 0; hgcd_values[i].res >= 0; i++)
  61.     {
  62.       mp_size_t res;
  63.       mpz_set_str (op1, hgcd_values[i].a, 0);
  64.       mpz_set_str (op2, hgcd_values[i].b, 0);
  65.       res = one_test (op1, op2, -1-i);
  66.       if (res != hgcd_values[i].res)
  67. {
  68.   fprintf (stderr, "ERROR in test %dn", -1-i);
  69.   fprintf (stderr, "Bad return code from hgcdn");
  70.   fprintf (stderr, "op1=");                 debug_mp (op1, -16);
  71.   fprintf (stderr, "op2=");                 debug_mp (op2, -16);
  72.   fprintf (stderr, "expected: %dn", hgcd_values[i].res);
  73.   fprintf (stderr, "hgcd:     %dn", (int) res);
  74.   abort ();
  75. }
  76.     }
  77.   for (i = 0; i < 15; i++)
  78.     {
  79.       /* Generate plain operands with unknown gcd.  These types of operands
  80.  have proven to trigger certain bugs in development versions of the
  81.  gcd code.  The "hgcd->row[3].rsize > M" ASSERT is not triggered by
  82.  the division chain code below, but that is most likely just a result
  83.  of that other ASSERTs are triggered before it.  */
  84.       mpz_urandomb (bs, rands, 32);
  85.       size_range = mpz_get_ui (bs) % 13 + 2;
  86.       mpz_urandomb (bs, rands, size_range);
  87.       mpz_urandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
  88.       mpz_urandomb (bs, rands, size_range);
  89.       mpz_urandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
  90.       if (mpz_cmp (op1, op2) < 0)
  91. mpz_swap (op1, op2);
  92.       if (mpz_size (op1) > 0)
  93. one_test (op1, op2, i);
  94.       /* Generate a division chain backwards, allowing otherwise
  95.  unlikely huge quotients.  */
  96.       mpz_set_ui (op1, 0);
  97.       mpz_urandomb (bs, rands, 32);
  98.       mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1);
  99.       mpz_rrandomb (op2, rands, mpz_get_ui (bs));
  100.       mpz_add_ui (op2, op2, 1);
  101. #if 0
  102.       chain_len = 1000000;
  103. #else
  104.       mpz_urandomb (bs, rands, 32);
  105.       chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256);
  106. #endif
  107.       for (j = 0; j < chain_len; j++)
  108. {
  109.   mpz_urandomb (bs, rands, 32);
  110.   mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
  111.   mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
  112.   mpz_add_ui (temp2, temp2, 1);
  113.   mpz_mul (temp1, op2, temp2);
  114.   mpz_add (op1, op1, temp1);
  115.   /* Don't generate overly huge operands.  */
  116.   if (SIZ (op1) > 3 * GCD_DC_THRESHOLD)
  117.     break;
  118.   mpz_urandomb (bs, rands, 32);
  119.   mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
  120.   mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
  121.   mpz_add_ui (temp2, temp2, 1);
  122.   mpz_mul (temp1, op1, temp2);
  123.   mpz_add (op2, op2, temp1);
  124.   /* Don't generate overly huge operands.  */
  125.   if (SIZ (op2) > 3 * GCD_DC_THRESHOLD)
  126.     break;
  127. }
  128.       if (mpz_cmp (op1, op2) < 0)
  129. mpz_swap (op1, op2);
  130.       if (mpz_size (op1) > 0)
  131. one_test (op1, op2, i);
  132.     }
  133.   mpz_clear (bs);
  134.   mpz_clear (op1);
  135.   mpz_clear (op2);
  136.   mpz_clear (temp1);
  137.   mpz_clear (temp2);
  138.   tests_end ();
  139.   exit (0);
  140. }
  141. static void
  142. debug_mp (mpz_t x, int base)
  143. {
  144.   mpz_out_str (stderr, base, x); fputc ('n', stderr);
  145. }
  146. static int
  147. mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize);
  148. static mp_size_t
  149. one_test (mpz_t a, mpz_t b, int i)
  150. {
  151.   struct hgcd_matrix hgcd;
  152.   struct hgcd_ref ref;
  153.   mpz_t ref_r0;
  154.   mpz_t ref_r1;
  155.   mpz_t hgcd_r0;
  156.   mpz_t hgcd_r1;
  157.   mp_size_t res[2];
  158.   mp_size_t asize;
  159.   mp_size_t bsize;
  160.   mp_size_t hgcd_init_scratch;
  161.   mp_size_t hgcd_scratch;
  162.   mp_ptr hgcd_init_tp;
  163.   mp_ptr hgcd_tp;
  164.   asize = a->_mp_size;
  165.   bsize = b->_mp_size;
  166.   ASSERT (asize >= bsize);
  167.   hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize);
  168.   hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch);
  169.   mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp);
  170.   hgcd_scratch = mpn_hgcd_itch (asize);
  171.   hgcd_tp = refmpn_malloc_limbs (hgcd_scratch);
  172. #if 0
  173.   fprintf (stderr,
  174.    "one_test: i = %d asize = %d, bsize = %dn",
  175.    i, a->_mp_size, b->_mp_size);
  176.   gmp_fprintf (stderr,
  177.        "one_test: i = %dn"
  178.        "  a = %Zxn"
  179.        "  b = %Zxn",
  180.        i, a, b);
  181. #endif
  182.   hgcd_ref_init (&ref);
  183.   mpz_init_set (ref_r0, a);
  184.   mpz_init_set (ref_r1, b);
  185.   res[0] = hgcd_ref (&ref, ref_r0, ref_r1);
  186.   mpz_init_set (hgcd_r0, a);
  187.   mpz_init_set (hgcd_r1, b);
  188.   if (bsize < asize)
  189.     {
  190.       _mpz_realloc (hgcd_r1, asize);
  191.       MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize);
  192.     }
  193.   res[1] = mpn_hgcd (hgcd_r0->_mp_d,
  194.      hgcd_r1->_mp_d,
  195.      asize,
  196.      &hgcd, hgcd_tp);
  197.   if (res[0] != res[1])
  198.     {
  199.       fprintf (stderr, "ERROR in test %dn", i);
  200.       fprintf (stderr, "Different return value from hgcd and hgcd_refn");
  201.       fprintf (stderr, "op1=");                 debug_mp (a, -16);
  202.       fprintf (stderr, "op2=");                 debug_mp (b, -16);
  203.       fprintf (stderr, "hgcd_ref: %ldn", (long) res[0]);
  204.       fprintf (stderr, "mpn_hgcd: %ldn", (long) res[1]);
  205.       abort ();
  206.     }
  207.   if (res[0] > 0)
  208.     {
  209.       if (!hgcd_ref_equal (&hgcd, &ref)
  210.   || !mpz_mpn_equal (ref_r0, hgcd_r0->_mp_d, res[1])
  211.   || !mpz_mpn_equal (ref_r1, hgcd_r1->_mp_d, res[1]))
  212. {
  213.   fprintf (stderr, "ERROR in test %dn", i);
  214.   fprintf (stderr, "mpn_hgcd and hgcd_ref returned different valuesn");
  215.   fprintf (stderr, "op1=");                 debug_mp (a, -16);
  216.   fprintf (stderr, "op2=");                 debug_mp (b, -16);
  217.   abort ();
  218. }
  219.     }
  220.   refmpn_free_limbs (hgcd_init_tp);
  221.   refmpn_free_limbs (hgcd_tp);
  222.   hgcd_ref_clear (&ref);
  223.   mpz_clear (ref_r0);
  224.   mpz_clear (ref_r1);
  225.   mpz_clear (hgcd_r0);
  226.   mpz_clear (hgcd_r1);
  227.   return res[0];
  228. }
  229. static void
  230. hgcd_ref_init (struct hgcd_ref *hgcd)
  231. {
  232.   unsigned i;
  233.   for (i = 0; i<2; i++)
  234.     {
  235.       unsigned j;
  236.       for (j = 0; j<2; j++)
  237. mpz_init (hgcd->m[i][j]);
  238.     }
  239. }
  240. static void
  241. hgcd_ref_clear (struct hgcd_ref *hgcd)
  242. {
  243.   unsigned i;
  244.   for (i = 0; i<2; i++)
  245.     {
  246.       unsigned j;
  247.       for (j = 0; j<2; j++)
  248. mpz_clear (hgcd->m[i][j]);
  249.     }
  250. }
  251. static int
  252. sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b)
  253. {
  254.   mpz_fdiv_qr (q, r, a, b);
  255.   if (mpz_size (r) <= s)
  256.     {
  257.       mpz_add (r, r, b);
  258.       mpz_sub_ui (q, q, 1);
  259.     }
  260.   return (mpz_sgn (q) > 0);
  261. }
  262. static int
  263. hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b)
  264. {
  265.   mp_size_t n = MAX (mpz_size (a), mpz_size (b));
  266.   mp_size_t s = n/2 + 1;
  267.   mp_size_t asize;
  268.   mp_size_t bsize;
  269.   mpz_t q;
  270.   int res;
  271.   if (mpz_size (a) <= s || mpz_size (b) <= s)
  272.     return 0;
  273.   res = mpz_cmp (a, b);
  274.   if (res < 0)
  275.     {
  276.       mpz_sub (b, b, a);
  277.       if (mpz_size (b) <= s)
  278. return 0;
  279.       mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0);
  280.       mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1);
  281.     }
  282.   else if (res > 0)
  283.     {
  284.       mpz_sub (a, a, b);
  285.       if (mpz_size (a) <= s)
  286. return 0;
  287.       mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1);
  288.       mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1);
  289.     }
  290.   else
  291.     return 0;
  292.   mpz_init (q);
  293.   for (;;)
  294.     {
  295.       ASSERT (mpz_size (a) > s);
  296.       ASSERT (mpz_size (b) > s);
  297.       if (mpz_cmp (a, b) > 0)
  298. {
  299.   if (!sdiv_qr (q, a, s, a, b))
  300.     break;
  301.   mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]);
  302.   mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]);
  303. }
  304.       else
  305. {
  306.   if (!sdiv_qr (q, b, s, b, a))
  307.     break;
  308.   mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]);
  309.   mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]);
  310. }
  311.     }
  312.   mpz_clear (q);
  313.   asize = mpz_size (a);
  314.   bsize = mpz_size (b);
  315.   return MAX (asize, bsize);
  316. }
  317. static int
  318. mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize)
  319. {
  320.   mp_srcptr ap = a->_mp_d;
  321.   mp_size_t asize = a->_mp_size;
  322.   MPN_NORMALIZE (bp, bsize);
  323.   return asize == bsize && mpn_cmp (ap, bp, asize) == 0;
  324. }
  325. static int
  326. hgcd_ref_equal (const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref)
  327. {
  328.   unsigned i;
  329.   for (i = 0; i<2; i++)
  330.     {
  331.       unsigned j;
  332.       for (j = 0; j<2; j++)
  333. if (!mpz_mpn_equal (ref->m[i][j], hgcd->p[i][j], hgcd->n))
  334.   return 0;
  335.     }
  336.   return 1;
  337. }