1 /* Test mpn_hgcd. 2 3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004 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 24 #include "gmp.h" 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 /* Fixed values, for regression testing of mpn_hgcd. */ 34 struct value { int res; const char *a; const char *b; }; 35 static const struct value hgcd_values[] = { 36 #if GMP_NUMB_BITS == 32 37 { 5, 38 "0x1bddff867272a9296ac493c251d7f46f09a5591fe", 39 "0xb55930a2a68a916450a7de006031068c5ddb0e5c" }, 40 { 4, 41 "0x2f0ece5b1ee9c15e132a01d55768dc13", 42 "0x1c6f4fd9873cdb24466e6d03e1cc66e7" }, 43 { 3, "0x7FFFFC003FFFFFFFFFC5", "0x3FFFFE001FFFFFFFFFE3"}, 44 #endif 45 { -1, NULL, NULL } 46 }; 47 48 struct hgcd_ref 49 { 50 mpz_t m[2][2]; 51 }; 52 53 static void hgcd_ref_init (struct hgcd_ref *); 54 static void hgcd_ref_clear (struct hgcd_ref *); 55 static int hgcd_ref (struct hgcd_ref *, mpz_t, mpz_t); 56 static int hgcd_ref_equal (const struct hgcd_matrix *, const struct hgcd_ref *); 57 58 int 59 main (int argc, char **argv) 60 { 61 mpz_t op1, op2, temp1, temp2; 62 int i, j, chain_len; 63 gmp_randstate_ptr rands; 64 mpz_t bs; 65 unsigned long size_range; 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; hgcd_values[i].res >= 0; i++) 77 { 78 mp_size_t res; 79 80 mpz_set_str (op1, hgcd_values[i].a, 0); 81 mpz_set_str (op2, hgcd_values[i].b, 0); 82 83 res = one_test (op1, op2, -1-i); 84 if (res != hgcd_values[i].res) 85 { 86 fprintf (stderr, "ERROR in test %d\n", -1-i); 87 fprintf (stderr, "Bad return code from hgcd\n"); 88 fprintf (stderr, "op1="); debug_mp (op1, -16); 89 fprintf (stderr, "op2="); debug_mp (op2, -16); 90 fprintf (stderr, "expected: %d\n", hgcd_values[i].res); 91 fprintf (stderr, "hgcd: %d\n", (int) res); 92 abort (); 93 } 94 } 95 96 for (i = 0; i < 15; i++) 97 { 98 /* Generate plain operands with unknown gcd. These types of operands 99 have proven to trigger certain bugs in development versions of the 100 gcd code. */ 101 102 mpz_urandomb (bs, rands, 32); 103 size_range = mpz_get_ui (bs) % 13 + 2; 104 105 mpz_urandomb (bs, rands, size_range); 106 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); 107 mpz_urandomb (bs, rands, size_range); 108 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); 109 110 if (mpz_cmp (op1, op2) < 0) 111 mpz_swap (op1, op2); 112 113 if (mpz_size (op1) > 0) 114 one_test (op1, op2, i); 115 116 /* Generate a division chain backwards, allowing otherwise 117 unlikely huge quotients. */ 118 119 mpz_set_ui (op1, 0); 120 mpz_urandomb (bs, rands, 32); 121 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1); 122 mpz_rrandomb (op2, rands, mpz_get_ui (bs)); 123 mpz_add_ui (op2, op2, 1); 124 125 #if 0 126 chain_len = 1000000; 127 #else 128 mpz_urandomb (bs, rands, 32); 129 chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256); 130 #endif 131 132 for (j = 0; j < chain_len; j++) 133 { 134 mpz_urandomb (bs, rands, 32); 135 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); 136 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); 137 mpz_add_ui (temp2, temp2, 1); 138 mpz_mul (temp1, op2, temp2); 139 mpz_add (op1, op1, temp1); 140 141 /* Don't generate overly huge operands. */ 142 if (SIZ (op1) > 3 * GCD_DC_THRESHOLD) 143 break; 144 145 mpz_urandomb (bs, rands, 32); 146 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); 147 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); 148 mpz_add_ui (temp2, temp2, 1); 149 mpz_mul (temp1, op1, temp2); 150 mpz_add (op2, op2, temp1); 151 152 /* Don't generate overly huge operands. */ 153 if (SIZ (op2) > 3 * GCD_DC_THRESHOLD) 154 break; 155 } 156 if (mpz_cmp (op1, op2) < 0) 157 mpz_swap (op1, op2); 158 159 if (mpz_size (op1) > 0) 160 one_test (op1, op2, i); 161 } 162 163 mpz_clear (bs); 164 mpz_clear (op1); 165 mpz_clear (op2); 166 mpz_clear (temp1); 167 mpz_clear (temp2); 168 169 tests_end (); 170 exit (0); 171 } 172 173 static void 174 debug_mp (mpz_t x, int base) 175 { 176 mpz_out_str (stderr, base, x); fputc ('\n', stderr); 177 } 178 179 static int 180 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize); 181 182 static mp_size_t 183 one_test (mpz_t a, mpz_t b, int i) 184 { 185 struct hgcd_matrix hgcd; 186 struct hgcd_ref ref; 187 188 mpz_t ref_r0; 189 mpz_t ref_r1; 190 mpz_t hgcd_r0; 191 mpz_t hgcd_r1; 192 193 mp_size_t res[2]; 194 mp_size_t asize; 195 mp_size_t bsize; 196 197 mp_size_t hgcd_init_scratch; 198 mp_size_t hgcd_scratch; 199 200 mp_ptr hgcd_init_tp; 201 mp_ptr hgcd_tp; 202 203 asize = a->_mp_size; 204 bsize = b->_mp_size; 205 206 ASSERT (asize >= bsize); 207 208 hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize); 209 hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch); 210 mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp); 211 212 hgcd_scratch = mpn_hgcd_itch (asize); 213 hgcd_tp = refmpn_malloc_limbs (hgcd_scratch); 214 215 #if 0 216 fprintf (stderr, 217 "one_test: i = %d asize = %d, bsize = %d\n", 218 i, a->_mp_size, b->_mp_size); 219 220 gmp_fprintf (stderr, 221 "one_test: i = %d\n" 222 " a = %Zx\n" 223 " b = %Zx\n", 224 i, a, b); 225 #endif 226 hgcd_ref_init (&ref); 227 228 mpz_init_set (ref_r0, a); 229 mpz_init_set (ref_r1, b); 230 res[0] = hgcd_ref (&ref, ref_r0, ref_r1); 231 232 mpz_init_set (hgcd_r0, a); 233 mpz_init_set (hgcd_r1, b); 234 if (bsize < asize) 235 { 236 _mpz_realloc (hgcd_r1, asize); 237 MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize); 238 } 239 res[1] = mpn_hgcd (hgcd_r0->_mp_d, 240 hgcd_r1->_mp_d, 241 asize, 242 &hgcd, hgcd_tp); 243 244 if (res[0] != res[1]) 245 { 246 fprintf (stderr, "ERROR in test %d\n", i); 247 fprintf (stderr, "Different return value from hgcd and hgcd_ref\n"); 248 fprintf (stderr, "op1="); debug_mp (a, -16); 249 fprintf (stderr, "op2="); debug_mp (b, -16); 250 fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]); 251 fprintf (stderr, "mpn_hgcd: %ld\n", (long) res[1]); 252 abort (); 253 } 254 if (res[0] > 0) 255 { 256 if (!hgcd_ref_equal (&hgcd, &ref) 257 || !mpz_mpn_equal (ref_r0, hgcd_r0->_mp_d, res[1]) 258 || !mpz_mpn_equal (ref_r1, hgcd_r1->_mp_d, res[1])) 259 { 260 fprintf (stderr, "ERROR in test %d\n", i); 261 fprintf (stderr, "mpn_hgcd and hgcd_ref returned different values\n"); 262 fprintf (stderr, "op1="); debug_mp (a, -16); 263 fprintf (stderr, "op2="); debug_mp (b, -16); 264 abort (); 265 } 266 } 267 268 refmpn_free_limbs (hgcd_init_tp); 269 refmpn_free_limbs (hgcd_tp); 270 hgcd_ref_clear (&ref); 271 mpz_clear (ref_r0); 272 mpz_clear (ref_r1); 273 mpz_clear (hgcd_r0); 274 mpz_clear (hgcd_r1); 275 276 return res[0]; 277 } 278 279 static void 280 hgcd_ref_init (struct hgcd_ref *hgcd) 281 { 282 unsigned i; 283 for (i = 0; i<2; i++) 284 { 285 unsigned j; 286 for (j = 0; j<2; j++) 287 mpz_init (hgcd->m[i][j]); 288 } 289 } 290 291 static void 292 hgcd_ref_clear (struct hgcd_ref *hgcd) 293 { 294 unsigned i; 295 for (i = 0; i<2; i++) 296 { 297 unsigned j; 298 for (j = 0; j<2; j++) 299 mpz_clear (hgcd->m[i][j]); 300 } 301 } 302 303 304 static int 305 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b) 306 { 307 mpz_fdiv_qr (q, r, a, b); 308 if (mpz_size (r) <= s) 309 { 310 mpz_add (r, r, b); 311 mpz_sub_ui (q, q, 1); 312 } 313 314 return (mpz_sgn (q) > 0); 315 } 316 317 static int 318 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b) 319 { 320 mp_size_t n = MAX (mpz_size (a), mpz_size (b)); 321 mp_size_t s = n/2 + 1; 322 mp_size_t asize; 323 mp_size_t bsize; 324 mpz_t q; 325 int res; 326 327 if (mpz_size (a) <= s || mpz_size (b) <= s) 328 return 0; 329 330 res = mpz_cmp (a, b); 331 if (res < 0) 332 { 333 mpz_sub (b, b, a); 334 if (mpz_size (b) <= s) 335 return 0; 336 337 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0); 338 mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1); 339 } 340 else if (res > 0) 341 { 342 mpz_sub (a, a, b); 343 if (mpz_size (a) <= s) 344 return 0; 345 346 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1); 347 mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1); 348 } 349 else 350 return 0; 351 352 mpz_init (q); 353 354 for (;;) 355 { 356 ASSERT (mpz_size (a) > s); 357 ASSERT (mpz_size (b) > s); 358 359 if (mpz_cmp (a, b) > 0) 360 { 361 if (!sdiv_qr (q, a, s, a, b)) 362 break; 363 mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]); 364 mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]); 365 } 366 else 367 { 368 if (!sdiv_qr (q, b, s, b, a)) 369 break; 370 mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]); 371 mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]); 372 } 373 } 374 375 mpz_clear (q); 376 377 asize = mpz_size (a); 378 bsize = mpz_size (b); 379 return MAX (asize, bsize); 380 } 381 382 static int 383 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize) 384 { 385 mp_srcptr ap = a->_mp_d; 386 mp_size_t asize = a->_mp_size; 387 388 MPN_NORMALIZE (bp, bsize); 389 return asize == bsize && mpn_cmp (ap, bp, asize) == 0; 390 } 391 392 static int 393 hgcd_ref_equal (const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref) 394 { 395 unsigned i; 396 397 for (i = 0; i<2; i++) 398 { 399 unsigned j; 400 401 for (j = 0; j<2; j++) 402 if (!mpz_mpn_equal (ref->m[i][j], hgcd->p[i][j], hgcd->n)) 403 return 0; 404 } 405 406 return 1; 407 } 408