tests/mpn/t-hgcd.c

   1 /* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui.
   2
   3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004 Free
   4 Software Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://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 __GMP_PROTO ((mpz_t, mpz_t, int));
  29 static void debug_mp __GMP_PROTO ((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 __GMP_PROTO ((struct hgcd_ref *hgcd));
  54 static void hgcd_ref_clear __GMP_PROTO ((struct hgcd_ref *hgcd));
  55 static int hgcd_ref __GMP_PROTO ((struct hgcd_ref *hgcd, mpz_t a, mpz_t b));
  56 static int hgcd_ref_equal __GMP_PROTO ((const struct hgcd_matrix *hgcd, const struct hgcd_ref *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.  The "hgcd->row[3].rsize > M" ASSERT is not triggered by
 101          the division chain code below, but that is most likely just a result
 102          of that other ASSERTs are triggered before it.  */
 103
 104       mpz_urandomb (bs, rands, 32);
 105       size_range = mpz_get_ui (bs) % 13 + 2;
 106
 107       mpz_urandomb (bs, rands, size_range);
 108       mpz_urandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
 109       mpz_urandomb (bs, rands, size_range);
 110       mpz_urandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
 111
 112       if (mpz_cmp (op1, op2) < 0)
 113         mpz_swap (op1, op2);
 114
 115       if (mpz_size (op1) > 0)
 116         one_test (op1, op2, i);
 117
 118       /* Generate a division chain backwards, allowing otherwise
 119          unlikely huge quotients.  */
 120
 121       mpz_set_ui (op1, 0);
 122       mpz_urandomb (bs, rands, 32);
 123       mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1);
 124       mpz_rrandomb (op2, rands, mpz_get_ui (bs));
 125       mpz_add_ui (op2, op2, 1);
 126
 127 #if 0
 128       chain_len = 1000000;
 129 #else
 130       mpz_urandomb (bs, rands, 32);
 131       chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256);
 132 #endif
 133
 134       for (j = 0; j < chain_len; j++)
 135         {
 136           mpz_urandomb (bs, rands, 32);
 137           mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
 138           mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
 139           mpz_add_ui (temp2, temp2, 1);
 140           mpz_mul (temp1, op2, temp2);
 141           mpz_add (op1, op1, temp1);
 142
 143           /* Don't generate overly huge operands.  */
 144           if (SIZ (op1) > 3 * GCD_DC_THRESHOLD)
 145             break;
 146
 147           mpz_urandomb (bs, rands, 32);
 148           mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
 149           mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
 150           mpz_add_ui (temp2, temp2, 1);
 151           mpz_mul (temp1, op1, temp2);
 152           mpz_add (op2, op2, temp1);
 153
 154           /* Don't generate overly huge operands.  */
 155           if (SIZ (op2) > 3 * GCD_DC_THRESHOLD)
 156             break;
 157         }
 158       if (mpz_cmp (op1, op2) < 0)
 159         mpz_swap (op1, op2);
 160
 161       if (mpz_size (op1) > 0)
 162         one_test (op1, op2, i);
 163     }
 164
 165   mpz_clear (bs);
 166   mpz_clear (op1);
 167   mpz_clear (op2);
 168   mpz_clear (temp1);
 169   mpz_clear (temp2);
 170
 171   tests_end ();
 172   exit (0);
 173 }
 174
 175 static void
 176 debug_mp (mpz_t x, int base)
 177 {
 178   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
 179 }
 180
 181 static int
 182 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize);
 183
 184 static mp_size_t
 185 one_test (mpz_t a, mpz_t b, int i)
 186 {
 187   struct hgcd_matrix hgcd;
 188   struct hgcd_ref ref;
 189
 190   mpz_t ref_r0;
 191   mpz_t ref_r1;
 192   mpz_t hgcd_r0;
 193   mpz_t hgcd_r1;
 194
 195   mp_size_t res[2];
 196   mp_size_t asize;
 197   mp_size_t bsize;
 198
 199   mp_size_t hgcd_init_scratch;
 200   mp_size_t hgcd_scratch;
 201
 202   mp_ptr hgcd_init_tp;
 203   mp_ptr hgcd_tp;
 204
 205   asize = a->_mp_size;
 206   bsize = b->_mp_size;
 207
 208   ASSERT (asize >= bsize);
 209
 210   hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize);
 211   hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch);
 212   mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp);
 213
 214   hgcd_scratch = mpn_hgcd_itch (asize);
 215   hgcd_tp = refmpn_malloc_limbs (hgcd_scratch);
 216
 217 #if 0
 218   fprintf (stderr,
 219            "one_test: i = %d asize = %d, bsize = %d\n",
 220            i, a->_mp_size, b->_mp_size);
 221
 222   gmp_fprintf (stderr,
 223                "one_test: i = %d\n"
 224                "  a = %Zx\n"
 225                "  b = %Zx\n",
 226                i, a, b);
 227 #endif
 228   hgcd_ref_init (&ref);
 229
 230   mpz_init_set (ref_r0, a);
 231   mpz_init_set (ref_r1, b);
 232   res[0] = hgcd_ref (&ref, ref_r0, ref_r1);
 233
 234   mpz_init_set (hgcd_r0, a);
 235   mpz_init_set (hgcd_r1, b);
 236   if (bsize < asize)
 237     {
 238       _mpz_realloc (hgcd_r1, asize);
 239       MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize);
 240     }
 241   res[1] = mpn_hgcd (hgcd_r0->_mp_d,
 242                      hgcd_r1->_mp_d,
 243                      asize,
 244                      &hgcd, hgcd_tp);
 245
 246   if (res[0] != res[1])
 247     {
 248       fprintf (stderr, "ERROR in test %d\n", i);
 249       fprintf (stderr, "Different return value from hgcd and hgcd_ref\n");
 250       fprintf (stderr, "op1=");                 debug_mp (a, -16);
 251       fprintf (stderr, "op2=");                 debug_mp (b, -16);
 252       fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]);
 253       fprintf (stderr, "mpn_hgcd: %ld\n", (long) res[1]);
 254       abort ();
 255     }
 256   if (res[0] > 0)
 257     {
 258       if (!hgcd_ref_equal (&hgcd, &ref)
 259           || !mpz_mpn_equal (ref_r0, hgcd_r0->_mp_d, res[1])
 260           || !mpz_mpn_equal (ref_r1, hgcd_r1->_mp_d, res[1]))
 261         {
 262           fprintf (stderr, "ERROR in test %d\n", i);
 263           fprintf (stderr, "mpn_hgcd and hgcd_ref returned different values\n");
 264           fprintf (stderr, "op1=");                 debug_mp (a, -16);
 265           fprintf (stderr, "op2=");                 debug_mp (b, -16);
 266           abort ();
 267         }
 268     }
 269
 270   refmpn_free_limbs (hgcd_init_tp);
 271   refmpn_free_limbs (hgcd_tp);
 272   hgcd_ref_clear (&ref);
 273   mpz_clear (ref_r0);
 274   mpz_clear (ref_r1);
 275   mpz_clear (hgcd_r0);
 276   mpz_clear (hgcd_r1);
 277
 278   return res[0];
 279 }
 280
 281 static void
 282 hgcd_ref_init (struct hgcd_ref *hgcd)
 283 {
 284   unsigned i;
 285   for (i = 0; i<2; i++)
 286     {
 287       unsigned j;
 288       for (j = 0; j<2; j++)
 289         mpz_init (hgcd->m[i][j]);
 290     }
 291 }
 292
 293 static void
 294 hgcd_ref_clear (struct hgcd_ref *hgcd)
 295 {
 296   unsigned i;
 297   for (i = 0; i<2; i++)
 298     {
 299       unsigned j;
 300       for (j = 0; j<2; j++)
 301         mpz_clear (hgcd->m[i][j]);
 302     }
 303 }
 304
 305
 306 static int
 307 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b)
 308 {
 309   mpz_fdiv_qr (q, r, a, b);
 310   if (mpz_size (r) <= s)
 311     {
 312       mpz_add (r, r, b);
 313       mpz_sub_ui (q, q, 1);
 314     }
 315
 316   return (mpz_sgn (q) > 0);
 317 }
 318
 319 static int
 320 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b)
 321 {
 322   mp_size_t n = MAX (mpz_size (a), mpz_size (b));
 323   mp_size_t s = n/2 + 1;
 324   mp_size_t asize;
 325   mp_size_t bsize;
 326   mpz_t q;
 327   int res;
 328
 329   if (mpz_size (a) <= s || mpz_size (b) <= s)
 330     return 0;
 331
 332   res = mpz_cmp (a, b);
 333   if (res < 0)
 334     {
 335       mpz_sub (b, b, a);
 336       if (mpz_size (b) <= s)
 337         return 0;
 338
 339       mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0);
 340       mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1);
 341     }
 342   else if (res > 0)
 343     {
 344       mpz_sub (a, a, b);
 345       if (mpz_size (a) <= s)
 346         return 0;
 347
 348       mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1);
 349       mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1);
 350     }
 351   else
 352     return 0;
 353
 354   mpz_init (q);
 355
 356   for (;;)
 357     {
 358       ASSERT (mpz_size (a) > s);
 359       ASSERT (mpz_size (b) > s);
 360
 361       if (mpz_cmp (a, b) > 0)
 362         {
 363           if (!sdiv_qr (q, a, s, a, b))
 364             break;
 365           mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]);
 366           mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]);
 367         }
 368       else
 369         {
 370           if (!sdiv_qr (q, b, s, b, a))
 371             break;
 372           mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]);
 373           mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]);
 374         }
 375     }
 376
 377   mpz_clear (q);
 378
 379   asize = mpz_size (a);
 380   bsize = mpz_size (b);
 381   return MAX (asize, bsize);
 382 }
 383
 384 static int
 385 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize)
 386 {
 387   mp_srcptr ap = a->_mp_d;
 388   mp_size_t asize = a->_mp_size;
 389
 390   MPN_NORMALIZE (bp, bsize);
 391   return asize == bsize && mpn_cmp (ap, bp, asize) == 0;
 392 }
 393
 394 static int
 395 hgcd_ref_equal (const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref)
 396 {
 397   unsigned i;
 398
 399   for (i = 0; i<2; i++)
 400     {
 401       unsigned j;
 402
 403       for (j = 0; j<2; j++)
 404         if (!mpz_mpn_equal (ref->m[i][j], hgcd->p[i][j], hgcd->n))
 405           return 0;
 406     }
 407
 408   return 1;
 409 }