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