tests/mpz/t-mul.c

   1 /* Test mpz_cmp, mpz_mul.
   2
   3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004 Free
   4 Software Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27 #include "tests.h"
  28
  29 void debug_mp __GMP_PROTO ((mpz_t));
  30 static void refmpz_mul __GMP_PROTO ((mpz_t, const mpz_t, const mpz_t));
  31 void dump_abort __GMP_PROTO ((int, char *, mpz_t, mpz_t, mpz_t, mpz_t));
  32
  33 #define FFT_MIN_BITSIZE 100000
  34
  35 char *extra_fft;
  36
  37 void
  38 one (int i, mpz_t multiplicand, mpz_t multiplier)
  39 {
  40   mpz_t product, ref_product;
  41
  42   mpz_init (product);
  43   mpz_init (ref_product);
  44
  45   /* Test plain multiplication comparing results against reference code.  */
  46   mpz_mul (product, multiplier, multiplicand);
  47   refmpz_mul (ref_product, multiplier, multiplicand);
  48   if (mpz_cmp (product, ref_product))
  49     dump_abort (i, "incorrect plain product",
  50                 multiplier, multiplicand, product, ref_product);
  51
  52   /* Test squaring, comparing results against plain multiplication  */
  53   mpz_mul (product, multiplier, multiplier);
  54   mpz_set (multiplicand, multiplier);
  55   mpz_mul (ref_product, multiplier, multiplicand);
  56   if (mpz_cmp (product, ref_product))
  57     dump_abort (i, "incorrect square product",
  58                 multiplier, multiplier, product, ref_product);
  59
  60   mpz_clear (product);
  61   mpz_clear (ref_product);
  62 }
  63
  64 int
  65 main (int argc, char **argv)
  66 {
  67   mpz_t op1, op2;
  68   int i;
  69   int fft_max_2exp;
  70
  71   gmp_randstate_ptr rands;
  72   mpz_t bs;
  73   unsigned long bsi, size_range, fsize_range;
  74
  75   tests_start ();
  76   rands = RANDS;
  77
  78   extra_fft = getenv ("GMP_CHECK_FFT");
  79   fft_max_2exp = 0;
  80   if (extra_fft != NULL)
  81     fft_max_2exp = atoi (extra_fft);
  82
  83   if (fft_max_2exp <= 1)        /* compat with old use of GMP_CHECK_FFT */
  84     fft_max_2exp = 22;          /* default limit, good for any machine */
  85
  86   mpz_init (bs);
  87   mpz_init (op1);
  88   mpz_init (op2);
  89
  90   fsize_range = 4 << 8;         /* a fraction 1/256 of size_range */
  91   for (i = 0;; i++)
  92     {
  93       size_range = fsize_range >> 8;
  94       fsize_range = fsize_range * 33 / 32;
  95
  96       if (size_range > fft_max_2exp)
  97         break;
  98
  99       mpz_urandomb (bs, rands, size_range);
 100       mpz_rrandomb (op1, rands, mpz_get_ui (bs));
 101       if (i & 1)
 102         mpz_urandomb (bs, rands, size_range);
 103       mpz_rrandomb (op2, rands, mpz_get_ui (bs));
 104
 105       mpz_urandomb (bs, rands, 4);
 106       bsi = mpz_get_ui (bs);
 107       if ((bsi & 0x3) == 0)
 108         mpz_neg (op1, op1);
 109       if ((bsi & 0xC) == 0)
 110         mpz_neg (op2, op2);
 111
 112       /* printf ("%d %d\n", SIZ (op1), SIZ (op2)); */
 113       one (i, op2, op1);
 114     }
 115
 116   for (i = -50; i < 0; i++)
 117     {
 118       mpz_urandomb (bs, rands, 32);
 119       size_range = mpz_get_ui (bs) % fft_max_2exp;
 120
 121       mpz_urandomb (bs, rands, size_range);
 122       mpz_rrandomb (op1, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
 123       mpz_urandomb (bs, rands, size_range);
 124       mpz_rrandomb (op2, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
 125
 126       /* printf ("%d: %d %d\n", i, SIZ (op1), SIZ (op2)); */
 127       fflush (stdout);
 128       one (-1, op2, op1);
 129     }
 130
 131   mpz_clear (bs);
 132   mpz_clear (op1);
 133   mpz_clear (op2);
 134
 135   tests_end ();
 136   exit (0);
 137 }
 138
 139 static void
 140 refmpz_mul (mpz_t w, const mpz_t u, const mpz_t v)
 141 {
 142   mp_size_t usize = u->_mp_size;
 143   mp_size_t vsize = v->_mp_size;
 144   mp_size_t wsize;
 145   mp_size_t sign_product;
 146   mp_ptr up, vp;
 147   mp_ptr wp;
 148   mp_size_t talloc;
 149
 150   sign_product = usize ^ vsize;
 151   usize = ABS (usize);
 152   vsize = ABS (vsize);
 153
 154   if (usize == 0 || vsize == 0)
 155     {
 156       SIZ (w) = 0;
 157       return;
 158     }
 159
 160   talloc = usize + vsize;
 161
 162   up = u->_mp_d;
 163   vp = v->_mp_d;
 164
 165   wp = __GMP_ALLOCATE_FUNC_LIMBS (talloc);
 166
 167   if (usize > vsize)
 168     refmpn_mul (wp, up, usize, vp, vsize);
 169   else
 170     refmpn_mul (wp, vp, vsize, up, usize);
 171   wsize = usize + vsize;
 172   wsize -= wp[wsize - 1] == 0;
 173   MPZ_REALLOC (w, wsize);
 174   MPN_COPY (PTR(w), wp, wsize);
 175
 176   SIZ(w) = sign_product < 0 ? -wsize : wsize;
 177   __GMP_FREE_FUNC_LIMBS (wp, talloc);
 178 }
 179
 180 void
 181 dump_abort (int i, char *s,
 182             mpz_t op1, mpz_t op2, mpz_t product, mpz_t ref_product)
 183 {
 184   mp_size_t b, e;
 185   fprintf (stderr, "ERROR: %s in test %d\n", s, i);
 186   fprintf (stderr, "op1          = "); debug_mp (op1);
 187   fprintf (stderr, "op2          = "); debug_mp (op2);
 188   fprintf (stderr, "    product  = "); debug_mp (product);
 189   fprintf (stderr, "ref_product  = "); debug_mp (ref_product);
 190   for (b = 0; b < ABSIZ(ref_product); b++)
 191     if (PTR(ref_product)[b] != PTR(product)[b])
 192       break;
 193   for (e = ABSIZ(ref_product) - 1; e >= 0; e--)
 194     if (PTR(ref_product)[e] != PTR(product)[e])
 195       break;
 196   printf ("ERRORS in %ld--%ld\n", b, e);
 197   abort();
 198 }
 199
 200 void
 201 debug_mp (mpz_t x)
 202 {
 203   size_t siz = mpz_sizeinbase (x, 16);
 204
 205   if (siz > 65)
 206     {
 207       mpz_t q;
 208       mpz_init (q);
 209       mpz_tdiv_q_2exp (q, x, 4 * (mpz_sizeinbase (x, 16) - 25));
 210       gmp_fprintf (stderr, "%ZX...", q);
 211       mpz_tdiv_r_2exp (q, x, 4 * 25);
 212       gmp_fprintf (stderr, "%025ZX [%d]\n", q, (int) siz);
 213       mpz_clear (q);
 214     }
 215   else
 216     {
 217       gmp_fprintf (stderr, "%ZX\n", x);
 218     }
 219 }