1/* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui. 2 3Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004, 2005, 42008, 2009 Free Software Foundation, Inc. 5 6This file is part of the GNU MP Library. 7 8The GNU MP Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MP Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along 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 28void one_test __GMP_PROTO ((mpz_t, mpz_t, mpz_t, int)); 29void debug_mp __GMP_PROTO ((mpz_t, int)); 30 31static int gcdext_valid_p __GMP_PROTO ((const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s)); 32 33void 34check_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. */ 83mpz_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 103int 104main (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_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE); 158 mpz_urandomb (bs, rands, size_range); 159 mpz_rrandomb (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 236void 237debug_mp (mpz_t x, int base) 238{ 239 mpz_out_str (stderr, base, x); fputc ('\n', stderr); 240} 241 242void 243one_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. */ 338static int 339gcdext_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