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 } 394