github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/tests/mpn/t-div.c (about) 1 /* Copyright 2006, 2007, 2009, 2010, 2013-2015 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 https://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 (!(nn == dn || 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 int 126 main (int argc, char **argv) 127 { 128 gmp_randstate_ptr rands; 129 unsigned long maxnbits, maxdbits, nbits, dbits; 130 mpz_t n, d, q, r, tz, junk; 131 mp_size_t maxnn, maxdn, nn, dn, clearn, i; 132 mp_ptr np, dup, dnp, qp, rp, junkp; 133 mp_limb_t t; 134 gmp_pi1_t dinv; 135 long count = COUNT; 136 mp_ptr scratch; 137 mp_limb_t ran; 138 mp_size_t alloc, itch; 139 mp_limb_t rran0, rran1, qran0, qran1; 140 TMP_DECL; 141 142 if (argc > 1) 143 { 144 char *end; 145 count = strtol (argv[1], &end, 0); 146 if (*end || count <= 0) 147 { 148 fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 149 return 1; 150 } 151 } 152 153 maxdbits = MAX_DN; 154 maxnbits = MAX_NN; 155 156 tests_start (); 157 rands = RANDS; 158 159 mpz_init (n); 160 mpz_init (d); 161 mpz_init (q); 162 mpz_init (r); 163 mpz_init (tz); 164 mpz_init (junk); 165 166 maxnn = maxnbits / GMP_NUMB_BITS + 1; 167 maxdn = maxdbits / GMP_NUMB_BITS + 1; 168 169 TMP_MARK; 170 171 qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 172 rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 173 dnp = TMP_ALLOC_LIMBS (maxdn); 174 175 alloc = 1; 176 scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc); 177 178 for (test = -300; test < count; test++) 179 { 180 nbits = urandom () % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS; 181 182 if (test < 0) 183 dbits = (test + 300) % (nbits - 1) + 1; 184 else 185 dbits = urandom () % (nbits - 1) % maxdbits + 1; 186 187 #if RAND_UNIFORM 188 #define RANDFUNC mpz_urandomb 189 #else 190 #define RANDFUNC mpz_rrandomb 191 #endif 192 193 do 194 RANDFUNC (d, rands, dbits); 195 while (mpz_sgn (d) == 0); 196 dn = SIZ (d); 197 dup = PTR (d); 198 MPN_COPY (dnp, dup, dn); 199 dnp[dn - 1] |= GMP_NUMB_HIGHBIT; 200 201 if (test % 2 == 0) 202 { 203 RANDFUNC (n, rands, nbits); 204 nn = SIZ (n); 205 ASSERT_ALWAYS (nn >= dn); 206 } 207 else 208 { 209 do 210 { 211 RANDFUNC (q, rands, urandom () % (nbits - dbits + 1)); 212 RANDFUNC (r, rands, urandom () % mpz_sizeinbase (d, 2)); 213 mpz_mul (n, q, d); 214 mpz_add (n, n, r); 215 nn = SIZ (n); 216 } 217 while (nn > maxnn || nn < dn); 218 } 219 220 ASSERT_ALWAYS (nn <= maxnn); 221 ASSERT_ALWAYS (dn <= maxdn); 222 223 mpz_urandomb (junk, rands, nbits); 224 junkp = PTR (junk); 225 226 np = PTR (n); 227 228 mpz_urandomb (tz, rands, 32); 229 t = mpz_get_ui (tz); 230 231 if (t % 17 == 0) 232 { 233 dnp[dn - 1] = GMP_NUMB_MAX; 234 dup[dn - 1] = GMP_NUMB_MAX; 235 } 236 237 switch ((int) t % 16) 238 { 239 case 0: 240 clearn = urandom () % nn; 241 for (i = clearn; i < nn; i++) 242 np[i] = 0; 243 break; 244 case 1: 245 mpn_sub_1 (np + nn - dn, dnp, dn, urandom ()); 246 break; 247 case 2: 248 mpn_add_1 (np + nn - dn, dnp, dn, urandom ()); 249 break; 250 } 251 252 if (dn >= 2) 253 invert_pi1 (dinv, dnp[dn - 1], dnp[dn - 2]); 254 255 rran0 = urandom (); 256 rran1 = urandom (); 257 qran0 = urandom (); 258 qran1 = urandom (); 259 260 qp[-1] = qran0; 261 qp[nn - dn + 1] = qran1; 262 rp[-1] = rran0; 263 264 ran = urandom (); 265 266 if ((double) (nn - dn) * dn < 1e5) 267 { 268 /* Test mpn_sbpi1_div_qr */ 269 if (dn > 2) 270 { 271 MPN_COPY (rp, np, nn); 272 if (nn > dn) 273 MPN_COPY (qp, junkp, nn - dn); 274 qp[nn - dn] = mpn_sbpi1_div_qr (qp, rp, nn, dnp, dn, dinv.inv32); 275 check_one (qp, rp, np, nn, dnp, dn, "mpn_sbpi1_div_qr", 0); 276 } 277 278 /* Test mpn_sbpi1_divappr_q */ 279 if (dn > 2) 280 { 281 MPN_COPY (rp, np, nn); 282 if (nn > dn) 283 MPN_COPY (qp, junkp, nn - dn); 284 qp[nn - dn] = mpn_sbpi1_divappr_q (qp, rp, nn, dnp, dn, dinv.inv32); 285 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_divappr_q", 1); 286 } 287 288 /* Test mpn_sbpi1_div_q */ 289 if (dn > 2) 290 { 291 MPN_COPY (rp, np, nn); 292 if (nn > dn) 293 MPN_COPY (qp, junkp, nn - dn); 294 qp[nn - dn] = mpn_sbpi1_div_q (qp, rp, nn, dnp, dn, dinv.inv32); 295 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_div_q", 0); 296 } 297 298 /* Test mpn_sec_div_qr */ 299 itch = mpn_sec_div_qr_itch (nn, dn); 300 if (itch + 1 > alloc) 301 { 302 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 303 alloc = itch + 1; 304 } 305 scratch[itch] = ran; 306 MPN_COPY (rp, np, nn); 307 if (nn >= dn) 308 MPN_COPY (qp, junkp, nn - dn + 1); 309 qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dup, dn, scratch); 310 ASSERT_ALWAYS (ran == scratch[itch]); 311 check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_qr (unnorm)", 0); 312 313 /* Test mpn_sec_div_r */ 314 itch = mpn_sec_div_r_itch (nn, dn); 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 mpn_sec_div_r (rp, nn, dup, dn, scratch); 323 ASSERT_ALWAYS (ran == scratch[itch]); 324 /* Note: Since check_one cannot cope with remainder-only functions, we 325 pass qp[] from the previous function, mpn_sec_div_qr. */ 326 check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_r (unnorm)", 0); 327 328 /* Normalised case, mpn_sec_div_qr */ 329 itch = mpn_sec_div_qr_itch (nn, dn); 330 scratch[itch] = ran; 331 332 MPN_COPY (rp, np, nn); 333 if (nn >= dn) 334 MPN_COPY (qp, junkp, nn - dn + 1); 335 qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dnp, dn, scratch); 336 ASSERT_ALWAYS (ran == scratch[itch]); 337 check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_qr (norm)", 0); 338 339 /* Normalised case, mpn_sec_div_r */ 340 itch = mpn_sec_div_r_itch (nn, dn); 341 scratch[itch] = ran; 342 MPN_COPY (rp, np, nn); 343 mpn_sec_div_r (rp, nn, dnp, dn, scratch); 344 ASSERT_ALWAYS (ran == scratch[itch]); 345 /* Note: Since check_one cannot cope with remainder-only functions, we 346 pass qp[] from the previous function, mpn_sec_div_qr. */ 347 check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_r (norm)", 0); 348 } 349 350 /* Test mpn_dcpi1_div_qr */ 351 if (dn >= 6 && nn - dn >= 3) 352 { 353 MPN_COPY (rp, np, nn); 354 if (nn > dn) 355 MPN_COPY (qp, junkp, nn - dn); 356 qp[nn - dn] = mpn_dcpi1_div_qr (qp, rp, nn, dnp, dn, &dinv); 357 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 358 ASSERT_ALWAYS (rp[-1] == rran0); 359 check_one (qp, rp, np, nn, dnp, dn, "mpn_dcpi1_div_qr", 0); 360 } 361 362 /* Test mpn_dcpi1_divappr_q */ 363 if (dn >= 6 && nn - dn >= 3) 364 { 365 MPN_COPY (rp, np, nn); 366 if (nn > dn) 367 MPN_COPY (qp, junkp, nn - dn); 368 qp[nn - dn] = mpn_dcpi1_divappr_q (qp, rp, nn, dnp, dn, &dinv); 369 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 370 ASSERT_ALWAYS (rp[-1] == rran0); 371 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_divappr_q", 1); 372 } 373 374 /* Test mpn_dcpi1_div_q */ 375 if (dn >= 6 && nn - dn >= 3) 376 { 377 MPN_COPY (rp, np, nn); 378 if (nn > dn) 379 MPN_COPY (qp, junkp, nn - dn); 380 qp[nn - dn] = mpn_dcpi1_div_q (qp, rp, nn, dnp, dn, &dinv); 381 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 382 ASSERT_ALWAYS (rp[-1] == rran0); 383 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_div_q", 0); 384 } 385 386 /* Test mpn_mu_div_qr */ 387 if (nn - dn > 2 && dn >= 2) 388 { 389 itch = mpn_mu_div_qr_itch (nn, dn, 0); 390 if (itch + 1 > alloc) 391 { 392 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 393 alloc = itch + 1; 394 } 395 scratch[itch] = ran; 396 MPN_COPY (qp, junkp, nn - dn); 397 MPN_ZERO (rp, dn); 398 rp[dn] = rran1; 399 qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dnp, dn, scratch); 400 ASSERT_ALWAYS (ran == scratch[itch]); 401 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 402 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1); 403 check_one (qp, rp, np, nn, dnp, dn, "mpn_mu_div_qr", 0); 404 } 405 406 /* Test mpn_mu_divappr_q */ 407 if (nn - dn > 2 && dn >= 2) 408 { 409 itch = mpn_mu_divappr_q_itch (nn, dn, 0); 410 if (itch + 1 > alloc) 411 { 412 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 413 alloc = itch + 1; 414 } 415 scratch[itch] = ran; 416 MPN_COPY (qp, junkp, nn - dn); 417 qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dnp, dn, scratch); 418 ASSERT_ALWAYS (ran == scratch[itch]); 419 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 420 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_divappr_q", 4); 421 } 422 423 /* Test mpn_mu_div_q */ 424 if (nn - dn > 2 && dn >= 2) 425 { 426 itch = mpn_mu_div_q_itch (nn, dn, 0); 427 if (itch + 1> alloc) 428 { 429 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 430 alloc = itch + 1; 431 } 432 scratch[itch] = ran; 433 MPN_COPY (qp, junkp, nn - dn); 434 qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dnp, dn, scratch); 435 ASSERT_ALWAYS (ran == scratch[itch]); 436 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 437 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_div_q", 0); 438 } 439 440 if (1) 441 { 442 itch = nn + 1; 443 if (itch + 1> alloc) 444 { 445 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 446 alloc = itch + 1; 447 } 448 scratch[itch] = ran; 449 mpn_div_q (qp, np, nn, dup, dn, scratch); 450 ASSERT_ALWAYS (ran == scratch[itch]); 451 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 452 check_one (qp, NULL, np, nn, dup, dn, "mpn_div_q", 0); 453 } 454 455 if (dn >= 2 && nn >= 2) 456 { 457 mp_limb_t qh; 458 459 /* mpn_divrem_2 */ 460 MPN_COPY (rp, np, nn); 461 qp[nn - 2] = qp[nn-1] = qran1; 462 463 qh = mpn_divrem_2 (qp, 0, rp, nn, dnp + dn - 2); 464 ASSERT_ALWAYS (qp[nn - 2] == qran1); 465 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1); 466 qp[nn - 2] = qh; 467 check_one (qp, rp, np, nn, dnp + dn - 2, 2, "mpn_divrem_2", 0); 468 469 /* Missing: divrem_2 with fraction limbs. */ 470 471 /* mpn_div_qr_2 */ 472 qp[nn - 2] = qran1; 473 474 qh = mpn_div_qr_2 (qp, rp, np, nn, dup + dn - 2); 475 ASSERT_ALWAYS (qp[nn - 2] == qran1); 476 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1); 477 qp[nn - 2] = qh; 478 check_one (qp, rp, np, nn, dup + dn - 2, 2, "mpn_div_qr_2", 0); 479 } 480 if (dn >= 1 && nn >= 1) 481 { 482 /* mpn_div_qr_1 */ 483 mp_limb_t qh; 484 qp[nn-1] = qran1; 485 rp[0] = mpn_div_qr_1 (qp, &qh, np, nn, dnp[dn - 1]); 486 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1); 487 qp[nn - 1] = qh; 488 check_one (qp, rp, np, nn, dnp + dn - 1, 1, "mpn_div_qr_1", 0); 489 } 490 } 491 492 __GMP_FREE_FUNC_LIMBS (scratch, alloc); 493 494 TMP_FREE; 495 496 mpz_clear (n); 497 mpz_clear (d); 498 mpz_clear (q); 499 mpz_clear (r); 500 mpz_clear (tz); 501 mpz_clear (junk); 502 503 tests_end (); 504 return 0; 505 }