1 /* Test mpn_hgcd_appr. 2 3 Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004, 2011 Free Software 4 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 #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 mp_size_t 161 one_test (mpz_t a, mpz_t b, int i) 162 { 163 struct hgcd_matrix hgcd; 164 struct hgcd_ref ref; 165 166 mpz_t ref_r0; 167 mpz_t ref_r1; 168 mpz_t hgcd_r0; 169 mpz_t hgcd_r1; 170 171 int res[2]; 172 mp_size_t asize; 173 mp_size_t bsize; 174 175 mp_size_t hgcd_init_scratch; 176 mp_size_t hgcd_scratch; 177 178 mp_ptr hgcd_init_tp; 179 mp_ptr hgcd_tp; 180 mp_limb_t marker[4]; 181 182 asize = a->_mp_size; 183 bsize = b->_mp_size; 184 185 ASSERT (asize >= bsize); 186 187 hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize); 188 hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch + 2) + 1; 189 mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp); 190 191 hgcd_scratch = mpn_hgcd_appr_itch (asize); 192 hgcd_tp = refmpn_malloc_limbs (hgcd_scratch + 2) + 1; 193 194 mpn_random (marker, 4); 195 196 hgcd_init_tp[-1] = marker[0]; 197 hgcd_init_tp[hgcd_init_scratch] = marker[1]; 198 hgcd_tp[-1] = marker[2]; 199 hgcd_tp[hgcd_scratch] = marker[3]; 200 201 #if 0 202 fprintf (stderr, 203 "one_test: i = %d asize = %d, bsize = %d\n", 204 i, a->_mp_size, b->_mp_size); 205 206 gmp_fprintf (stderr, 207 "one_test: i = %d\n" 208 " a = %Zx\n" 209 " b = %Zx\n", 210 i, a, b); 211 #endif 212 hgcd_ref_init (&ref); 213 214 mpz_init_set (ref_r0, a); 215 mpz_init_set (ref_r1, b); 216 res[0] = hgcd_ref (&ref, ref_r0, ref_r1); 217 218 mpz_init_set (hgcd_r0, a); 219 mpz_init_set (hgcd_r1, b); 220 if (bsize < asize) 221 { 222 _mpz_realloc (hgcd_r1, asize); 223 MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize); 224 } 225 res[1] = mpn_hgcd_appr (hgcd_r0->_mp_d, 226 hgcd_r1->_mp_d, 227 asize, 228 &hgcd, hgcd_tp); 229 230 if (hgcd_init_tp[-1] != marker[0] 231 || hgcd_init_tp[hgcd_init_scratch] != marker[1] 232 || hgcd_tp[-1] != marker[2] 233 || hgcd_tp[hgcd_scratch] != marker[3]) 234 { 235 fprintf (stderr, "ERROR in test %d\n", i); 236 fprintf (stderr, "scratch space overwritten!\n"); 237 238 if (hgcd_init_tp[-1] != marker[0]) 239 gmp_fprintf (stderr, 240 "before init_tp: %Mx\n" 241 "expected: %Mx\n", 242 hgcd_init_tp[-1], marker[0]); 243 if (hgcd_init_tp[hgcd_init_scratch] != marker[1]) 244 gmp_fprintf (stderr, 245 "after init_tp: %Mx\n" 246 "expected: %Mx\n", 247 hgcd_init_tp[hgcd_init_scratch], marker[1]); 248 if (hgcd_tp[-1] != marker[2]) 249 gmp_fprintf (stderr, 250 "before tp: %Mx\n" 251 "expected: %Mx\n", 252 hgcd_tp[-1], marker[2]); 253 if (hgcd_tp[hgcd_scratch] != marker[3]) 254 gmp_fprintf (stderr, 255 "after tp: %Mx\n" 256 "expected: %Mx\n", 257 hgcd_tp[hgcd_scratch], marker[3]); 258 259 abort (); 260 } 261 262 if (!hgcd_appr_valid_p (a, b, res[0], &ref, ref_r0, ref_r1, 263 res[1], &hgcd)) 264 { 265 fprintf (stderr, "ERROR in test %d\n", i); 266 fprintf (stderr, "Invalid results for hgcd and hgcd_ref\n"); 267 fprintf (stderr, "op1="); debug_mp (a, -16); 268 fprintf (stderr, "op2="); debug_mp (b, -16); 269 fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]); 270 fprintf (stderr, "mpn_hgcd_appr: %ld\n", (long) res[1]); 271 abort (); 272 } 273 274 refmpn_free_limbs (hgcd_init_tp - 1); 275 refmpn_free_limbs (hgcd_tp - 1); 276 hgcd_ref_clear (&ref); 277 mpz_clear (ref_r0); 278 mpz_clear (ref_r1); 279 mpz_clear (hgcd_r0); 280 mpz_clear (hgcd_r1); 281 282 return res[0]; 283 } 284 285 static void 286 hgcd_ref_init (struct hgcd_ref *hgcd) 287 { 288 unsigned i; 289 for (i = 0; i<2; i++) 290 { 291 unsigned j; 292 for (j = 0; j<2; j++) 293 mpz_init (hgcd->m[i][j]); 294 } 295 } 296 297 static void 298 hgcd_ref_clear (struct hgcd_ref *hgcd) 299 { 300 unsigned i; 301 for (i = 0; i<2; i++) 302 { 303 unsigned j; 304 for (j = 0; j<2; j++) 305 mpz_clear (hgcd->m[i][j]); 306 } 307 } 308 309 static int 310 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b) 311 { 312 mpz_fdiv_qr (q, r, a, b); 313 if (mpz_size (r) <= s) 314 { 315 mpz_add (r, r, b); 316 mpz_sub_ui (q, q, 1); 317 } 318 319 return (mpz_sgn (q) > 0); 320 } 321 322 static int 323 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b) 324 { 325 mp_size_t n = MAX (mpz_size (a), mpz_size (b)); 326 mp_size_t s = n/2 + 1; 327 mpz_t q; 328 int res; 329 330 if (mpz_size (a) <= s || mpz_size (b) <= s) 331 return 0; 332 333 res = mpz_cmp (a, b); 334 if (res < 0) 335 { 336 mpz_sub (b, b, a); 337 if (mpz_size (b) <= s) 338 return 0; 339 340 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0); 341 mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1); 342 } 343 else if (res > 0) 344 { 345 mpz_sub (a, a, b); 346 if (mpz_size (a) <= s) 347 return 0; 348 349 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1); 350 mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1); 351 } 352 else 353 return 0; 354 355 mpz_init (q); 356 357 for (;;) 358 { 359 ASSERT (mpz_size (a) > s); 360 ASSERT (mpz_size (b) > s); 361 362 if (mpz_cmp (a, b) > 0) 363 { 364 if (!sdiv_qr (q, a, s, a, b)) 365 break; 366 mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]); 367 mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]); 368 } 369 else 370 { 371 if (!sdiv_qr (q, b, s, b, a)) 372 break; 373 mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]); 374 mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]); 375 } 376 } 377 378 mpz_clear (q); 379 380 return 1; 381 } 382 383 static int 384 hgcd_ref_equal (const struct hgcd_ref *A, const struct hgcd_ref *B) 385 { 386 unsigned i; 387 388 for (i = 0; i<2; i++) 389 { 390 unsigned j; 391 392 for (j = 0; j<2; j++) 393 if (mpz_cmp (A->m[i][j], B->m[i][j]) != 0) 394 return 0; 395 } 396 397 return 1; 398 } 399 400 static int 401 hgcd_appr_valid_p (mpz_t a, mpz_t b, mp_size_t res0, 402 struct hgcd_ref *ref, mpz_t ref_r0, mpz_t ref_r1, 403 mp_size_t res1, struct hgcd_matrix *hgcd) 404 { 405 mp_size_t n = MAX (mpz_size (a), mpz_size (b)); 406 mp_size_t s = n/2 + 1; 407 408 mp_bitcnt_t dbits, abits, margin; 409 mpz_t appr_r0, appr_r1, t, q; 410 struct hgcd_ref appr; 411 412 if (!res0) 413 { 414 if (!res1) 415 return 1; 416 417 fprintf (stderr, "mpn_hgcd_appr returned 1 when no reduction possible.\n"); 418 return 0; 419 } 420 421 /* NOTE: No *_clear calls on error return, since we're going to 422 abort anyway. */ 423 mpz_init (t); 424 mpz_init (q); 425 hgcd_ref_init (&appr); 426 mpz_init (appr_r0); 427 mpz_init (appr_r1); 428 429 if (mpz_size (ref_r0) <= s) 430 { 431 fprintf (stderr, "ref_r0 too small!!!: "); debug_mp (ref_r0, 16); 432 return 0; 433 } 434 if (mpz_size (ref_r1) <= s) 435 { 436 fprintf (stderr, "ref_r1 too small!!!: "); debug_mp (ref_r1, 16); 437 return 0; 438 } 439 440 mpz_sub (t, ref_r0, ref_r1); 441 dbits = mpz_sizeinbase (t, 2); 442 if (dbits > s*GMP_NUMB_BITS) 443 { 444 fprintf (stderr, "ref |r0 - r1| too large!!!: "); debug_mp (t, 16); 445 return 0; 446 } 447 448 if (!res1) 449 { 450 mpz_set (appr_r0, a); 451 mpz_set (appr_r1, b); 452 } 453 else 454 { 455 unsigned i; 456 457 for (i = 0; i<2; i++) 458 { 459 unsigned j; 460 461 for (j = 0; j<2; j++) 462 { 463 mp_size_t mn = hgcd->n; 464 MPN_NORMALIZE (hgcd->p[i][j], mn); 465 mpz_realloc (appr.m[i][j], mn); 466 MPN_COPY (PTR (appr.m[i][j]), hgcd->p[i][j], mn); 467 SIZ (appr.m[i][j]) = mn; 468 } 469 } 470 mpz_mul (appr_r0, appr.m[1][1], a); 471 mpz_mul (t, appr.m[0][1], b); 472 mpz_sub (appr_r0, appr_r0, t); 473 if (mpz_sgn (appr_r0) <= 0 474 || mpz_size (appr_r0) <= s) 475 { 476 fprintf (stderr, "appr_r0 too small: "); debug_mp (appr_r0, 16); 477 return 0; 478 } 479 480 mpz_mul (appr_r1, appr.m[1][0], a); 481 mpz_mul (t, appr.m[0][0], b); 482 mpz_sub (appr_r1, t, appr_r1); 483 if (mpz_sgn (appr_r1) <= 0 484 || mpz_size (appr_r1) <= s) 485 { 486 fprintf (stderr, "appr_r1 too small: "); debug_mp (appr_r1, 16); 487 return 0; 488 } 489 } 490 491 mpz_sub (t, appr_r0, appr_r1); 492 abits = mpz_sizeinbase (t, 2); 493 if (abits < dbits) 494 { 495 fprintf (stderr, "|r0 - r1| too small: "); debug_mp (t, 16); 496 return 0; 497 } 498 499 /* We lose one bit each time we discard the least significant limbs. 500 For the lehmer code, that can happen at most s * (GMP_NUMB_BITS) 501 / (GMP_NUMB_BITS - 1) times. For the dc code, we lose an entire 502 limb (or more?) for each level of recursion. */ 503 504 margin = (n/2+1) * GMP_NUMB_BITS / (GMP_NUMB_BITS - 1); 505 { 506 mp_size_t rn; 507 for (rn = n; ABOVE_THRESHOLD (rn, HGCD_APPR_THRESHOLD); rn = (rn + 1)/2) 508 margin += GMP_NUMB_BITS; 509 } 510 511 if (verbose_flag && abits > dbits) 512 fprintf (stderr, "n = %u: sbits = %u: ref #(r0-r1): %u, appr #(r0-r1): %u excess: %d, margin: %u\n", 513 (unsigned) n, (unsigned) s*GMP_NUMB_BITS, 514 (unsigned) dbits, (unsigned) abits, 515 (int) (abits - s * GMP_NUMB_BITS), (unsigned) margin); 516 517 if (abits > s*GMP_NUMB_BITS + margin) 518 { 519 fprintf (stderr, "appr |r0 - r1| much larger than minimal (by %u bits, margin = %u bits)\n", 520 (unsigned) (abits - s*GMP_NUMB_BITS), (unsigned) margin); 521 return 0; 522 } 523 524 while (mpz_cmp (appr_r0, ref_r0) > 0 || mpz_cmp (appr_r1, ref_r1) > 0) 525 { 526 ASSERT (mpz_size (appr_r0) > s); 527 ASSERT (mpz_size (appr_r1) > s); 528 529 if (mpz_cmp (appr_r0, appr_r1) > 0) 530 { 531 if (!sdiv_qr (q, appr_r0, s, appr_r0, appr_r1)) 532 break; 533 mpz_addmul (appr.m[0][1], q, appr.m[0][0]); 534 mpz_addmul (appr.m[1][1], q, appr.m[1][0]); 535 } 536 else 537 { 538 if (!sdiv_qr (q, appr_r1, s, appr_r1, appr_r0)) 539 break; 540 mpz_addmul (appr.m[0][0], q, appr.m[0][1]); 541 mpz_addmul (appr.m[1][0], q, appr.m[1][1]); 542 } 543 } 544 545 if (mpz_cmp (appr_r0, ref_r0) != 0 546 || mpz_cmp (appr_r1, ref_r1) != 0 547 || !hgcd_ref_equal (ref, &appr)) 548 { 549 fprintf (stderr, "appr_r0: "); debug_mp (appr_r0, 16); 550 fprintf (stderr, "ref_r0: "); debug_mp (ref_r0, 16); 551 552 fprintf (stderr, "appr_r1: "); debug_mp (appr_r1, 16); 553 fprintf (stderr, "ref_r1: "); debug_mp (ref_r1, 16); 554 555 return 0; 556 } 557 mpz_clear (t); 558 mpz_clear (q); 559 hgcd_ref_clear (&appr); 560 mpz_clear (appr_r0); 561 mpz_clear (appr_r1); 562 563 return 1; 564 } 565