1 /* Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
2
3 This file is part of the GNU MP Library test suite.
4
5 The GNU MP Library test suite is free software; you can redistribute it
6 and/or modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 3 of the License,
8 or (at your option) any later version.
9
10 The GNU MP Library test suite is distributed in the hope that it will be
11 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
13 Public License for more details.
14
15 You should have received a copy of the GNU General Public License along with
16 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
17
18
19 #include <stdlib.h> /* for strtol */
20 #include <stdio.h> /* for printf */
21
22 #include "gmp.h"
23 #include "gmp-impl.h"
24 #include "longlong.h"
25 #include "tests/tests.h"
26
27
28 static void
dumpy(mp_srcptr p,mp_size_t n)29 dumpy (mp_srcptr p, mp_size_t n)
30 {
31 mp_size_t i;
32 if (n > 20)
33 {
34 for (i = n - 1; i >= n - 4; i--)
35 {
36 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
37 printf (" ");
38 }
39 printf ("... ");
40 for (i = 3; i >= 0; i--)
41 {
42 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
43 printf (" " + (i == 0));
44 }
45 }
46 else
47 {
48 for (i = n - 1; i >= 0; i--)
49 {
50 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
51 printf (" " + (i == 0));
52 }
53 }
54 puts ("");
55 }
56
57 static unsigned long test;
58
59 void
check_one(mp_ptr qp,mp_srcptr rp,mp_limb_t rh,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,const char * fname)60 check_one (mp_ptr qp, mp_srcptr rp, mp_limb_t rh,
61 mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, const char *fname)
62 {
63 mp_size_t qn;
64 int cmp;
65 mp_ptr tp;
66 mp_limb_t cy = 4711; /* silence warnings */
67 TMP_DECL;
68
69 qn = nn - dn;
70
71 if (qn == 0)
72 return;
73
74 TMP_MARK;
75
76 tp = TMP_ALLOC_LIMBS (nn + 1);
77
78 if (dn >= qn)
79 mpn_mul (tp, dp, dn, qp, qn);
80 else
81 mpn_mul (tp, qp, qn, dp, dn);
82
83 if (rp != NULL)
84 {
85 cy = mpn_add_n (tp + qn, tp + qn, rp, dn);
86 cmp = cy != rh || mpn_cmp (tp, np, nn) != 0;
87 }
88 else
89 cmp = mpn_cmp (tp, np, nn - dn) != 0;
90
91 if (cmp != 0)
92 {
93 printf ("\r*******************************************************************************\n");
94 printf ("%s inconsistent in test %lu\n", fname, test);
95 printf ("N= "); dumpy (np, nn);
96 printf ("D= "); dumpy (dp, dn);
97 printf ("Q= "); dumpy (qp, qn);
98 if (rp != NULL)
99 {
100 printf ("R= "); dumpy (rp, dn);
101 printf ("Rb= %d, Cy=%d\n", (int) cy, (int) rh);
102 }
103 printf ("T= "); dumpy (tp, nn);
104 printf ("nn = %ld, dn = %ld, qn = %ld", nn, dn, qn);
105 printf ("\n*******************************************************************************\n");
106 abort ();
107 }
108
109 TMP_FREE;
110 }
111
112
113 /* These are *bit* sizes. */
114 #define SIZE_LOG 16
115 #define MAX_DN (1L << SIZE_LOG)
116 #define MAX_NN (1L << (SIZE_LOG + 1))
117
118 #define COUNT 500
119
120 mp_limb_t
random_word(gmp_randstate_ptr rs)121 random_word (gmp_randstate_ptr rs)
122 {
123 mpz_t x;
124 mp_limb_t r;
125 TMP_DECL;
126 TMP_MARK;
127
128 MPZ_TMP_INIT (x, 2);
129 mpz_urandomb (x, rs, 32);
130 r = mpz_get_ui (x);
131 TMP_FREE;
132 return r;
133 }
134
135 int
main(int argc,char ** argv)136 main (int argc, char **argv)
137 {
138 gmp_randstate_ptr rands;
139 unsigned long maxnbits, maxdbits, nbits, dbits;
140 mpz_t n, d, tz;
141 mp_size_t maxnn, maxdn, nn, dn, clearn, i;
142 mp_ptr np, dp, qp, rp;
143 mp_limb_t rh;
144 mp_limb_t t;
145 mp_limb_t dinv;
146 int count = COUNT;
147 mp_ptr scratch;
148 mp_limb_t ran;
149 mp_size_t alloc, itch;
150 mp_limb_t rran0, rran1, qran0, qran1;
151 TMP_DECL;
152
153 if (argc > 1)
154 {
155 char *end;
156 count = strtol (argv[1], &end, 0);
157 if (*end || count <= 0)
158 {
159 fprintf (stderr, "Invalid test count: %s.\n", argv[1]);
160 return 1;
161 }
162 }
163
164
165 maxdbits = MAX_DN;
166 maxnbits = MAX_NN;
167
168 tests_start ();
169 rands = RANDS;
170
171 mpz_init (n);
172 mpz_init (d);
173 mpz_init (tz);
174
175 maxnn = maxnbits / GMP_NUMB_BITS + 1;
176 maxdn = maxdbits / GMP_NUMB_BITS + 1;
177
178 TMP_MARK;
179
180 qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
181 rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
182
183 alloc = 1;
184 scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);
185
186 for (test = 0; test < count;)
187 {
188 nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;
189 if (maxdbits > nbits)
190 dbits = random_word (rands) % nbits + 1;
191 else
192 dbits = random_word (rands) % maxdbits + 1;
193
194 #if RAND_UNIFORM
195 #define RANDFUNC mpz_urandomb
196 #else
197 #define RANDFUNC mpz_rrandomb
198 #endif
199
200 do
201 {
202 RANDFUNC (n, rands, nbits);
203 do
204 {
205 RANDFUNC (d, rands, dbits);
206 }
207 while (mpz_sgn (d) == 0);
208
209 np = PTR (n);
210 dp = PTR (d);
211 nn = SIZ (n);
212 dn = SIZ (d);
213 }
214 while (nn < dn);
215
216 dp[0] |= 1;
217
218 mpz_urandomb (tz, rands, 32);
219 t = mpz_get_ui (tz);
220
221 if (t % 17 == 0)
222 dp[0] = GMP_NUMB_MAX;
223
224 switch ((int) t % 16)
225 {
226 case 0:
227 clearn = random_word (rands) % nn;
228 for (i = 0; i <= clearn; i++)
229 np[i] = 0;
230 break;
231 case 1:
232 mpn_sub_1 (np + nn - dn, dp, dn, random_word (rands));
233 break;
234 case 2:
235 mpn_add_1 (np + nn - dn, dp, dn, random_word (rands));
236 break;
237 }
238
239 test++;
240
241 binvert_limb (dinv, dp[0]);
242
243 rran0 = random_word (rands);
244 rran1 = random_word (rands);
245 qran0 = random_word (rands);
246 qran1 = random_word (rands);
247
248 qp[-1] = qran0;
249 qp[nn - dn + 1] = qran1;
250 rp[-1] = rran0;
251
252 ran = random_word (rands);
253
254 if ((double) (nn - dn) * dn < 1e5)
255 {
256 if (nn > dn)
257 {
258 /* Test mpn_sbpi1_bdiv_qr */
259 MPN_ZERO (qp, nn - dn);
260 MPN_ZERO (rp, dn);
261 MPN_COPY (rp, np, nn);
262 rh = mpn_sbpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
263 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
264 ASSERT_ALWAYS (rp[-1] == rran0);
265 check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_sbpi1_bdiv_qr");
266 }
267
268 if (nn > dn)
269 {
270 /* Test mpn_sbpi1_bdiv_q */
271 MPN_COPY (rp, np, nn);
272 MPN_ZERO (qp, nn - dn);
273 mpn_sbpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
274 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
275 ASSERT_ALWAYS (rp[-1] == rran0);
276 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_sbpi1_bdiv_q");
277 }
278 }
279
280 if (dn >= 4 && nn - dn >= 2)
281 {
282 /* Test mpn_dcpi1_bdiv_qr */
283 MPN_COPY (rp, np, nn);
284 MPN_ZERO (qp, nn - dn);
285 rh = mpn_dcpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv);
286 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
287 ASSERT_ALWAYS (rp[-1] == rran0);
288 check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_dcpi1_bdiv_qr");
289 }
290
291 if (dn >= 4 && nn - dn >= 2)
292 {
293 /* Test mpn_dcpi1_bdiv_q */
294 MPN_COPY (rp, np, nn);
295 MPN_ZERO (qp, nn - dn);
296 mpn_dcpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv);
297 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
298 ASSERT_ALWAYS (rp[-1] == rran0);
299 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_dcpi1_bdiv_q");
300 }
301
302 if (nn > dn)
303 {
304 /* Test mpn_bdiv_qr */
305 itch = mpn_bdiv_qr_itch (nn, dn);
306 if (itch + 1 > alloc)
307 {
308 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
309 alloc = itch + 1;
310 }
311 scratch[itch] = ran;
312 MPN_ZERO (qp, nn - dn);
313 MPN_ZERO (rp, dn);
314 rp[dn] = rran1;
315 rh = mpn_bdiv_qr (qp, rp, np, nn, dp, dn, scratch);
316 ASSERT_ALWAYS (ran == scratch[itch]);
317 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
318 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1);
319
320 check_one (qp, rp, rh, np, nn, dp, dn, "mpn_bdiv_qr");
321 }
322
323 if (nn - dn < 2 || dn < 2)
324 continue;
325
326 /* Test mpn_mu_bdiv_qr */
327 itch = mpn_mu_bdiv_qr_itch (nn, dn);
328 if (itch + 1 > alloc)
329 {
330 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
331 alloc = itch + 1;
332 }
333 scratch[itch] = ran;
334 MPN_ZERO (qp, nn - dn);
335 MPN_ZERO (rp, dn);
336 rp[dn] = rran1;
337 rh = mpn_mu_bdiv_qr (qp, rp, np, nn, dp, dn, scratch);
338 ASSERT_ALWAYS (ran == scratch[itch]);
339 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
340 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1);
341 check_one (qp, rp, rh, np, nn, dp, dn, "mpn_mu_bdiv_qr");
342
343 /* Test mpn_mu_bdiv_q */
344 itch = mpn_mu_bdiv_q_itch (nn, dn);
345 if (itch + 1 > alloc)
346 {
347 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
348 alloc = itch + 1;
349 }
350 scratch[itch] = ran;
351 MPN_ZERO (qp, nn - dn + 1);
352 mpn_mu_bdiv_q (qp, np, nn - dn, dp, dn, scratch);
353 ASSERT_ALWAYS (ran == scratch[itch]);
354 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
355 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_mu_bdiv_q");
356 }
357
358 __GMP_FREE_FUNC_LIMBS (scratch, alloc);
359
360 TMP_FREE;
361
362 mpz_clear (n);
363 mpz_clear (d);
364 mpz_clear (tz);
365
366 tests_end ();
367 return 0;
368 }
369