1 /* Copyright 2006, 2007, 2009, 2010, 2013 Free Software Foundation, Inc. 2 3 This file is part of the GNU MP Library test suite. 4 5 The GNU MP Library test suite is free software; you can redistribute it 6 and/or modify it under the terms of the GNU General Public License as 7 published by the Free Software Foundation; either version 3 of the License, 8 or (at your option) any later version. 9 10 The GNU MP Library test suite is distributed in the hope that it will be 11 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 13 Public License for more details. 14 15 You should have received a copy of the GNU General Public License along with 16 the GNU MP Library test suite. If not, see http://www.gnu.org/licenses/. */ 17 18 19 #include <stdlib.h> /* for strtol */ 20 #include <stdio.h> /* for printf */ 21 22 #include "gmp.h" 23 #include "gmp-impl.h" 24 #include "longlong.h" 25 #include "tests/tests.h" 26 27 28 static void 29 dumpy (mp_srcptr p, mp_size_t n) 30 { 31 mp_size_t i; 32 if (n > 20) 33 { 34 for (i = n - 1; i >= n - 4; i--) 35 { 36 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 37 printf (" "); 38 } 39 printf ("... "); 40 for (i = 3; i >= 0; i--) 41 { 42 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 43 printf (" " + (i == 0)); 44 } 45 } 46 else 47 { 48 for (i = n - 1; i >= 0; i--) 49 { 50 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 51 printf (" " + (i == 0)); 52 } 53 } 54 puts (""); 55 } 56 57 static signed long test; 58 59 static void 60 check_one (mp_ptr qp, mp_srcptr rp, 61 mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, 62 const char *fname, mp_limb_t q_allowed_err) 63 { 64 mp_size_t qn = nn - dn + 1; 65 mp_ptr tp; 66 const char *msg; 67 const char *tvalue; 68 mp_limb_t i; 69 TMP_DECL; 70 TMP_MARK; 71 72 tp = TMP_ALLOC_LIMBS (nn + 1); 73 if (dn >= qn) 74 refmpn_mul (tp, dp, dn, qp, qn); 75 else 76 refmpn_mul (tp, qp, qn, dp, dn); 77 78 for (i = 0; i < q_allowed_err && (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0); i++) 79 ASSERT_NOCARRY (refmpn_sub (tp, tp, nn+1, dp, dn)); 80 81 if (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0) 82 { 83 msg = "q too large"; 84 tvalue = "Q*D"; 85 error: 86 printf ("\r*******************************************************************************\n"); 87 printf ("%s failed test %ld: %s\n", fname, test, msg); 88 printf ("N= "); dumpy (np, nn); 89 printf ("D= "); dumpy (dp, dn); 90 printf ("Q= "); dumpy (qp, qn); 91 if (rp) 92 { printf ("R= "); dumpy (rp, dn); } 93 printf ("%5s=", tvalue); dumpy (tp, nn+1); 94 printf ("nn = %ld, dn = %ld, qn = %ld\n", nn, dn, qn); 95 abort (); 96 } 97 98 ASSERT_NOCARRY (refmpn_sub_n (tp, np, tp, nn)); 99 tvalue = "N-Q*D"; 100 if (!mpn_zero_p (tp + dn, nn - dn) || mpn_cmp (tp, dp, dn) >= 0) 101 { 102 msg = "q too small"; 103 goto error; 104 } 105 106 if (rp && mpn_cmp (rp, tp, dn) != 0) 107 { 108 msg = "r incorrect"; 109 goto error; 110 } 111 112 TMP_FREE; 113 } 114 115 116 /* These are *bit* sizes. */ 117 #ifndef SIZE_LOG 118 #define SIZE_LOG 17 119 #endif 120 #define MAX_DN (1L << SIZE_LOG) 121 #define MAX_NN (1L << (SIZE_LOG + 1)) 122 123 #define COUNT 200 124 125 mp_limb_t 126 random_word (gmp_randstate_ptr rs) 127 { 128 mpz_t x; 129 mp_limb_t r; 130 TMP_DECL; 131 TMP_MARK; 132 133 MPZ_TMP_INIT (x, 2); 134 mpz_urandomb (x, rs, 32); 135 r = mpz_get_ui (x); 136 TMP_FREE; 137 return r; 138 } 139 140 int 141 main (int argc, char **argv) 142 { 143 gmp_randstate_ptr rands; 144 unsigned long maxnbits, maxdbits, nbits, dbits; 145 mpz_t n, d, q, r, tz, junk; 146 mp_size_t maxnn, maxdn, nn, dn, clearn, i; 147 mp_ptr np, dup, dnp, qp, rp, junkp; 148 mp_limb_t t; 149 gmp_pi1_t dinv; 150 long count = COUNT; 151 mp_ptr scratch; 152 mp_limb_t ran; 153 mp_size_t alloc, itch; 154 mp_limb_t rran0, rran1, qran0, qran1; 155 TMP_DECL; 156 157 if (argc > 1) 158 { 159 char *end; 160 count = strtol (argv[1], &end, 0); 161 if (*end || count <= 0) 162 { 163 fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 164 return 1; 165 } 166 } 167 168 maxdbits = MAX_DN; 169 maxnbits = MAX_NN; 170 171 tests_start (); 172 rands = RANDS; 173 174 mpz_init (n); 175 mpz_init (d); 176 mpz_init (q); 177 mpz_init (r); 178 mpz_init (tz); 179 mpz_init (junk); 180 181 maxnn = maxnbits / GMP_NUMB_BITS + 1; 182 maxdn = maxdbits / GMP_NUMB_BITS + 1; 183 184 TMP_MARK; 185 186 qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 187 rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 188 dnp = TMP_ALLOC_LIMBS (maxdn); 189 190 alloc = 1; 191 scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc); 192 193 for (test = -300; test < count; test++) 194 { 195 nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS; 196 197 if (test < 0) 198 dbits = (test + 300) % (nbits - 1) + 1; 199 else 200 dbits = random_word (rands) % (nbits - 1) % maxdbits + 1; 201 202 #if RAND_UNIFORM 203 #define RANDFUNC mpz_urandomb 204 #else 205 #define RANDFUNC mpz_rrandomb 206 #endif 207 208 do 209 RANDFUNC (d, rands, dbits); 210 while (mpz_sgn (d) == 0); 211 dn = SIZ (d); 212 dup = PTR (d); 213 MPN_COPY (dnp, dup, dn); 214 dnp[dn - 1] |= GMP_NUMB_HIGHBIT; 215 216 if (test % 2 == 0) 217 { 218 RANDFUNC (n, rands, nbits); 219 nn = SIZ (n); 220 ASSERT_ALWAYS (nn >= dn); 221 } 222 else 223 { 224 do 225 { 226 RANDFUNC (q, rands, random_word (rands) % (nbits - dbits + 1)); 227 RANDFUNC (r, rands, random_word (rands) % mpz_sizeinbase (d, 2)); 228 mpz_mul (n, q, d); 229 mpz_add (n, n, r); 230 nn = SIZ (n); 231 } 232 while (nn > maxnn || nn < dn); 233 } 234 235 ASSERT_ALWAYS (nn <= maxnn); 236 ASSERT_ALWAYS (dn <= maxdn); 237 238 mpz_urandomb (junk, rands, nbits); 239 junkp = PTR (junk); 240 241 np = PTR (n); 242 243 mpz_urandomb (tz, rands, 32); 244 t = mpz_get_ui (tz); 245 246 if (t % 17 == 0) 247 { 248 dnp[dn - 1] = GMP_NUMB_MAX; 249 dup[dn - 1] = GMP_NUMB_MAX; 250 } 251 252 switch ((int) t % 16) 253 { 254 case 0: 255 clearn = random_word (rands) % nn; 256 for (i = clearn; i < nn; i++) 257 np[i] = 0; 258 break; 259 case 1: 260 mpn_sub_1 (np + nn - dn, dnp, dn, random_word (rands)); 261 break; 262 case 2: 263 mpn_add_1 (np + nn - dn, dnp, dn, random_word (rands)); 264 break; 265 } 266 267 if (dn >= 2) 268 invert_pi1 (dinv, dnp[dn - 1], dnp[dn - 2]); 269 270 rran0 = random_word (rands); 271 rran1 = random_word (rands); 272 qran0 = random_word (rands); 273 qran1 = random_word (rands); 274 275 qp[-1] = qran0; 276 qp[nn - dn + 1] = qran1; 277 rp[-1] = rran0; 278 279 ran = random_word (rands); 280 281 if ((double) (nn - dn) * dn < 1e5) 282 { 283 /* Test mpn_sbpi1_div_qr */ 284 if (dn > 2) 285 { 286 MPN_COPY (rp, np, nn); 287 if (nn > dn) 288 MPN_COPY (qp, junkp, nn - dn); 289 qp[nn - dn] = mpn_sbpi1_div_qr (qp, rp, nn, dnp, dn, dinv.inv32); 290 check_one (qp, rp, np, nn, dnp, dn, "mpn_sbpi1_div_qr", 0); 291 } 292 293 /* Test mpn_sbpi1_divappr_q */ 294 if (dn > 2) 295 { 296 MPN_COPY (rp, np, nn); 297 if (nn > dn) 298 MPN_COPY (qp, junkp, nn - dn); 299 qp[nn - dn] = mpn_sbpi1_divappr_q (qp, rp, nn, dnp, dn, dinv.inv32); 300 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_divappr_q", 1); 301 } 302 303 /* Test mpn_sbpi1_div_q */ 304 if (dn > 2) 305 { 306 MPN_COPY (rp, np, nn); 307 if (nn > dn) 308 MPN_COPY (qp, junkp, nn - dn); 309 qp[nn - dn] = mpn_sbpi1_div_q (qp, rp, nn, dnp, dn, dinv.inv32); 310 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_div_q", 0); 311 } 312 313 /* Test mpn_sb_div_qr_sec */ 314 itch = 3 * nn + 4; 315 if (itch + 1 > alloc) 316 { 317 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 318 alloc = itch + 1; 319 } 320 scratch[itch] = ran; 321 MPN_COPY (rp, np, nn); 322 if (nn >= dn) 323 MPN_COPY (qp, junkp, nn - dn + 1); 324 mpn_sb_div_qr_sec (qp, rp, nn, dup, dn, scratch); 325 ASSERT_ALWAYS (ran == scratch[itch]); 326 check_one (qp, rp, np, nn, dup, dn, "mpn_sb_div_qr_sec", 0); 327 328 /* Test mpn_sb_div_r_sec */ 329 itch = nn + 2 * dn + 2; 330 if (itch + 1 > alloc) 331 { 332 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 333 alloc = itch + 1; 334 } 335 scratch[itch] = ran; 336 MPN_COPY (rp, np, nn); 337 mpn_sb_div_r_sec (rp, nn, dup, dn, scratch); 338 ASSERT_ALWAYS (ran == scratch[itch]); 339 /* Note: Since check_one cannot cope with random-only functions, we 340 pass qp[] from the previous function, mpn_sb_div_qr_sec. */ 341 check_one (qp, rp, np, nn, dup, dn, "mpn_sb_div_r_sec", 0); 342 } 343 344 /* Test mpn_dcpi1_div_qr */ 345 if (dn >= 6 && nn - dn >= 3) 346 { 347 MPN_COPY (rp, np, nn); 348 if (nn > dn) 349 MPN_COPY (qp, junkp, nn - dn); 350 qp[nn - dn] = mpn_dcpi1_div_qr (qp, rp, nn, dnp, dn, &dinv); 351 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 352 ASSERT_ALWAYS (rp[-1] == rran0); 353 check_one (qp, rp, np, nn, dnp, dn, "mpn_dcpi1_div_qr", 0); 354 } 355 356 /* Test mpn_dcpi1_divappr_q */ 357 if (dn >= 6 && nn - dn >= 3) 358 { 359 MPN_COPY (rp, np, nn); 360 if (nn > dn) 361 MPN_COPY (qp, junkp, nn - dn); 362 qp[nn - dn] = mpn_dcpi1_divappr_q (qp, rp, nn, dnp, dn, &dinv); 363 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 364 ASSERT_ALWAYS (rp[-1] == rran0); 365 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_divappr_q", 1); 366 } 367 368 /* Test mpn_dcpi1_div_q */ 369 if (dn >= 6 && nn - dn >= 3) 370 { 371 MPN_COPY (rp, np, nn); 372 if (nn > dn) 373 MPN_COPY (qp, junkp, nn - dn); 374 qp[nn - dn] = mpn_dcpi1_div_q (qp, rp, nn, dnp, dn, &dinv); 375 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 376 ASSERT_ALWAYS (rp[-1] == rran0); 377 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_div_q", 0); 378 } 379 380 /* Test mpn_mu_div_qr */ 381 if (nn - dn > 2 && dn >= 2) 382 { 383 itch = mpn_mu_div_qr_itch (nn, dn, 0); 384 if (itch + 1 > alloc) 385 { 386 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 387 alloc = itch + 1; 388 } 389 scratch[itch] = ran; 390 MPN_COPY (qp, junkp, nn - dn); 391 MPN_ZERO (rp, dn); 392 rp[dn] = rran1; 393 qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dnp, dn, scratch); 394 ASSERT_ALWAYS (ran == scratch[itch]); 395 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 396 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1); 397 check_one (qp, rp, np, nn, dnp, dn, "mpn_mu_div_qr", 0); 398 } 399 400 /* Test mpn_mu_divappr_q */ 401 if (nn - dn > 2 && dn >= 2) 402 { 403 itch = mpn_mu_divappr_q_itch (nn, dn, 0); 404 if (itch + 1 > alloc) 405 { 406 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 407 alloc = itch + 1; 408 } 409 scratch[itch] = ran; 410 MPN_COPY (qp, junkp, nn - dn); 411 qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dnp, dn, scratch); 412 ASSERT_ALWAYS (ran == scratch[itch]); 413 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 414 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_divappr_q", 4); 415 } 416 417 /* Test mpn_mu_div_q */ 418 if (nn - dn > 2 && dn >= 2) 419 { 420 itch = mpn_mu_div_q_itch (nn, dn, 0); 421 if (itch + 1> alloc) 422 { 423 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 424 alloc = itch + 1; 425 } 426 scratch[itch] = ran; 427 MPN_COPY (qp, junkp, nn - dn); 428 qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dnp, dn, scratch); 429 ASSERT_ALWAYS (ran == scratch[itch]); 430 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 431 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_div_q", 0); 432 } 433 434 if (1) 435 { 436 itch = nn + 1; 437 if (itch + 1> alloc) 438 { 439 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 440 alloc = itch + 1; 441 } 442 scratch[itch] = ran; 443 mpn_div_q (qp, np, nn, dup, dn, scratch); 444 ASSERT_ALWAYS (ran == scratch[itch]); 445 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 446 check_one (qp, NULL, np, nn, dup, dn, "mpn_div_q", 0); 447 } 448 449 if (dn >= 2 && nn >= 2) 450 { 451 mp_limb_t qh; 452 453 /* mpn_divrem_2 */ 454 MPN_COPY (rp, np, nn); 455 qp[nn - 2] = qp[nn-1] = qran1; 456 457 qh = mpn_divrem_2 (qp, 0, rp, nn, dnp + dn - 2); 458 ASSERT_ALWAYS (qp[nn - 2] == qran1); 459 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1); 460 qp[nn - 2] = qh; 461 check_one (qp, rp, np, nn, dnp + dn - 2, 2, "mpn_divrem_2", 0); 462 463 /* Missing: divrem_2 with fraction limbs. */ 464 465 /* mpn_div_qr_2 */ 466 qp[nn - 2] = qran1; 467 468 qh = mpn_div_qr_2 (qp, rp, np, nn, dup + dn - 2); 469 ASSERT_ALWAYS (qp[nn - 2] == qran1); 470 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1); 471 qp[nn - 2] = qh; 472 check_one (qp, rp, np, nn, dup + dn - 2, 2, "mpn_div_qr_2", 0); 473 } 474 } 475 476 __GMP_FREE_FUNC_LIMBS (scratch, alloc); 477 478 TMP_FREE; 479 480 mpz_clear (n); 481 mpz_clear (d); 482 mpz_clear (q); 483 mpz_clear (r); 484 mpz_clear (tz); 485 mpz_clear (junk); 486 487 tests_end (); 488 return 0; 489 } 490