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