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