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