tests/mpn/t-hgcd.c

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