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