xref: /netbsd-src/external/lgpl3/gmp/dist/tests/mpz/t-mul.c (revision 154bfe8e089c1a0a4e9ed8414f08d3da90949162)
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.h"
25 #include "gmp-impl.h"
26 #include "longlong.h"
27 #include "tests.h"
28 
29 void debug_mp (mpz_t);
30 static void refmpz_mul (mpz_t, const mpz_t, const mpz_t);
31 void dump_abort (int, const 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, const 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 }
220