1 /* Copyright 2006, 2007, 2009, 2010, 2013-2015, 2018 Free Software
2 Foundation, Inc.
3
4 This file is part of the GNU MP Library test suite.
5
6 The GNU MP Library test suite is free software; you can redistribute it
7 and/or modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 3 of the License,
9 or (at your option) any later version.
10
11 The GNU MP Library test suite is distributed in the hope that it will be
12 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
14 Public License for more details.
15
16 You should have received a copy of the GNU General Public License along with
17 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
18
19
20 #include <stdlib.h> /* for strtol */
21 #include <stdio.h> /* for printf */
22
23 #include "gmp-impl.h"
24 #include "longlong.h"
25 #include "tests/tests.h"
26
27
28 static void
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 signed long test;
58
59 static void
60 check_one (mp_ptr qp, mp_srcptr rp,
61 mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn,
62 const char *fname, mp_limb_t q_allowed_err)
63 {
64 mp_size_t qn = nn - dn + 1;
65 mp_ptr tp;
66 const char *msg;
67 const char *tvalue;
68 mp_limb_t i;
69 TMP_DECL;
70 TMP_MARK;
71
72 tp = TMP_ALLOC_LIMBS (nn + 1);
73 if (dn >= qn)
74 refmpn_mul (tp, dp, dn, qp, qn);
75 else
76 refmpn_mul (tp, qp, qn, dp, dn);
77
78 for (i = 0; i < q_allowed_err && (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0); i++)
79 ASSERT_NOCARRY (refmpn_sub (tp, tp, nn+1, dp, dn));
80
81 if (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0)
82 {
83 msg = "q too large";
84 tvalue = "Q*D";
85 error:
86 printf ("\r*******************************************************************************\n");
87 printf ("%s failed test %ld: %s\n", fname, test, msg);
88 printf ("N= "); dumpy (np, nn);
89 printf ("D= "); dumpy (dp, dn);
90 printf ("Q= "); dumpy (qp, qn);
91 if (rp)
92 { printf ("R= "); dumpy (rp, dn); }
93 printf ("%5s=", tvalue); dumpy (tp, nn+1);
94 printf ("nn = %ld, dn = %ld, qn = %ld\n", nn, dn, qn);
95 abort ();
96 }
97
98 ASSERT_NOCARRY (refmpn_sub_n (tp, np, tp, nn));
99 tvalue = "N-Q*D";
100 if (!(nn == dn || mpn_zero_p (tp + dn, nn - dn)) || mpn_cmp (tp, dp, dn) >= 0)
101 {
102 msg = "q too small";
103 goto error;
104 }
105
106 if (rp && mpn_cmp (rp, tp, dn) != 0)
107 {
108 msg = "r incorrect";
109 goto error;
110 }
111
112 TMP_FREE;
113 }
114
115
116 /* These are *bit* sizes. */
117 #ifndef SIZE_LOG
118 #define SIZE_LOG 17
119 #endif
120 #define MAX_DN (1L << SIZE_LOG)
121 #define MAX_NN (1L << (SIZE_LOG + 1))
122
123 #define COUNT 200
124
125 int
126 main (int argc, char **argv)
127 {
128 gmp_randstate_ptr rands;
129 unsigned long maxnbits, maxdbits, nbits, dbits;
130 mpz_t n, d, q, r, tz, junk;
131 mp_size_t maxnn, maxdn, nn, dn, clearn, i;
132 mp_ptr np, dup, dnp, qp, rp, junkp;
133 mp_limb_t t;
134 gmp_pi1_t dinv;
135 long count = COUNT;
136 mp_ptr scratch;
137 mp_limb_t ran;
138 mp_size_t alloc, itch;
139 mp_limb_t rran0, rran1, qran0, qran1;
140 TMP_DECL;
141
142 TESTS_REPS (count, argv, argc);
143
144 maxdbits = MAX_DN;
145 maxnbits = MAX_NN;
146
147 tests_start ();
148 rands = RANDS;
149
150 mpz_init (n);
151 mpz_init (d);
152 mpz_init (q);
153 mpz_init (r);
154 mpz_init (tz);
155 mpz_init (junk);
156
157 maxnn = maxnbits / GMP_NUMB_BITS + 1;
158 maxdn = maxdbits / GMP_NUMB_BITS + 1;
159
160 TMP_MARK;
161
162 qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
163 rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
164 dnp = TMP_ALLOC_LIMBS (maxdn);
165
166 alloc = 1;
167 scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);
168
169 for (test = -300; test < count; test++)
170 {
171 nbits = urandom () % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;
172
173 if (test < 0)
174 dbits = (test + 300) % (nbits - 1) + 1;
175 else
176 dbits = urandom () % (nbits - 1) % maxdbits + 1;
177
178 #if RAND_UNIFORM
179 #define RANDFUNC mpz_urandomb
180 #else
181 #define RANDFUNC mpz_rrandomb
182 #endif
183
184 do
185 RANDFUNC (d, rands, dbits);
186 while (mpz_sgn (d) == 0);
187 dn = SIZ (d);
188 dup = PTR (d);
189 MPN_COPY (dnp, dup, dn);
190 dnp[dn - 1] |= GMP_NUMB_HIGHBIT;
191
192 if (test % 2 == 0)
193 {
194 RANDFUNC (n, rands, nbits);
195 nn = SIZ (n);
196 ASSERT_ALWAYS (nn >= dn);
197 }
198 else
199 {
200 do
201 {
202 RANDFUNC (q, rands, urandom () % (nbits - dbits + 1));
203 RANDFUNC (r, rands, urandom () % mpz_sizeinbase (d, 2));
204 mpz_mul (n, q, d);
205 mpz_add (n, n, r);
206 nn = SIZ (n);
207 }
208 while (nn > maxnn || nn < dn);
209 }
210
211 ASSERT_ALWAYS (nn <= maxnn);
212 ASSERT_ALWAYS (dn <= maxdn);
213
214 mpz_urandomb (junk, rands, nbits);
215 junkp = PTR (junk);
216
217 np = PTR (n);
218
219 mpz_urandomb (tz, rands, 32);
220 t = mpz_get_ui (tz);
221
222 if (t % 17 == 0)
223 {
224 dnp[dn - 1] = GMP_NUMB_MAX;
225 dup[dn - 1] = GMP_NUMB_MAX;
226 }
227
228 switch ((int) t % 16)
229 {
230 case 0:
231 clearn = urandom () % nn;
232 for (i = clearn; i < nn; i++)
233 np[i] = 0;
234 break;
235 case 1:
236 mpn_sub_1 (np + nn - dn, dnp, dn, urandom ());
237 break;
238 case 2:
239 mpn_add_1 (np + nn - dn, dnp, dn, urandom ());
240 break;
241 }
242
243 if (dn >= 2)
244 invert_pi1 (dinv, dnp[dn - 1], dnp[dn - 2]);
245
246 rran0 = urandom ();
247 rran1 = urandom ();
248 qran0 = urandom ();
249 qran1 = urandom ();
250
251 qp[-1] = qran0;
252 qp[nn - dn + 1] = qran1;
253 rp[-1] = rran0;
254
255 ran = urandom ();
256
257 if ((double) (nn - dn) * dn < 1e5)
258 {
259 /* Test mpn_sbpi1_div_qr */
260 if (dn > 2)
261 {
262 MPN_COPY (rp, np, nn);
263 if (nn > dn)
264 MPN_COPY (qp, junkp, nn - dn);
265 qp[nn - dn] = mpn_sbpi1_div_qr (qp, rp, nn, dnp, dn, dinv.inv32);
266 check_one (qp, rp, np, nn, dnp, dn, "mpn_sbpi1_div_qr", 0);
267 }
268
269 /* Test mpn_sbpi1_divappr_q */
270 if (dn > 2)
271 {
272 MPN_COPY (rp, np, nn);
273 if (nn > dn)
274 MPN_COPY (qp, junkp, nn - dn);
275 qp[nn - dn] = mpn_sbpi1_divappr_q (qp, rp, nn, dnp, dn, dinv.inv32);
276 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_divappr_q", 1);
277 }
278
279 /* Test mpn_sbpi1_div_q */
280 if (dn > 2)
281 {
282 MPN_COPY (rp, np, nn);
283 if (nn > dn)
284 MPN_COPY (qp, junkp, nn - dn);
285 qp[nn - dn] = mpn_sbpi1_div_q (qp, rp, nn, dnp, dn, dinv.inv32);
286 check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_div_q", 0);
287 }
288
289 /* Test mpn_sec_div_qr */
290 itch = mpn_sec_div_qr_itch (nn, dn);
291 if (itch + 1 > alloc)
292 {
293 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
294 alloc = itch + 1;
295 }
296 scratch[itch] = ran;
297 MPN_COPY (rp, np, nn);
298 if (nn >= dn)
299 MPN_COPY (qp, junkp, nn - dn + 1);
300 qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dup, dn, scratch);
301 ASSERT_ALWAYS (ran == scratch[itch]);
302 check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_qr (unnorm)", 0);
303
304 /* Test mpn_sec_div_r */
305 itch = mpn_sec_div_r_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_COPY (rp, np, nn);
313 mpn_sec_div_r (rp, nn, dup, dn, scratch);
314 ASSERT_ALWAYS (ran == scratch[itch]);
315 /* Note: Since check_one cannot cope with remainder-only functions, we
316 pass qp[] from the previous function, mpn_sec_div_qr. */
317 check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_r (unnorm)", 0);
318
319 /* Normalised case, mpn_sec_div_qr */
320 itch = mpn_sec_div_qr_itch (nn, dn);
321 scratch[itch] = ran;
322
323 MPN_COPY (rp, np, nn);
324 if (nn >= dn)
325 MPN_COPY (qp, junkp, nn - dn + 1);
326 qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dnp, dn, scratch);
327 ASSERT_ALWAYS (ran == scratch[itch]);
328 check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_qr (norm)", 0);
329
330 /* Normalised case, mpn_sec_div_r */
331 itch = mpn_sec_div_r_itch (nn, dn);
332 scratch[itch] = ran;
333 MPN_COPY (rp, np, nn);
334 mpn_sec_div_r (rp, nn, dnp, dn, scratch);
335 ASSERT_ALWAYS (ran == scratch[itch]);
336 /* Note: Since check_one cannot cope with remainder-only functions, we
337 pass qp[] from the previous function, mpn_sec_div_qr. */
338 check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_r (norm)", 0);
339 }
340
341 /* Test mpn_dcpi1_div_qr */
342 if (dn >= 6 && nn - dn >= 3)
343 {
344 MPN_COPY (rp, np, nn);
345 if (nn > dn)
346 MPN_COPY (qp, junkp, nn - dn);
347 qp[nn - dn] = mpn_dcpi1_div_qr (qp, rp, nn, dnp, dn, &dinv);
348 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
349 ASSERT_ALWAYS (rp[-1] == rran0);
350 check_one (qp, rp, np, nn, dnp, dn, "mpn_dcpi1_div_qr", 0);
351 }
352
353 /* Test mpn_dcpi1_divappr_q */
354 if (dn >= 6 && nn - dn >= 3)
355 {
356 MPN_COPY (rp, np, nn);
357 if (nn > dn)
358 MPN_COPY (qp, junkp, nn - dn);
359 qp[nn - dn] = mpn_dcpi1_divappr_q (qp, rp, nn, dnp, dn, &dinv);
360 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
361 ASSERT_ALWAYS (rp[-1] == rran0);
362 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_divappr_q", 1);
363 }
364
365 /* Test mpn_dcpi1_div_q */
366 if (dn >= 6 && nn - dn >= 3)
367 {
368 MPN_COPY (rp, np, nn);
369 if (nn > dn)
370 MPN_COPY (qp, junkp, nn - dn);
371 qp[nn - dn] = mpn_dcpi1_div_q (qp, rp, nn, dnp, dn, &dinv);
372 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
373 ASSERT_ALWAYS (rp[-1] == rran0);
374 check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_div_q", 0);
375 }
376
377 /* Test mpn_mu_div_qr */
378 if (nn - dn > 2 && dn >= 2)
379 {
380 itch = mpn_mu_div_qr_itch (nn, dn, 0);
381 if (itch + 1 > alloc)
382 {
383 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
384 alloc = itch + 1;
385 }
386 scratch[itch] = ran;
387 MPN_COPY (qp, junkp, nn - dn);
388 MPN_ZERO (rp, dn);
389 rp[dn] = rran1;
390 qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dnp, dn, scratch);
391 ASSERT_ALWAYS (ran == scratch[itch]);
392 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
393 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1);
394 check_one (qp, rp, np, nn, dnp, dn, "mpn_mu_div_qr", 0);
395 }
396
397 /* Test mpn_mu_divappr_q */
398 if (nn - dn > 2 && dn >= 2)
399 {
400 itch = mpn_mu_divappr_q_itch (nn, dn, 0);
401 if (itch + 1 > alloc)
402 {
403 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
404 alloc = itch + 1;
405 }
406 scratch[itch] = ran;
407 MPN_COPY (qp, junkp, nn - dn);
408 qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dnp, dn, scratch);
409 ASSERT_ALWAYS (ran == scratch[itch]);
410 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
411 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_divappr_q", 4);
412 }
413
414 /* Test mpn_mu_div_q */
415 if (nn - dn > 2 && dn >= 2)
416 {
417 itch = mpn_mu_div_q_itch (nn, dn, 0);
418 if (itch + 1> alloc)
419 {
420 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
421 alloc = itch + 1;
422 }
423 scratch[itch] = ran;
424 MPN_COPY (qp, junkp, nn - dn);
425 qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dnp, dn, scratch);
426 ASSERT_ALWAYS (ran == scratch[itch]);
427 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
428 check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_div_q", 0);
429 }
430
431 if (1)
432 {
433 itch = nn + 1;
434 if (itch + 1> alloc)
435 {
436 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
437 alloc = itch + 1;
438 }
439 scratch[itch] = ran;
440 MPN_COPY (qp, junkp, nn - dn + 1);
441 mpn_div_q (qp, np, nn, dup, dn, scratch);
442 ASSERT_ALWAYS (ran == scratch[itch]);
443 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
444 check_one (qp, NULL, np, nn, dup, dn, "mpn_div_q", 0);
445 }
446
447 if (dn >= 2 && nn >= 2)
448 {
449 mp_limb_t qh;
450
451 /* mpn_divrem_2 */
452 MPN_COPY (rp, np, nn);
453 qp[nn - 2] = qp[nn-1] = qran1;
454
455 qh = mpn_divrem_2 (qp, 0, rp, nn, dnp + dn - 2);
456 ASSERT_ALWAYS (qp[nn - 2] == qran1);
457 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
458 qp[nn - 2] = qh;
459 check_one (qp, rp, np, nn, dnp + dn - 2, 2, "mpn_divrem_2", 0);
460
461 /* Missing: divrem_2 with fraction limbs. */
462
463 /* mpn_div_qr_2 */
464 qp[nn - 2] = qran1;
465
466 qh = mpn_div_qr_2 (qp, rp, np, nn, dup + dn - 2);
467 ASSERT_ALWAYS (qp[nn - 2] == qran1);
468 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
469 qp[nn - 2] = qh;
470 check_one (qp, rp, np, nn, dup + dn - 2, 2, "mpn_div_qr_2", 0);
471 }
472 if (dn >= 1 && nn >= 1)
473 {
474 /* mpn_div_qr_1 */
475 mp_limb_t qh;
476 qp[nn-1] = qran1;
477 rp[0] = mpn_div_qr_1 (qp, &qh, np, nn, dnp[dn - 1]);
478 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
479 qp[nn - 1] = qh;
480 check_one (qp, rp, np, nn, dnp + dn - 1, 1, "mpn_div_qr_1", 0);
481 }
482 }
483
484 __GMP_FREE_FUNC_LIMBS (scratch, alloc);
485
486 TMP_FREE;
487
488 mpz_clear (n);
489 mpz_clear (d);
490 mpz_clear (q);
491 mpz_clear (r);
492 mpz_clear (tz);
493 mpz_clear (junk);
494
495 tests_end ();
496 return 0;
497 }