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