1 /* Test mpn_hgcd_appr. 2 3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004, 2011 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 http://www.gnu.org/licenses/. */ 20 21 #include <stdio.h> 22 #include <stdlib.h> 23 #include <string.h> 24 25 #include "gmp.h" 26 #include "gmp-impl.h" 27 #include "tests.h" 28 29 static mp_size_t one_test (mpz_t, mpz_t, int); 30 static void debug_mp (mpz_t, int); 31 32 #define MIN_OPERAND_SIZE 2 33 34 struct hgcd_ref 35 { 36 mpz_t m[2][2]; 37 }; 38 39 static void hgcd_ref_init (struct hgcd_ref *hgcd); 40 static void hgcd_ref_clear (struct hgcd_ref *hgcd); 41 static int hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b); 42 static int hgcd_ref_equal (const struct hgcd_ref *, const struct hgcd_ref *); 43 static int hgcd_appr_valid_p (mpz_t, mpz_t, mp_size_t, struct hgcd_ref *, 44 mpz_t, mpz_t, mp_size_t, struct hgcd_matrix *); 45 46 static int verbose_flag = 0; 47 48 int 49 main (int argc, char **argv) 50 { 51 mpz_t op1, op2, temp1, temp2; 52 int i, j, chain_len; 53 gmp_randstate_ptr rands; 54 mpz_t bs; 55 unsigned long size_range; 56 57 if (argc > 1) 58 { 59 if (strcmp (argv[1], "-v") == 0) 60 verbose_flag = 1; 61 else 62 { 63 fprintf (stderr, "Invalid argument.\n"); 64 return 1; 65 } 66 } 67 68 tests_start (); 69 rands = RANDS; 70 71 mpz_init (bs); 72 mpz_init (op1); 73 mpz_init (op2); 74 mpz_init (temp1); 75 mpz_init (temp2); 76 77 for (i = 0; i < 15; i++) 78 { 79 /* Generate plain operands with unknown gcd. These types of operands 80 have proven to trigger certain bugs in development versions of the 81 gcd code. */ 82 83 mpz_urandomb (bs, rands, 32); 84 size_range = mpz_get_ui (bs) % 13 + 2; 85 86 mpz_urandomb (bs, rands, size_range); 87 mpz_urandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); 88 mpz_urandomb (bs, rands, size_range); 89 mpz_urandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); 90 91 if (mpz_cmp (op1, op2) < 0) 92 mpz_swap (op1, op2); 93 94 if (mpz_size (op1) > 0) 95 one_test (op1, op2, i); 96 97 /* Generate a division chain backwards, allowing otherwise 98 unlikely huge quotients. */ 99 100 mpz_set_ui (op1, 0); 101 mpz_urandomb (bs, rands, 32); 102 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1); 103 mpz_rrandomb (op2, rands, mpz_get_ui (bs)); 104 mpz_add_ui (op2, op2, 1); 105 106 #if 0 107 chain_len = 1000000; 108 #else 109 mpz_urandomb (bs, rands, 32); 110 chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256); 111 #endif 112 113 for (j = 0; j < chain_len; j++) 114 { 115 mpz_urandomb (bs, rands, 32); 116 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); 117 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); 118 mpz_add_ui (temp2, temp2, 1); 119 mpz_mul (temp1, op2, temp2); 120 mpz_add (op1, op1, temp1); 121 122 /* Don't generate overly huge operands. */ 123 if (SIZ (op1) > 3 * GCD_DC_THRESHOLD) 124 break; 125 126 mpz_urandomb (bs, rands, 32); 127 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); 128 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); 129 mpz_add_ui (temp2, temp2, 1); 130 mpz_mul (temp1, op1, temp2); 131 mpz_add (op2, op2, temp1); 132 133 /* Don't generate overly huge operands. */ 134 if (SIZ (op2) > 3 * GCD_DC_THRESHOLD) 135 break; 136 } 137 if (mpz_cmp (op1, op2) < 0) 138 mpz_swap (op1, op2); 139 140 if (mpz_size (op1) > 0) 141 one_test (op1, op2, i); 142 } 143 144 mpz_clear (bs); 145 mpz_clear (op1); 146 mpz_clear (op2); 147 mpz_clear (temp1); 148 mpz_clear (temp2); 149 150 tests_end (); 151 exit (0); 152 } 153 154 static void 155 debug_mp (mpz_t x, int base) 156 { 157 mpz_out_str (stderr, base, x); fputc ('\n', stderr); 158 } 159 160 static int 161 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize); 162 163 static mp_size_t 164 one_test (mpz_t a, mpz_t b, int i) 165 { 166 struct hgcd_matrix hgcd; 167 struct hgcd_ref ref; 168 169 mpz_t ref_r0; 170 mpz_t ref_r1; 171 mpz_t hgcd_r0; 172 mpz_t hgcd_r1; 173 174 int res[2]; 175 mp_size_t asize; 176 mp_size_t bsize; 177 178 mp_size_t hgcd_init_scratch; 179 mp_size_t hgcd_scratch; 180 181 mp_ptr hgcd_init_tp; 182 mp_ptr hgcd_tp; 183 mp_limb_t marker[4]; 184 185 asize = a->_mp_size; 186 bsize = b->_mp_size; 187 188 ASSERT (asize >= bsize); 189 190 hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize); 191 hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch + 2) + 1; 192 mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp); 193 194 hgcd_scratch = mpn_hgcd_appr_itch (asize); 195 hgcd_tp = refmpn_malloc_limbs (hgcd_scratch + 2) + 1; 196 197 mpn_random (marker, 4); 198 199 hgcd_init_tp[-1] = marker[0]; 200 hgcd_init_tp[hgcd_init_scratch] = marker[1]; 201 hgcd_tp[-1] = marker[2]; 202 hgcd_tp[hgcd_scratch] = marker[3]; 203 204 #if 0 205 fprintf (stderr, 206 "one_test: i = %d asize = %d, bsize = %d\n", 207 i, a->_mp_size, b->_mp_size); 208 209 gmp_fprintf (stderr, 210 "one_test: i = %d\n" 211 " a = %Zx\n" 212 " b = %Zx\n", 213 i, a, b); 214 #endif 215 hgcd_ref_init (&ref); 216 217 mpz_init_set (ref_r0, a); 218 mpz_init_set (ref_r1, b); 219 res[0] = hgcd_ref (&ref, ref_r0, ref_r1); 220 221 mpz_init_set (hgcd_r0, a); 222 mpz_init_set (hgcd_r1, b); 223 if (bsize < asize) 224 { 225 _mpz_realloc (hgcd_r1, asize); 226 MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize); 227 } 228 res[1] = mpn_hgcd_appr (hgcd_r0->_mp_d, 229 hgcd_r1->_mp_d, 230 asize, 231 &hgcd, hgcd_tp); 232 233 if (hgcd_init_tp[-1] != marker[0] 234 || hgcd_init_tp[hgcd_init_scratch] != marker[1] 235 || hgcd_tp[-1] != marker[2] 236 || hgcd_tp[hgcd_scratch] != marker[3]) 237 { 238 fprintf (stderr, "ERROR in test %d\n", i); 239 fprintf (stderr, "scratch space overwritten!\n"); 240 241 if (hgcd_init_tp[-1] != marker[0]) 242 gmp_fprintf (stderr, 243 "before init_tp: %Mx\n" 244 "expected: %Mx\n", 245 hgcd_init_tp[-1], marker[0]); 246 if (hgcd_init_tp[hgcd_init_scratch] != marker[1]) 247 gmp_fprintf (stderr, 248 "after init_tp: %Mx\n" 249 "expected: %Mx\n", 250 hgcd_init_tp[hgcd_init_scratch], marker[1]); 251 if (hgcd_tp[-1] != marker[2]) 252 gmp_fprintf (stderr, 253 "before tp: %Mx\n" 254 "expected: %Mx\n", 255 hgcd_tp[-1], marker[2]); 256 if (hgcd_tp[hgcd_scratch] != marker[3]) 257 gmp_fprintf (stderr, 258 "after tp: %Mx\n" 259 "expected: %Mx\n", 260 hgcd_tp[hgcd_scratch], marker[3]); 261 262 abort (); 263 } 264 265 if (!hgcd_appr_valid_p (a, b, res[0], &ref, ref_r0, ref_r1, 266 res[1], &hgcd)) 267 { 268 fprintf (stderr, "ERROR in test %d\n", i); 269 fprintf (stderr, "Invalid results for hgcd and hgcd_ref\n"); 270 fprintf (stderr, "op1="); debug_mp (a, -16); 271 fprintf (stderr, "op2="); debug_mp (b, -16); 272 fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]); 273 fprintf (stderr, "mpn_hgcd_appr: %ld\n", (long) res[1]); 274 abort (); 275 } 276 277 refmpn_free_limbs (hgcd_init_tp - 1); 278 refmpn_free_limbs (hgcd_tp - 1); 279 hgcd_ref_clear (&ref); 280 mpz_clear (ref_r0); 281 mpz_clear (ref_r1); 282 mpz_clear (hgcd_r0); 283 mpz_clear (hgcd_r1); 284 285 return res[0]; 286 } 287 288 static void 289 hgcd_ref_init (struct hgcd_ref *hgcd) 290 { 291 unsigned i; 292 for (i = 0; i<2; i++) 293 { 294 unsigned j; 295 for (j = 0; j<2; j++) 296 mpz_init (hgcd->m[i][j]); 297 } 298 } 299 300 static void 301 hgcd_ref_clear (struct hgcd_ref *hgcd) 302 { 303 unsigned i; 304 for (i = 0; i<2; i++) 305 { 306 unsigned j; 307 for (j = 0; j<2; j++) 308 mpz_clear (hgcd->m[i][j]); 309 } 310 } 311 312 static int 313 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b) 314 { 315 mpz_fdiv_qr (q, r, a, b); 316 if (mpz_size (r) <= s) 317 { 318 mpz_add (r, r, b); 319 mpz_sub_ui (q, q, 1); 320 } 321 322 return (mpz_sgn (q) > 0); 323 } 324 325 static int 326 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b) 327 { 328 mp_size_t n = MAX (mpz_size (a), mpz_size (b)); 329 mp_size_t s = n/2 + 1; 330 mpz_t q; 331 int res; 332 333 if (mpz_size (a) <= s || mpz_size (b) <= s) 334 return 0; 335 336 res = mpz_cmp (a, b); 337 if (res < 0) 338 { 339 mpz_sub (b, b, a); 340 if (mpz_size (b) <= s) 341 return 0; 342 343 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0); 344 mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1); 345 } 346 else if (res > 0) 347 { 348 mpz_sub (a, a, b); 349 if (mpz_size (a) <= s) 350 return 0; 351 352 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1); 353 mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1); 354 } 355 else 356 return 0; 357 358 mpz_init (q); 359 360 for (;;) 361 { 362 ASSERT (mpz_size (a) > s); 363 ASSERT (mpz_size (b) > s); 364 365 if (mpz_cmp (a, b) > 0) 366 { 367 if (!sdiv_qr (q, a, s, a, b)) 368 break; 369 mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]); 370 mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]); 371 } 372 else 373 { 374 if (!sdiv_qr (q, b, s, b, a)) 375 break; 376 mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]); 377 mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]); 378 } 379 } 380 381 mpz_clear (q); 382 383 return 1; 384 } 385 386 static int 387 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize) 388 { 389 mp_srcptr ap = a->_mp_d; 390 mp_size_t asize = a->_mp_size; 391 392 MPN_NORMALIZE (bp, bsize); 393 return asize == bsize && mpn_cmp (ap, bp, asize) == 0; 394 } 395 396 static int 397 hgcd_ref_equal (const struct hgcd_ref *A, const struct hgcd_ref *B) 398 { 399 unsigned i; 400 401 for (i = 0; i<2; i++) 402 { 403 unsigned j; 404 405 for (j = 0; j<2; j++) 406 if (mpz_cmp (A->m[i][j], B->m[i][j]) != 0) 407 return 0; 408 } 409 410 return 1; 411 } 412 413 static int 414 hgcd_appr_valid_p (mpz_t a, mpz_t b, mp_size_t res0, 415 struct hgcd_ref *ref, mpz_t ref_r0, mpz_t ref_r1, 416 mp_size_t res1, struct hgcd_matrix *hgcd) 417 { 418 mp_size_t n = MAX (mpz_size (a), mpz_size (b)); 419 mp_size_t s = n/2 + 1; 420 421 mp_bitcnt_t dbits, abits, margin; 422 mpz_t appr_r0, appr_r1, t, q; 423 struct hgcd_ref appr; 424 425 if (!res0) 426 { 427 if (!res1) 428 return 1; 429 430 fprintf (stderr, "mpn_hgcd_appr returned 1 when no reduction possible.\n"); 431 return 0; 432 } 433 434 /* NOTE: No *_clear calls on error return, since we're going to 435 abort anyway. */ 436 mpz_init (t); 437 mpz_init (q); 438 hgcd_ref_init (&appr); 439 mpz_init (appr_r0); 440 mpz_init (appr_r1); 441 442 if (mpz_size (ref_r0) <= s) 443 { 444 fprintf (stderr, "ref_r0 too small!!!: "); debug_mp (ref_r0, 16); 445 return 0; 446 } 447 if (mpz_size (ref_r1) <= s) 448 { 449 fprintf (stderr, "ref_r1 too small!!!: "); debug_mp (ref_r1, 16); 450 return 0; 451 } 452 453 mpz_sub (t, ref_r0, ref_r1); 454 dbits = mpz_sizeinbase (t, 2); 455 if (dbits > s*GMP_NUMB_BITS) 456 { 457 fprintf (stderr, "ref |r0 - r1| too large!!!: "); debug_mp (t, 16); 458 return 0; 459 } 460 461 if (!res1) 462 { 463 mpz_set (appr_r0, a); 464 mpz_set (appr_r1, b); 465 } 466 else 467 { 468 unsigned i; 469 470 for (i = 0; i<2; i++) 471 { 472 unsigned j; 473 474 for (j = 0; j<2; j++) 475 { 476 mp_size_t mn = hgcd->n; 477 MPN_NORMALIZE (hgcd->p[i][j], mn); 478 mpz_realloc (appr.m[i][j], mn); 479 MPN_COPY (PTR (appr.m[i][j]), hgcd->p[i][j], mn); 480 SIZ (appr.m[i][j]) = mn; 481 } 482 } 483 mpz_mul (appr_r0, appr.m[1][1], a); 484 mpz_mul (t, appr.m[0][1], b); 485 mpz_sub (appr_r0, appr_r0, t); 486 if (mpz_sgn (appr_r0) <= 0 487 || mpz_size (appr_r0) <= s) 488 { 489 fprintf (stderr, "appr_r0 too small: "); debug_mp (appr_r0, 16); 490 return 0; 491 } 492 493 mpz_mul (appr_r1, appr.m[1][0], a); 494 mpz_mul (t, appr.m[0][0], b); 495 mpz_sub (appr_r1, t, appr_r1); 496 if (mpz_sgn (appr_r1) <= 0 497 || mpz_size (appr_r1) <= s) 498 { 499 fprintf (stderr, "appr_r1 too small: "); debug_mp (appr_r1, 16); 500 return 0; 501 } 502 } 503 504 mpz_sub (t, appr_r0, appr_r1); 505 abits = mpz_sizeinbase (t, 2); 506 if (abits < dbits) 507 { 508 fprintf (stderr, "|r0 - r1| too small: "); debug_mp (t, 16); 509 return 0; 510 } 511 512 /* We lose one bit each time we discard the least significant limbs. 513 For the lehmer code, that can happen at most s * (GMP_NUMB_BITS) 514 / (GMP_NUMB_BITS - 1) times. For the dc code, we lose an entire 515 limb (or more?) for each level of recursion. */ 516 517 margin = (n/2+1) * GMP_NUMB_BITS / (GMP_NUMB_BITS - 1); 518 { 519 mp_size_t rn; 520 for (rn = n; ABOVE_THRESHOLD (rn, HGCD_APPR_THRESHOLD); rn = (rn + 1)/2) 521 margin += GMP_NUMB_BITS; 522 } 523 524 if (verbose_flag && abits > dbits) 525 fprintf (stderr, "n = %u: sbits = %u: ref #(r0-r1): %u, appr #(r0-r1): %u excess: %d, margin: %u\n", 526 (unsigned) n, (unsigned) s*GMP_NUMB_BITS, 527 (unsigned) dbits, (unsigned) abits, 528 (int) abits - s * GMP_NUMB_BITS, (unsigned) margin); 529 530 if (abits > s*GMP_NUMB_BITS + margin) 531 { 532 fprintf (stderr, "appr |r0 - r1| much larger than minimal (by %u bits, margin = %u bits)\n", 533 (unsigned) (abits - s*GMP_NUMB_BITS), (unsigned) margin); 534 return 0; 535 } 536 537 while (mpz_cmp (appr_r0, ref_r0) > 0 || mpz_cmp (appr_r1, ref_r1) > 0) 538 { 539 ASSERT (mpz_size (appr_r0) > s); 540 ASSERT (mpz_size (appr_r1) > s); 541 542 if (mpz_cmp (appr_r0, appr_r1) > 0) 543 { 544 if (!sdiv_qr (q, appr_r0, s, appr_r0, appr_r1)) 545 break; 546 mpz_addmul (appr.m[0][1], q, appr.m[0][0]); 547 mpz_addmul (appr.m[1][1], q, appr.m[1][0]); 548 } 549 else 550 { 551 if (!sdiv_qr (q, appr_r1, s, appr_r1, appr_r0)) 552 break; 553 mpz_addmul (appr.m[0][0], q, appr.m[0][1]); 554 mpz_addmul (appr.m[1][0], q, appr.m[1][1]); 555 } 556 } 557 558 if (mpz_cmp (appr_r0, ref_r0) != 0 559 || mpz_cmp (appr_r1, ref_r1) != 0 560 || !hgcd_ref_equal (ref, &appr)) 561 { 562 fprintf (stderr, "appr_r0: "); debug_mp (appr_r0, 16); 563 fprintf (stderr, "ref_r0: "); debug_mp (ref_r0, 16); 564 565 fprintf (stderr, "appr_r1: "); debug_mp (appr_r1, 16); 566 fprintf (stderr, "ref_r1: "); debug_mp (ref_r1, 16); 567 568 return 0; 569 } 570 mpz_clear (t); 571 mpz_clear (q); 572 hgcd_ref_clear (&appr); 573 mpz_clear (appr_r0); 574 mpz_clear (appr_r1); 575 576 return 1; 577 } 578