contrib/libgmp/mpz/tests/t-mul.c

   1 /* Test mpz_add, mpz_cmp, mpz_cmp_ui, mpz_divmod, mpz_mul.
   2
   3 Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Library General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Library General Public License
  18 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include <stdio.h>
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25 #include "longlong.h"
  26 #include "urandom.h"
  27
  28 void debug_mp ();
  29 mp_size_t _mpn_mul_classic ();
  30 void mpz_refmul ();
  31
  32 #ifndef SIZE
  33 #define SIZE 128
  34 #endif
  35
  36 main (argc, argv)
  37      int argc;
  38      char **argv;
  39 {
  40   mpz_t multiplier, multiplicand;
  41   mpz_t product, ref_product;
  42   mpz_t quotient, remainder;
  43   mp_size_t multiplier_size, multiplicand_size;
  44   int i;
  45   int reps = 10000;
  46
  47   if (argc == 2)
  48      reps = atoi (argv[1]);
  49
  50   mpz_init (multiplier);
  51   mpz_init (multiplicand);
  52   mpz_init (product);
  53   mpz_init (ref_product);
  54   mpz_init (quotient);
  55   mpz_init (remainder);
  56
  57   for (i = 0; i < reps; i++)
  58     {
  59       multiplier_size = urandom () % SIZE - SIZE/2;
  60       multiplicand_size = urandom () % SIZE - SIZE/2;
  61
  62       mpz_random2 (multiplier, multiplier_size);
  63       mpz_random2 (multiplicand, multiplicand_size);
  64
  65       mpz_mul (product, multiplier, multiplicand);
  66       mpz_refmul (ref_product, multiplier, multiplicand);
  67       if (mpz_cmp_ui (multiplicand, 0) != 0)
  68         mpz_divmod (quotient, remainder, product, multiplicand);
  69
  70       if (mpz_cmp (product, ref_product))
  71         dump_abort (multiplier, multiplicand);
  72
  73       if (mpz_cmp_ui (multiplicand, 0) != 0)
  74       if (mpz_cmp_ui (remainder, 0) || mpz_cmp (quotient, multiplier))
  75         dump_abort (multiplier, multiplicand);
  76     }
  77
  78   exit (0);
  79 }
  80
  81 void
  82 mpz_refmul (w, u, v)
  83      mpz_t w;
  84      const mpz_t u;
  85      const mpz_t v;
  86 {
  87   mp_size_t usize = u->_mp_size;
  88   mp_size_t vsize = v->_mp_size;
  89   mp_size_t wsize;
  90   mp_size_t sign_product;
  91   mp_ptr up, vp;
  92   mp_ptr wp;
  93   mp_ptr free_me = NULL;
  94   size_t free_me_size;
  95   TMP_DECL (marker);
  96
  97   TMP_MARK (marker);
  98   sign_product = usize ^ vsize;
  99   usize = ABS (usize);
 100   vsize = ABS (vsize);
 101
 102   if (usize < vsize)
 103     {
 104       /* Swap U and V.  */
 105       {const __mpz_struct *t = u; u = v; v = t;}
 106       {mp_size_t t = usize; usize = vsize; vsize = t;}
 107     }
 108
 109   up = u->_mp_d;
 110   vp = v->_mp_d;
 111   wp = w->_mp_d;
 112
 113   /* Ensure W has space enough to store the result.  */
 114   wsize = usize + vsize;
 115   if (w->_mp_alloc < wsize)
 116     {
 117       if (wp == up || wp == vp)
 118         {
 119           free_me = wp;
 120           free_me_size = w->_mp_alloc;
 121         }
 122       else
 123         (*_mp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB);
 124
 125       w->_mp_alloc = wsize;
 126       wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
 127       w->_mp_d = wp;
 128     }
 129   else
 130     {
 131       /* Make U and V not overlap with W.  */
 132       if (wp == up)
 133         {
 134           /* W and U are identical.  Allocate temporary space for U.  */
 135           up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
 136           /* Is V identical too?  Keep it identical with U.  */
 137           if (wp == vp)
 138             vp = up;
 139           /* Copy to the temporary space.  */
 140           MPN_COPY (up, wp, usize);
 141         }
 142       else if (wp == vp)
 143         {
 144           /* W and V are identical.  Allocate temporary space for V.  */
 145           vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
 146           /* Copy to the temporary space.  */
 147           MPN_COPY (vp, wp, vsize);
 148         }
 149     }
 150
 151   wsize = _mpn_mul_classic (wp, up, usize, vp, vsize);
 152   w->_mp_size = sign_product < 0 ? -wsize : wsize;
 153   if (free_me != NULL)
 154     (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
 155
 156   TMP_FREE (marker);
 157 }
 158
 159 mp_size_t
 160 _mpn_mul_classic (prodp, up, usize, vp, vsize)
 161      mp_ptr prodp;
 162      mp_srcptr up;
 163      mp_size_t usize;
 164      mp_srcptr vp;
 165      mp_size_t vsize;
 166 {
 167   mp_size_t i, j;
 168   mp_limb_t prod_low, prod_high;
 169   mp_limb_t cy_dig;
 170   mp_limb_t v_limb, c;
 171
 172   if (vsize == 0)
 173     return 0;
 174
 175   /* Offset UP and PRODP so that the inner loop can be faster.  */
 176   up += usize;
 177   prodp += usize;
 178
 179   /* Multiply by the first limb in V separately, as the result can
 180      be stored (not added) to PROD.  We also avoid a loop for zeroing.  */
 181   v_limb = vp[0];
 182   cy_dig = 0;
 183   j = -usize;
 184   do
 185     {
 186       umul_ppmm (prod_high, prod_low, up[j], v_limb);
 187       add_ssaaaa (cy_dig, prodp[j], prod_high, prod_low, 0, cy_dig);
 188       j++;
 189     }
 190   while (j < 0);
 191
 192   prodp[j] = cy_dig;
 193   prodp++;
 194
 195   /* For each iteration in the outer loop, multiply one limb from
 196      U with one limb from V, and add it to PROD.  */
 197   for (i = 1; i < vsize; i++)
 198     {
 199       v_limb = vp[i];
 200       cy_dig = 0;
 201       j = -usize;
 202
 203       /* Inner loops.  Simulate the carry flag by jumping between
 204          these loops.  The first is used when there was no carry
 205          in the previois iteration; the second when there was carry.  */
 206
 207       do
 208         {
 209           umul_ppmm (prod_high, prod_low, up[j], v_limb);
 210           add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
 211           c = prodp[j];
 212           prod_low += c;
 213           prodp[j] = prod_low;
 214           if (prod_low < c)
 215             goto cy_loop;
 216         ncy_loop:
 217           j++;
 218         }
 219       while (j < 0);
 220
 221       prodp[j] = cy_dig;
 222       prodp++;
 223       continue;
 224
 225       do
 226         {
 227           umul_ppmm (prod_high, prod_low, up[j], v_limb);
 228           add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
 229           c = prodp[j];
 230           prod_low += c + 1;
 231           prodp[j] = prod_low;
 232           if (prod_low > c)
 233             goto ncy_loop;
 234         cy_loop:
 235           j++;
 236         }
 237       while (j < 0);
 238
 239       cy_dig += 1;
 240       prodp[j] = cy_dig;
 241       prodp++;
 242     }
 243
 244   return usize + vsize - (cy_dig == 0);
 245 }
 246
 247 dump_abort (multiplier, multiplicand)
 248      mpz_t multiplier, multiplicand;
 249 {
 250   fprintf (stderr, "ERROR\n");
 251   fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
 252   fprintf (stderr, "multiplicand  = "); debug_mp (multiplicand, -16);
 253   abort();
 254 }
 255
 256 void
 257 debug_mp (x, base)
 258      mpz_t x;
 259 {
 260   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
 261 }