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