11844Swollman/* Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc. 250476Speter 31638SrgrimesThis program is free software; you can redistribute it and/or modify it under 494940Sruthe terms of the GNU General Public License as published by the Free Software 51638SrgrimesFoundation; either version 3 of the License, or (at your option) any later 6103713Smarkmversion. 71638Srgrimes 8133653SruThis program is distributed in the hope that it will be useful, but WITHOUT ANY 9119607SruWARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 10119607SruPARTICULAR PURPOSE. See the GNU General Public License for more details. 11119607Sru 12119607SruYou should have received a copy of the GNU General Public License along with 13162210Simpthis program. If not, see http://www.gnu.org/licenses/. */ 14162210Simp 15162293Sobrien 16162210Simp#include <stdlib.h> /* for strtol */ 17162210Simp#include <stdio.h> /* for printf */ 18206082Snetchild 19206082Snetchild#include "gmp.h" 20206082Snetchild#include "gmp-impl.h" 21206082Snetchild#include "longlong.h" 22206082Snetchild#include "tests/tests.h" 23119607Sru 24119607Sru 25204027Smarcelstatic void 26179184Sjbdumpy (mp_srcptr p, mp_size_t n) 27179184Sjb{ 28179184Sjb mp_size_t i; 29119607Sru if (n > 20) 30179184Sjb { 31119607Sru for (i = n - 1; i >= n - 4; i--) 32119607Sru { 33117034Sgordon printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 34119607Sru printf (" "); 351638Srgrimes } 362827Sjkh printf ("... "); 371638Srgrimes for (i = 3; i >= 0; i--) 382827Sjkh { 391638Srgrimes printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 40139112Sru printf (" " + (i == 0)); 411844Swollman } 421844Swollman } 431638Srgrimes else 4494424Sru { 4594424Sru for (i = n - 1; i >= 0; i--) 4694424Sru { 4794424Sru printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 481638Srgrimes printf (" " + (i == 0)); 491638Srgrimes } 501638Srgrimes } 5136054Sbde puts (""); 52125119Sru} 53172401Sru 5436054Sbdestatic unsigned long test; 55172401Sru 56172401Srustatic void 57172401Srucheck_one (mp_ptr qp, mp_srcptr rp, 5836054Sbde mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, 59172401Sru char *fname, mp_limb_t q_allowed_err) 6036054Sbde{ 611844Swollman mp_size_t qn = nn - dn + 1; 621638Srgrimes mp_ptr tp; 63212423Srpaulo const char *msg; 64212423Srpaulo const char *tvalue; 65212423Srpaulo mp_limb_t i; 66212423Srpaulo TMP_DECL; 67212423Srpaulo TMP_MARK; 6894518Sru 6994518Sru tp = TMP_ALLOC_LIMBS (nn + 1); 7094518Sru if (dn >= qn) 7194518Sru refmpn_mul (tp, dp, dn, qp, qn); 7294518Sru else 73210612Srpaulo refmpn_mul (tp, qp, qn, dp, dn); 74210636Srpaulo 75210636Srpaulo for (i = 0; i < q_allowed_err && (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0); i++) 7694518Sru ASSERT_NOCARRY (refmpn_sub (tp, tp, nn+1, dp, dn)); 77144893Sharti 781844Swollman if (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0) 7994518Sru { 8094424Sru msg = "q too large"; 8194424Sru tvalue = "Q*D"; 8294424Sru error: 832351Sbde printf ("\r*******************************************************************************\n"); 8494424Sru printf ("%s failed test %lu: %s\n", fname, test, msg); 851638Srgrimes printf ("N= "); dumpy (np, nn); 862351Sbde printf ("D= "); dumpy (dp, dn); 872351Sbde printf ("Q= "); dumpy (qp, qn); 882351Sbde if (rp) 892351Sbde { printf ("R= "); dumpy (rp, dn); } 902351Sbde printf ("%5s=", tvalue); dumpy (tp, nn+1); 912351Sbde printf ("nn = %ld, dn = %ld, qn = %ld\n", nn, dn, qn); 9233624Seivind abort (); 93212423Srpaulo } 94212423Srpaulo 95212423Srpaulo ASSERT_NOCARRY (refmpn_sub_n (tp, np, tp, nn)); 9634081Sbde tvalue = "N-Q*D"; 97212423Srpaulo if (!mpn_zero_p (tp + dn, nn - dn) || mpn_cmp (tp, dp, dn) >= 0) 9894497Sru { 9994497Sru msg = "q too small"; 10094497Sru goto error; 10194410Sru } 10233624Seivind 103210612Srpaulo if (rp && mpn_cmp (rp, tp, dn) != 0) 104210656Srpaulo { 105210636Srpaulo msg = "r incorrect"; 10694518Sru goto error; 1071638Srgrimes } 1081638Srgrimes 1091638Srgrimes TMP_FREE; 110156813Sru} 11174806Sru 1121638Srgrimes 11358493Sru/* These are *bit* sizes. */ 11474806Sru#ifndef SIZE_LOG 11574941Sru#define SIZE_LOG 17 11674941Sru#endif 1171638Srgrimes#define MAX_DN (1L << SIZE_LOG) 1181638Srgrimes#define MAX_NN (1L << (SIZE_LOG + 1)) 1191638Srgrimes 12097769Sru#define COUNT 200 121156813Sru 12296164Srump_limb_t 12374806Srurandom_word (gmp_randstate_ptr rs) 1241638Srgrimes{ 125119607Sru mpz_t x; 126119607Sru mp_limb_t r; 127119607Sru TMP_DECL; 1281638Srgrimes TMP_MARK; 129119607Sru 130119607Sru MPZ_TMP_INIT (x, 2); 131119607Sru mpz_urandomb (x, rs, 32); 132119607Sru r = mpz_get_ui (x); 133125620Sru TMP_FREE; 134125620Sru return r; 13555670Sbde} 13624750Sbde 137125620Sruint 138125620Srumain (int argc, char **argv) 139125620Sru{ 140125620Sru gmp_randstate_ptr rands; 141125620Sru unsigned long maxnbits, maxdbits, nbits, dbits; 14246541Sbde mpz_t n, d, q, r, tz; 14394497Sru mp_size_t maxnn, maxdn, nn, dn, clearn, i; 14494497Sru mp_ptr np, dp, qp, rp; 14524750Sbde mp_limb_t t; 14628945Speter gmp_pi1_t dinv; 147125620Sru int count = COUNT; 14824750Sbde mp_ptr scratch; 1491638Srgrimes mp_limb_t ran; 1501638Srgrimes mp_size_t alloc, itch; 151137614Sru mp_limb_t rran0, rran1, qran0, qran1; 152139111Sru TMP_DECL; 153137164Sru 154137164Sru if (argc > 1) 155137614Sru { 156137614Sru char *end; 157137164Sru count = strtol (argv[1], &end, 0); 15849328Shoek if (*end || count <= 0) 15949328Shoek { 16049328Shoek fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 16149328Shoek return 1; 16249328Shoek } 163125620Sru } 16496163Sru 16599343Sru 16696163Sru maxdbits = MAX_DN; 1671638Srgrimes maxnbits = MAX_NN; 16875083Sru 169100872Sru tests_start (); 17075083Sru rands = RANDS; 17175083Sru 172100872Sru mpz_init (n); 17349328Shoek mpz_init (d); 1741638Srgrimes mpz_init (q); 17575083Sru mpz_init (r); 176144893Sharti mpz_init (tz); 1771638Srgrimes 17875284Sru maxnn = maxnbits / GMP_NUMB_BITS + 1; 17975284Sru maxdn = maxdbits / GMP_NUMB_BITS + 1; 18099343Sru 18175284Sru TMP_MARK; 18275284Sru 18375284Sru qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 18475284Sru rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 18575284Sru 18675284Sru alloc = 1; 18775284Sru scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc); 18875284Sru 18975284Sru for (test = 0; test < count;) 19075284Sru { 19175284Sru do 19275284Sru { 19375284Sru nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS; 19475284Sru if (maxdbits > nbits) 19575284Sru dbits = random_word (rands) % nbits + 1; 19675284Sru else 19788055Sru dbits = random_word (rands) % maxdbits + 1; 19888055Sru } 199100872Sru while (nbits < dbits); 20075284Sru 20194954Sru#if RAND_UNIFORM 20275284Sru#define RANDFUNC mpz_urandomb 20375284Sru#else 20475284Sru#define RANDFUNC mpz_rrandomb 20575284Sru#endif 20699257Sru 20799257Sru do 20899257Sru RANDFUNC (d, rands, dbits); 20997769Sru while (mpz_sgn (d) == 0); 21096668Sru dn = SIZ (d); 21199256Sru dp = PTR (d); 21296462Sru dp[dn - 1] |= GMP_NUMB_HIGHBIT; 213156813Sru 21496164Sru if (test % 2 == 0) 21599343Sru { 21696163Sru RANDFUNC (n, rands, nbits); 21796163Sru nn = SIZ (n); 2181844Swollman ASSERT_ALWAYS (nn >= dn); 2191638Srgrimes } 2201638Srgrimes else 221103713Smarkm { 2221638Srgrimes do 223103713Smarkm { 2241638Srgrimes RANDFUNC (q, rands, random_word (rands) % (nbits - dbits + 1)); 2251638Srgrimes RANDFUNC (r, rands, random_word (rands) % mpz_sizeinbase (d, 2)); 2261638Srgrimes mpz_mul (n, q, d); 227156813Sru mpz_add (n, n, r); 2281638Srgrimes nn = SIZ (n); 22974842Sru } 2301844Swollman while (nn > maxnn || nn < dn); 2311844Swollman } 23234081Sbde 23394113Sru ASSERT_ALWAYS (nn <= maxnn); 23434087Sbde ASSERT_ALWAYS (dn <= maxdn); 23534081Sbde 23634081Sbde np = PTR (n); 23716663Sjkh 23876861Skris mpz_urandomb (tz, rands, 32); 23976861Skris t = mpz_get_ui (tz); 240139068Spaul 241139068Spaul if (t % 17 == 0) 242139068Spaul dp[dn - 1] = GMP_NUMB_MAX; 243139068Spaul 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