tests/mpz/t-gcd.c

   1 /* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui.
   2
   3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004, 2005,
   4 2008, 2009 Free 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 void one_test __GMP_PROTO ((mpz_t, mpz_t, mpz_t, int));
  29 void debug_mp __GMP_PROTO ((mpz_t, int));
  30
  31 static int gcdext_valid_p __GMP_PROTO ((const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s));
  32
  33 void
  34 check_data (void)
  35 {
  36   static const struct {
  37     const char *a;
  38     const char *b;
  39     const char *want;
  40   } data[] = {
  41     /* This tickled a bug in gmp 4.1.2 mpn/x86/k6/gcd_finda.asm. */
  42     { "0x3FFC000007FFFFFFFFFF00000000003F83FFFFFFFFFFFFFFF80000000000000001",
  43       "0x1FFE0007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFC000000000000000000000001",
  44       "5" }
  45   };
  46
  47   mpz_t  a, b, got, want;
  48   int    i;
  49
  50   mpz_init (a);
  51   mpz_init (b);
  52   mpz_init (got);
  53   mpz_init (want);
  54
  55   for (i = 0; i < numberof (data); i++)
  56     {
  57       mpz_set_str_or_abort (a, data[i].a, 0);
  58       mpz_set_str_or_abort (b, data[i].b, 0);
  59       mpz_set_str_or_abort (want, data[i].want, 0);
  60       mpz_gcd (got, a, b);
  61       MPZ_CHECK_FORMAT (got);
  62       if (mpz_cmp (got, want) != 0)
  63         {
  64           printf    ("mpz_gcd wrong on data[%d]\n", i);
  65           printf    (" a  %s\n", data[i].a);
  66           printf    (" b  %s\n", data[i].b);
  67           mpz_trace (" a", a);
  68           mpz_trace (" b", b);
  69           mpz_trace (" want", want);
  70           mpz_trace (" got ", got);
  71           abort ();
  72         }
  73     }
  74
  75   mpz_clear (a);
  76   mpz_clear (b);
  77   mpz_clear (got);
  78   mpz_clear (want);
  79 }
  80
  81 /* Keep one_test's variables global, so that we don't need
  82    to reinitialize them for each test.  */
  83 mpz_t gcd1, gcd2, s, t, temp1, temp2, temp3;
  84
  85 #if GCD_DC_THRESHOLD > GCDEXT_DC_THRESHOLD
  86 #define MAX_SCHOENHAGE_THRESHOLD GCD_DC_THRESHOLD
  87 #else
  88 #define MAX_SCHOENHAGE_THRESHOLD GCDEXT_DC_THRESHOLD
  89 #endif
  90
  91 /* Define this to make all operands be large enough for Schoenhage gcd
  92    to be used.  */
  93 #ifndef WHACK_SCHOENHAGE
  94 #define WHACK_SCHOENHAGE 0
  95 #endif
  96
  97 #if WHACK_SCHOENHAGE
  98 #define MIN_OPERAND_BITSIZE (MAX_SCHOENHAGE_THRESHOLD * GMP_NUMB_BITS)
  99 #else
 100 #define MIN_OPERAND_BITSIZE 1
 101 #endif
 102
 103 int
 104 main (int argc, char **argv)
 105 {
 106   mpz_t op1, op2, ref;
 107   int i, j, chain_len;
 108   gmp_randstate_ptr rands;
 109   mpz_t bs;
 110   unsigned long bsi, size_range;
 111   int reps = 200;
 112
 113   tests_start ();
 114   TESTS_REPS (reps, argv, argc);
 115
 116   rands = RANDS;
 117
 118   check_data ();
 119
 120   mpz_init (bs);
 121   mpz_init (op1);
 122   mpz_init (op2);
 123   mpz_init (ref);
 124   mpz_init (gcd1);
 125   mpz_init (gcd2);
 126   mpz_init (temp1);
 127   mpz_init (temp2);
 128   mpz_init (temp3);
 129   mpz_init (s);
 130   mpz_init (t);
 131
 132   /* Testcase to exercise the u0 == u1 case in mpn_gcdext_lehmer_n. */
 133   mpz_set_ui (op2, GMP_NUMB_MAX);
 134   mpz_mul_2exp (op1, op2, 100);
 135   mpz_add (op1, op1, op2);
 136   mpz_mul_ui (op2, op2, 2);
 137   one_test (op1, op2, NULL, -1);
 138
 139 #if 0
 140   mpz_set_str (op1, "4da8e405e0d2f70d6d679d3de08a5100a81ec2cff40f97b313ae75e1183f1df2b244e194ebb02a4ece50d943640a301f0f6cc7f539117b783c3f3a3f91649f8a00d2e1444d52722810562bce02fccdbbc8fe3276646e306e723dd3b", 16);
 141   mpz_set_str (op2, "76429e12e4fdd8929d89c21657097fbac09d1dc08cf7f1323a34e78ca34226e1a7a29b86fee0fa7fe2cc2a183d46d50df1fe7029590974ad7da77605f35f902cb8b9b8d22dd881eaae5919675d49a337145a029c3b33fc2b0", 16);
 142   one_test (op1, op2, NULL, -1);
 143 #endif
 144
 145   for (i = 0; i < reps; i++)
 146     {
 147       /* Generate plain operands with unknown gcd.  These types of operands
 148          have proven to trigger certain bugs in development versions of the
 149          gcd code.  The "hgcd->row[3].rsize > M" ASSERT is not triggered by
 150          the division chain code below, but that is most likely just a result
 151          of that other ASSERTs are triggered before it.  */
 152
 153       mpz_urandomb (bs, rands, 32);
 154       size_range = mpz_get_ui (bs) % 17 + 2;
 155
 156       mpz_urandomb (bs, rands, size_range);
 157       mpz_urandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
 158       mpz_urandomb (bs, rands, size_range);
 159       mpz_urandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
 160
 161       mpz_urandomb (bs, rands, 8);
 162       bsi = mpz_get_ui (bs);
 163
 164       if ((bsi & 0x3c) == 4)
 165         mpz_mul (op1, op1, op2);        /* make op1 a multiple of op2 */
 166       else if ((bsi & 0x3c) == 8)
 167         mpz_mul (op2, op1, op2);        /* make op2 a multiple of op1 */
 168
 169       if ((bsi & 1) != 0)
 170         mpz_neg (op1, op1);
 171       if ((bsi & 2) != 0)
 172         mpz_neg (op2, op2);
 173
 174       one_test (op1, op2, NULL, i);
 175
 176       /* Generate a division chain backwards, allowing otherwise unlikely huge
 177          quotients.  */
 178
 179       mpz_set_ui (op1, 0);
 180       mpz_urandomb (bs, rands, 32);
 181       mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1);
 182       mpz_rrandomb (op2, rands, mpz_get_ui (bs));
 183       mpz_add_ui (op2, op2, 1);
 184       mpz_set (ref, op2);
 185
 186 #if WHACK_SCHOENHAGE
 187       chain_len = 1000000;
 188 #else
 189       mpz_urandomb (bs, rands, 32);
 190       chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * MAX_SCHOENHAGE_THRESHOLD / 256);
 191 #endif
 192
 193       for (j = 0; j < chain_len; j++)
 194         {
 195           mpz_urandomb (bs, rands, 32);
 196           mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
 197           mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
 198           mpz_add_ui (temp2, temp2, 1);
 199           mpz_mul (temp1, op2, temp2);
 200           mpz_add (op1, op1, temp1);
 201
 202           /* Don't generate overly huge operands.  */
 203           if (SIZ (op1) > 3 * MAX_SCHOENHAGE_THRESHOLD)
 204             break;
 205
 206           mpz_urandomb (bs, rands, 32);
 207           mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
 208           mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
 209           mpz_add_ui (temp2, temp2, 1);
 210           mpz_mul (temp1, op1, temp2);
 211           mpz_add (op2, op2, temp1);
 212
 213           /* Don't generate overly huge operands.  */
 214           if (SIZ (op2) > 3 * MAX_SCHOENHAGE_THRESHOLD)
 215             break;
 216         }
 217       one_test (op1, op2, ref, i);
 218     }
 219
 220   mpz_clear (bs);
 221   mpz_clear (op1);
 222   mpz_clear (op2);
 223   mpz_clear (ref);
 224   mpz_clear (gcd1);
 225   mpz_clear (gcd2);
 226   mpz_clear (temp1);
 227   mpz_clear (temp2);
 228   mpz_clear (temp3);
 229   mpz_clear (s);
 230   mpz_clear (t);
 231
 232   tests_end ();
 233   exit (0);
 234 }
 235
 236 void
 237 debug_mp (mpz_t x, int base)
 238 {
 239   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
 240 }
 241
 242 void
 243 one_test (mpz_t op1, mpz_t op2, mpz_t ref, int i)
 244 {
 245   /*
 246   printf ("%ld %ld %ld\n", SIZ (op1), SIZ (op2), SIZ (ref));
 247   fflush (stdout);
 248   */
 249
 250   /*
 251   fprintf (stderr, "op1=");  debug_mp (op1, -16);
 252   fprintf (stderr, "op2=");  debug_mp (op2, -16);
 253   */
 254
 255   mpz_gcdext (gcd1, s, NULL, op1, op2);
 256   MPZ_CHECK_FORMAT (gcd1);
 257   MPZ_CHECK_FORMAT (s);
 258
 259   if (ref && mpz_cmp (ref, gcd1) != 0)
 260     {
 261       fprintf (stderr, "ERROR in test %d\n", i);
 262       fprintf (stderr, "mpz_gcdext returned incorrect result\n");
 263       fprintf (stderr, "op1=");                 debug_mp (op1, -16);
 264       fprintf (stderr, "op2=");                 debug_mp (op2, -16);
 265       fprintf (stderr, "expected result:\n");   debug_mp (ref, -16);
 266       fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
 267       abort ();
 268     }
 269
 270   if (!gcdext_valid_p(op1, op2, gcd1, s))
 271     {
 272       fprintf (stderr, "ERROR in test %d\n", i);
 273       fprintf (stderr, "mpz_gcdext returned invalid result\n");
 274       fprintf (stderr, "op1=");                 debug_mp (op1, -16);
 275       fprintf (stderr, "op2=");                 debug_mp (op2, -16);
 276       fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
 277       fprintf (stderr, "s=");                   debug_mp (s, -16);
 278       abort ();
 279     }
 280
 281   mpz_gcd (gcd2, op1, op2);
 282   MPZ_CHECK_FORMAT (gcd2);
 283
 284   if (mpz_cmp (gcd2, gcd1) != 0)
 285     {
 286       fprintf (stderr, "ERROR in test %d\n", i);
 287       fprintf (stderr, "mpz_gcd returned incorrect result\n");
 288       fprintf (stderr, "op1=");                 debug_mp (op1, -16);
 289       fprintf (stderr, "op2=");                 debug_mp (op2, -16);
 290       fprintf (stderr, "expected result:\n");   debug_mp (gcd1, -16);
 291       fprintf (stderr, "mpz_gcd returns:\n");   debug_mp (gcd2, -16);
 292       abort ();
 293     }
 294
 295   /* This should probably move to t-gcd_ui.c */
 296   if (mpz_fits_ulong_p (op1) || mpz_fits_ulong_p (op2))
 297     {
 298       if (mpz_fits_ulong_p (op1))
 299         mpz_gcd_ui (gcd2, op2, mpz_get_ui (op1));
 300       else
 301         mpz_gcd_ui (gcd2, op1, mpz_get_ui (op2));
 302       if (mpz_cmp (gcd2, gcd1))
 303         {
 304           fprintf (stderr, "ERROR in test %d\n", i);
 305           fprintf (stderr, "mpz_gcd_ui returned incorrect result\n");
 306           fprintf (stderr, "op1=");                 debug_mp (op1, -16);
 307           fprintf (stderr, "op2=");                 debug_mp (op2, -16);
 308           fprintf (stderr, "expected result:\n");   debug_mp (gcd1, -16);
 309           fprintf (stderr, "mpz_gcd_ui returns:\n");   debug_mp (gcd2, -16);
 310           abort ();
 311         }
 312     }
 313
 314   mpz_gcdext (gcd2, temp1, temp2, op1, op2);
 315   MPZ_CHECK_FORMAT (gcd2);
 316   MPZ_CHECK_FORMAT (temp1);
 317   MPZ_CHECK_FORMAT (temp2);
 318
 319   mpz_mul (temp1, temp1, op1);
 320   mpz_mul (temp2, temp2, op2);
 321   mpz_add (temp1, temp1, temp2);
 322
 323   if (mpz_cmp (gcd1, gcd2) != 0
 324       || mpz_cmp (gcd2, temp1) != 0)
 325     {
 326       fprintf (stderr, "ERROR in test %d\n", i);
 327       fprintf (stderr, "mpz_gcdext returned incorrect result\n");
 328       fprintf (stderr, "op1=");                 debug_mp (op1, -16);
 329       fprintf (stderr, "op2=");                 debug_mp (op2, -16);
 330       fprintf (stderr, "expected result:\n");   debug_mp (gcd1, -16);
 331       fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd2, -16);
 332       abort ();
 333     }
 334 }
 335
 336 /* Called when g is supposed to be gcd(a,b), and g = s a + t b, for some t.
 337    Uses temp1, temp2 and temp3. */
 338 static int
 339 gcdext_valid_p (const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s)
 340 {
 341   /* It's not clear that gcd(0,0) is well defined, but we allow it and require that
 342      gcd(0,0) = 0. */
 343   if (mpz_sgn (g) < 0)
 344     return 0;
 345
 346   if (mpz_sgn (a) == 0)
 347     {
 348       /* Must have g == abs (b). Any value for s is in some sense "correct",
 349          but it makes sense to require that s == 0. */
 350       return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0;
 351     }
 352   else if (mpz_sgn (b) == 0)
 353     {
 354       /* Must have g == abs (a), s == sign (a) */
 355       return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0;
 356     }
 357
 358   if (mpz_sgn (g) <= 0)
 359     return 0;
 360
 361   mpz_tdiv_qr (temp1, temp3, a, g);
 362   if (mpz_sgn (temp3) != 0)
 363     return 0;
 364
 365   mpz_tdiv_qr (temp2, temp3, b, g);
 366   if (mpz_sgn (temp3) != 0)
 367     return 0;
 368
 369   /* Require that 2 |s| < |b/g|, or |s| == 1. */
 370   if (mpz_cmpabs_ui (s, 1) > 0)
 371     {
 372       mpz_mul_2exp (temp3, s, 1);
 373       if (mpz_cmpabs (temp3, temp2) > 0)
 374         return 0;
 375     }
 376
 377   /* Compute the other cofactor. */
 378   mpz_mul(temp2, s, a);
 379   mpz_sub(temp2, g, temp2);
 380   mpz_tdiv_qr(temp2, temp3, temp2, b);
 381
 382   if (mpz_sgn (temp3) != 0)
 383     return 0;
 384
 385   /* Require that 2 |t| < |a/g| or |t| == 1*/
 386   if (mpz_cmpabs_ui (temp2, 1) > 0)
 387     {
 388       mpz_mul_2exp (temp2, temp2, 1);
 389       if (mpz_cmpabs (temp2, temp1) > 0)
 390         return 0;
 391     }
 392   return 1;
 393 }