t-hgcd_appr.c revision 1.1.1.3
1/* Test mpn_hgcd_appr. 2 3Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004, 2011 Free Software 4Foundation, Inc. 5 6This file is part of the GNU MP Library test suite. 7 8The GNU MP Library test suite is free software; you can redistribute it 9and/or modify it under the terms of the GNU General Public License as 10published by the Free Software Foundation; either version 3 of the License, 11or (at your option) any later version. 12 13The GNU MP Library test suite is distributed in the hope that it will be 14useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 16Public License for more details. 17 18You should have received a copy of the GNU General Public License along with 19the 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 28static mp_size_t one_test (mpz_t, mpz_t, int); 29static void debug_mp (mpz_t, int); 30 31#define MIN_OPERAND_SIZE 2 32 33struct hgcd_ref 34{ 35 mpz_t m[2][2]; 36}; 37 38static void hgcd_ref_init (struct hgcd_ref *hgcd); 39static void hgcd_ref_clear (struct hgcd_ref *hgcd); 40static int hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b); 41static int hgcd_ref_equal (const struct hgcd_ref *, const struct hgcd_ref *); 42static 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 45static int verbose_flag = 0; 46 47int 48main (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 153static void 154debug_mp (mpz_t x, int base) 155{ 156 mpz_out_str (stderr, base, x); fputc ('\n', stderr); 157} 158 159static mp_size_t 160one_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 284static void 285hgcd_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 296static void 297hgcd_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 308static int 309sdiv_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 321static int 322hgcd_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 382static int 383hgcd_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 399static int 400hgcd_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