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