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 2.1 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; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "mpir.h"
27 #include "gmp-impl.h"
28 #include "longlong.h"
29 #include "tests.h"
30
31 void debug_mp _PROTO ((mpz_t));
32 static void ref_mpn_mul _PROTO ((mp_ptr,mp_srcptr,mp_size_t,mp_srcptr,mp_size_t));
33 static void ref_mpz_mul _PROTO ((mpz_t, const mpz_t, const mpz_t));
34 void dump_abort _PROTO ((int, char *, mpz_t, mpz_t, mpz_t, mpz_t));
35
36 #define FFT_MIN_BITSIZE 100000
37
38 char *extra_fft;
39
40 void
one(int i,mpz_t multiplicand,mpz_t multiplier)41 one (int i, mpz_t multiplicand, mpz_t multiplier)
42 {
43 mpz_t product, ref_product;
44 mpz_t quotient;
45
46 mpz_init (product);
47 mpz_init (ref_product);
48 mpz_init (quotient);
49
50 /* Test plain multiplication comparing results against reference code. */
51 mpz_mul (product, multiplier, multiplicand);
52 ref_mpz_mul (ref_product, multiplier, multiplicand);
53 if (mpz_cmp (product, ref_product))
54 dump_abort (i, "incorrect plain product",
55 multiplier, multiplicand, product, ref_product);
56
57 /* Test squaring, comparing results against plain multiplication */
58 mpz_mul (product, multiplier, multiplier);
59 mpz_set (multiplicand, multiplier);
60 mpz_mul (ref_product, multiplier, multiplicand);
61 if (mpz_cmp (product, ref_product))
62 dump_abort (i, "incorrect square product",
63 multiplier, multiplier, product, ref_product);
64
65 mpz_clear (product);
66 mpz_clear (ref_product);
67 mpz_clear (quotient);
68 }
69
70 int
main(int argc,char ** argv)71 main (int argc, char **argv)
72 {
73 mpz_t op1, op2;
74 int i;
75
76 gmp_randstate_t rands;
77 mpz_t bs;
78 unsigned long bsi, size_range, fsize_range;
79
80 tests_start ();
81 gmp_randinit_default(rands);
82
83 extra_fft = getenv ("GMP_CHECK_FFT");
84
85 mpz_init (bs);
86 mpz_init (op1);
87 mpz_init (op2);
88
89 fsize_range = 4 << 8; /* a fraction 1/256 of size_range */
90 for (i = 0; fsize_range >> 8 < (extra_fft ? 27 : 22); i++)
91 {
92 size_range = fsize_range >> 8;
93 fsize_range = fsize_range * 33 / 32;
94
95 mpz_urandomb (bs, rands, size_range);
96 mpz_rrandomb (op1, rands, mpz_get_ui (bs));
97 mpz_urandomb (bs, rands, size_range);
98 mpz_rrandomb (op2, rands, mpz_get_ui (bs));
99
100 mpz_urandomb (bs, rands, 4);
101 bsi = mpz_get_ui (bs);
102 if ((bsi & 0x3) == 0)
103 mpz_neg (op1, op1);
104 if ((bsi & 0xC) == 0)
105 mpz_neg (op2, op2);
106
107 /* printf ("%d %d\n", SIZ (op1), SIZ (op2)); */
108 one (i, op2, op1);
109 }
110
111 if (extra_fft)
112 for (i = -50; i < 0; i++)
113 {
114 mpz_urandomb (bs, rands, 32);
115 size_range = mpz_get_ui (bs) % 27;
116
117 mpz_urandomb (bs, rands, size_range);
118 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
119 mpz_urandomb (bs, rands, size_range);
120 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
121
122 /* printf ("%d: %d %d\n", i, SIZ (op1), SIZ (op2)); */
123 fflush (stdout);
124 one (-1, op2, op1);
125 }
126
127 mpz_clear (bs);
128 mpz_clear (op1);
129 mpz_clear (op2);
130 gmp_randclear(rands);
131 tests_end ();
132 exit (0);
133 }
134
135 static void
ref_mpz_mul(mpz_t w,const mpz_t u,const mpz_t v)136 ref_mpz_mul (mpz_t w, const mpz_t u, const mpz_t v)
137 {
138 mp_size_t usize = u->_mp_size;
139 mp_size_t vsize = v->_mp_size;
140 mp_size_t wsize;
141 mp_size_t sign_product;
142 mp_ptr up, vp;
143 mp_ptr wp;
144 mp_size_t talloc;
145
146 sign_product = usize ^ vsize;
147 usize = ABS (usize);
148 vsize = ABS (vsize);
149
150 if (usize == 0 || vsize == 0)
151 {
152 SIZ (w) = 0;
153 return;
154 }
155
156 talloc = usize + vsize;
157
158 up = u->_mp_d;
159 vp = v->_mp_d;
160
161 wp = __GMP_ALLOCATE_FUNC_LIMBS (talloc);
162
163 if (usize > vsize)
164 ref_mpn_mul (wp, up, usize, vp, vsize);
165 else
166 ref_mpn_mul (wp, vp, vsize, up, usize);
167 wsize = usize + vsize;
168 wsize -= wp[wsize - 1] == 0;
169 MPZ_REALLOC (w, wsize);
170 MPN_COPY (PTR(w), wp, wsize);
171
172 SIZ(w) = sign_product < 0 ? -wsize : wsize;
173 __GMP_FREE_FUNC_LIMBS (wp, talloc);
174 }
175
176 static void mul_basecase __GMP_PROTO ((mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t));
177
178 #define TOOM3_THRESHOLD (MAX (MUL_TOOM3_THRESHOLD, SQR_TOOM3_THRESHOLD))
179 #define FFT_THRESHOLD (MAX (MUL_FFT_FULL_THRESHOLD, SQR_FFT_FULL_THRESHOLD))
180
181 static void
ref_mpn_mul(mp_ptr wp,mp_srcptr up,mp_size_t un,mp_srcptr vp,mp_size_t vn)182 ref_mpn_mul (mp_ptr wp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
183 {
184 mp_ptr tp;
185 mp_size_t tn;
186 mp_limb_t cy;
187
188 if (vn < TOOM3_THRESHOLD)
189 {
190 /* In the mpn_mul_basecase and mpn_kara_mul_n range, use our own
191 mul_basecase. */
192 if (vn != 0)
193 mul_basecase (wp, up, un, vp, vn);
194 else
195 MPN_ZERO (wp, un);
196 return;
197 }
198
199 if (vn < FFT_THRESHOLD)
200 {
201 /* In the mpn_toom3_mul_n and mpn_toom4_mul_n range, use mpn_kara_mul_n. */
202 tn = 2 * vn + MPN_KARA_MUL_N_TSIZE (vn);
203 tp = __GMP_ALLOCATE_FUNC_LIMBS (tn);
204 mpn_kara_mul_n (tp, up, vp, vn, tp + 2 * vn);
205 }
206 else
207 {
208 /* Finally, for the largest operands, use mpn_toom3_mul_n. */
209 /* The "- 63 + 255" tweaks the allocation to allow for huge operands.
210 See the definition of this macro in gmp-impl.h to understand this. */
211 tn = 2 * vn + MPN_TOOM3_MUL_N_TSIZE (vn) - 63 + 255;
212 tp = __GMP_ALLOCATE_FUNC_LIMBS (tn);
213 mpn_toom3_mul_n (tp, up, vp, vn, tp + 2 * vn);
214 }
215
216 if (un != vn)
217 {
218 if (un - vn < vn)
219 ref_mpn_mul (wp + vn, vp, vn, up + vn, un - vn);
220 else
221 ref_mpn_mul (wp + vn, up + vn, un - vn, vp, vn);
222
223 MPN_COPY (wp, tp, vn);
224 cy = mpn_add_n (wp + vn, wp + vn, tp + vn, vn);
225 mpn_incr_u (wp + 2 * vn, cy);
226 }
227 else
228 {
229 MPN_COPY (wp, tp, 2 * vn);
230 }
231
232 __GMP_FREE_FUNC_LIMBS (tp, tn);
233 }
234
235 static void
mul_basecase(mp_ptr wp,mp_srcptr up,mp_size_t un,mp_srcptr vp,mp_size_t vn)236 mul_basecase (mp_ptr wp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
237 {
238 mp_size_t i, j;
239 mp_limb_t prod_low, prod_high;
240 mp_limb_t cy_dig;
241 mp_limb_t v_limb;
242
243 /* Multiply by the first limb in V separately, as the result can
244 be stored (not added) to PROD. We also avoid a loop for zeroing. */
245 v_limb = vp[0];
246 cy_dig = 0;
247 for (j = un; j > 0; j--)
248 {
249 mp_limb_t u_limb, w_limb;
250 u_limb = *up++;
251 umul_ppmm (prod_high, prod_low, u_limb, v_limb << GMP_NAIL_BITS);
252 add_ssaaaa (cy_dig, w_limb, prod_high, prod_low, 0, cy_dig << GMP_NAIL_BITS);
253 *wp++ = w_limb >> GMP_NAIL_BITS;
254 }
255
256 *wp++ = cy_dig;
257 wp -= un;
258 up -= un;
259
260 /* For each iteration in the outer loop, multiply one limb from
261 U with one limb from V, and add it to PROD. */
262 for (i = 1; i < vn; i++)
263 {
264 v_limb = vp[i];
265 cy_dig = 0;
266
267 for (j = un; j > 0; j--)
268 {
269 mp_limb_t u_limb, w_limb;
270 u_limb = *up++;
271 umul_ppmm (prod_high, prod_low, u_limb, v_limb << GMP_NAIL_BITS);
272 w_limb = *wp;
273 add_ssaaaa (prod_high, prod_low, prod_high, prod_low, 0, w_limb << GMP_NAIL_BITS);
274 prod_low >>= GMP_NAIL_BITS;
275 prod_low += cy_dig;
276 #if GMP_NAIL_BITS == 0
277 cy_dig = prod_high + (prod_low < cy_dig);
278 #else
279 cy_dig = prod_high;
280 cy_dig += prod_low >> GMP_NUMB_BITS;
281 #endif
282 *wp++ = prod_low & GMP_NUMB_MASK;
283 }
284
285 *wp++ = cy_dig;
286 wp -= un;
287 up -= un;
288 }
289 }
290
291 void
dump_abort(int i,char * s,mpz_t op1,mpz_t op2,mpz_t product,mpz_t ref_product)292 dump_abort (int i, char *s,
293 mpz_t op1, mpz_t op2, mpz_t product, mpz_t ref_product)
294 {
295 fprintf (stderr, "ERROR: %s in test %d\n", s, i);
296 fprintf (stderr, "op1 = "); debug_mp (op1);
297 fprintf (stderr, "op2 = "); debug_mp (op2);
298 fprintf (stderr, " product = "); debug_mp (product);
299 fprintf (stderr, "ref_product = "); debug_mp (ref_product);
300 abort();
301 }
302
303 void
debug_mp(mpz_t x)304 debug_mp (mpz_t x)
305 {
306 size_t siz = mpz_sizeinbase (x, 16);
307
308 if (siz > 65)
309 {
310 mpz_t q;
311 mpz_init (q);
312 mpz_tdiv_q_2exp (q, x, 4 * (mpz_sizeinbase (x, 16) - 25));
313 gmp_fprintf (stderr, "%ZX...", q);
314 mpz_tdiv_r_2exp (q, x, 4 * 25);
315 gmp_fprintf (stderr, "%025ZX [%d]\n", q, (int) siz);
316 mpz_clear (q);
317 }
318 else
319 {
320 gmp_fprintf (stderr, "%ZX\n", x);
321 }
322 }
323