tests/mpn/t-bdiv.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 void
  57 check_one (mp_ptr qp, mp_srcptr rp, mp_limb_t rh,
  58            mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, char *fname)
  59 {
  60   mp_size_t qn;
  61   int cmp;
  62   mp_ptr tp;
  63   mp_limb_t cy = 4711;          /* silence warnings */
  64   TMP_DECL;
  65
  66   qn = nn - dn;
  67
  68   if (qn == 0)
  69     return;
  70
  71   TMP_MARK;
  72
  73   tp = TMP_ALLOC_LIMBS (nn + 1);
  74
  75   if (dn >= qn)
  76     mpn_mul (tp, dp, dn, qp, qn);
  77   else
  78     mpn_mul (tp, qp, qn, dp, dn);
  79
  80   if (rp != NULL)
  81     {
  82       cy = mpn_add_n (tp + qn, tp + qn, rp, dn);
  83       cmp = cy != rh || mpn_cmp (tp, np, nn) != 0;
  84     }
  85   else
  86     cmp = mpn_cmp (tp, np, nn - dn) != 0;
  87
  88   if (cmp != 0)
  89     {
  90       printf ("\r*******************************************************************************\n");
  91       printf ("%s inconsistent in test %lu\n", fname, test);
  92       printf ("N=   "); dumpy (np, nn);
  93       printf ("D=   "); dumpy (dp, dn);
  94       printf ("Q=   "); dumpy (qp, qn);
  95       if (rp != NULL)
  96         {
  97           printf ("R=   "); dumpy (rp, dn);
  98           printf ("Rb=  %d, Cy=%d\n", (int) cy, (int) rh);
  99         }
 100       printf ("T=   "); dumpy (tp, nn);
 101       printf ("nn = %ld, dn = %ld, qn = %ld", nn, dn, qn);
 102       printf ("\n*******************************************************************************\n");
 103       abort ();
 104     }
 105
 106   TMP_FREE;
 107 }
 108
 109
 110 /* These are *bit* sizes. */
 111 #define SIZE_LOG 16
 112 #define MAX_DN (1L << SIZE_LOG)
 113 #define MAX_NN (1L << (SIZE_LOG + 1))
 114
 115 #define COUNT 500
 116
 117 mp_limb_t
 118 random_word (gmp_randstate_ptr rs)
 119 {
 120   mpz_t x;
 121   mp_limb_t r;
 122   TMP_DECL;
 123   TMP_MARK;
 124
 125   MPZ_TMP_INIT (x, 2);
 126   mpz_urandomb (x, rs, 32);
 127   r = mpz_get_ui (x);
 128   TMP_FREE;
 129   return r;
 130 }
 131
 132 int
 133 main (int argc, char **argv)
 134 {
 135   gmp_randstate_ptr rands;
 136   unsigned long maxnbits, maxdbits, nbits, dbits;
 137   mpz_t n, d, tz;
 138   mp_size_t maxnn, maxdn, nn, dn, clearn, i;
 139   mp_ptr np, dp, qp, rp;
 140   mp_limb_t rh;
 141   mp_limb_t t;
 142   mp_limb_t dinv;
 143   int count = COUNT;
 144   mp_ptr scratch;
 145   mp_limb_t ran;
 146   mp_size_t alloc, itch;
 147   mp_limb_t rran0, rran1, qran0, qran1;
 148   TMP_DECL;
 149
 150   if (argc > 1)
 151     {
 152       char *end;
 153       count = strtol (argv[1], &end, 0);
 154       if (*end || count <= 0)
 155         {
 156           fprintf (stderr, "Invalid test count: %s.\n", argv[1]);
 157           return 1;
 158         }
 159     }
 160
 161
 162   maxdbits = MAX_DN;
 163   maxnbits = MAX_NN;
 164
 165   tests_start ();
 166   rands = RANDS;
 167
 168   mpz_init (n);
 169   mpz_init (d);
 170   mpz_init (tz);
 171
 172   maxnn = maxnbits / GMP_NUMB_BITS + 1;
 173   maxdn = maxdbits / GMP_NUMB_BITS + 1;
 174
 175   TMP_MARK;
 176
 177   qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
 178   rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
 179
 180   alloc = 1;
 181   scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);
 182
 183   for (test = 0; test < count;)
 184     {
 185       nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;
 186       if (maxdbits > nbits)
 187         dbits = random_word (rands) % nbits + 1;
 188       else
 189         dbits = random_word (rands) % maxdbits + 1;
 190
 191 #if RAND_UNIFORM
 192 #define RANDFUNC mpz_urandomb
 193 #else
 194 #define RANDFUNC mpz_rrandomb
 195 #endif
 196
 197       do
 198         {
 199           RANDFUNC (n, rands, nbits);
 200           do
 201             {
 202               RANDFUNC (d, rands, dbits);
 203             }
 204           while (mpz_sgn (d) == 0);
 205
 206           np = PTR (n);
 207           dp = PTR (d);
 208           nn = SIZ (n);
 209           dn = SIZ (d);
 210         }
 211       while (nn < dn);
 212
 213       dp[0] |= 1;
 214
 215       mpz_urandomb (tz, rands, 32);
 216       t = mpz_get_ui (tz);
 217
 218       if (t % 17 == 0)
 219         dp[0] = GMP_NUMB_MAX;
 220
 221       switch (t % 16)
 222         {
 223         case 0:
 224           clearn = random_word (rands) % nn;
 225           for (i = 0; i <= clearn; i++)
 226             np[i] = 0;
 227           break;
 228         case 1:
 229           mpn_sub_1 (np + nn - dn, dp, dn, random_word (rands));
 230           break;
 231         case 2:
 232           mpn_add_1 (np + nn - dn, dp, dn, random_word (rands));
 233           break;
 234         }
 235
 236       test++;
 237
 238       binvert_limb (dinv, dp[0]);
 239
 240       rran0 = random_word (rands);
 241       rran1 = random_word (rands);
 242       qran0 = random_word (rands);
 243       qran1 = random_word (rands);
 244
 245       qp[-1] = qran0;
 246       qp[nn - dn + 1] = qran1;
 247       rp[-1] = rran0;
 248
 249       ran = random_word (rands);
 250
 251       if ((double) (nn - dn) * dn < 1e5)
 252         {
 253           if (nn > dn)
 254             {
 255               /* Test mpn_sbpi1_bdiv_qr */
 256               MPN_ZERO (qp, nn - dn);
 257               MPN_ZERO (rp, dn);
 258               MPN_COPY (rp, np, nn);
 259               rh = mpn_sbpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
 260               ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 261               ASSERT_ALWAYS (rp[-1] == rran0);
 262               check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_sbpi1_bdiv_qr");
 263             }
 264
 265           if (nn > dn)
 266             {
 267               /* Test mpn_sbpi1_bdiv_q */
 268               MPN_COPY (rp, np, nn);
 269               MPN_ZERO (qp, nn - dn);
 270               mpn_sbpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
 271               ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 272               ASSERT_ALWAYS (rp[-1] == rran0);
 273               check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_sbpi1_bdiv_q");
 274             }
 275         }
 276
 277       if (dn >= 4 && nn - dn >= 2)
 278         {
 279           /* Test mpn_dcpi1_bdiv_qr */
 280           MPN_COPY (rp, np, nn);
 281           MPN_ZERO (qp, nn - dn);
 282           rh = mpn_dcpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
 283           ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 284           ASSERT_ALWAYS (rp[-1] == rran0);
 285           check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_dcpi1_bdiv_qr");
 286         }
 287
 288       if (dn >= 4 && nn - dn >= 2)
 289         {
 290           /* Test mpn_dcpi1_bdiv_q */
 291           MPN_COPY (rp, np, nn);
 292           MPN_ZERO (qp, nn - dn);
 293           mpn_dcpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
 294           ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 295           ASSERT_ALWAYS (rp[-1] == rran0);
 296           check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_dcpi1_bdiv_q");
 297         }
 298
 299       if (nn - dn < 2 || dn < 2)
 300         continue;
 301
 302       /* Test mpn_mu_bdiv_qr */
 303       itch = mpn_mu_bdiv_qr_itch (nn, dn);
 304       if (itch + 1 > alloc)
 305         {
 306           scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
 307           alloc = itch + 1;
 308         }
 309       scratch[itch] = ran;
 310       MPN_ZERO (qp, nn - dn);
 311       MPN_ZERO (rp, dn);
 312       rp[dn] = rran1;
 313       rh = mpn_mu_bdiv_qr (qp, rp, np, nn, dp, dn, scratch);
 314       ASSERT_ALWAYS (ran == scratch[itch]);
 315       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 316       ASSERT_ALWAYS (rp[-1] == rran0);  ASSERT_ALWAYS (rp[dn] == rran1);
 317       check_one (qp, rp, rh, np, nn, dp, dn, "mpn_mu_bdiv_qr");
 318
 319       /* Test mpn_mu_bdiv_q */
 320       itch = mpn_mu_bdiv_q_itch (nn, dn);
 321       if (itch + 1 > alloc)
 322         {
 323           scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
 324           alloc = itch + 1;
 325         }
 326       scratch[itch] = ran;
 327       MPN_ZERO (qp, nn - dn + 1);
 328       mpn_mu_bdiv_q (qp, np, nn - dn, dp, dn, scratch);
 329       ASSERT_ALWAYS (ran == scratch[itch]);
 330       ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
 331       check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_mu_bdiv_q");
 332     }
 333
 334   __GMP_FREE_FUNC_LIMBS (scratch, alloc);
 335
 336   TMP_FREE;
 337
 338   mpz_clear (n);
 339   mpz_clear (d);
 340   mpz_clear (tz);
 341
 342   tests_end ();
 343   return 0;
 344 }