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
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
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
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
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
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 }