tests/mpn/t-div.c

   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 (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 }