1/* Test mpn_mul function for all sizes up to a selected limit.
2
3Copyright 2011, 2012 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
19
20
21#include <stdlib.h>
22#include <stdio.h>
23
24#include "gmp.h"
25#include "gmp-impl.h"
26#include "tests.h"
27
28static unsigned
29isqrt (unsigned t)
30{
31  unsigned s, b;
32
33  for (b = 0, s = t;  b++, s >>= 1; )
34    ;
35
36  s = 1 << (b >> 1);
37  if (b & 1)
38    s += s >> 1;
39
40  do
41    {
42      b = t / s;
43      s = (s + b) >> 1;
44    }
45  while (b < s);
46
47  return s;
48}
49
50int
51main (int argc, char **argv)
52{
53  mp_ptr ap, bp, rp, refp;
54  mp_size_t max_n, an, bn, rn;
55  gmp_randstate_ptr rands;
56  int reps;
57  TMP_DECL;
58  TMP_MARK;
59
60  reps = 1;
61
62  tests_start ();
63  TESTS_REPS (reps, argv, argc);
64
65  rands = RANDS;
66
67  /* Re-interpret reps argument as a size argument.  */
68  max_n = isqrt (reps * 25000);
69
70  ap = TMP_ALLOC_LIMBS (max_n + 1);
71  bp = TMP_ALLOC_LIMBS (max_n + 1);
72  rp = TMP_ALLOC_LIMBS (2 * max_n);
73  refp = TMP_ALLOC_LIMBS (2 * max_n);
74
75  for (an = 1; an <= max_n; an += 1)
76    {
77      for (bn = 1; bn <= an; bn += 1)
78	{
79	  mpn_random2 (ap, an + 1);
80	  mpn_random2 (bp, bn + 1);
81
82	  refmpn_mul (refp, ap, an, bp, bn);
83	  mpn_mul (rp, ap, an, bp, bn);
84
85	  rn = an + bn;
86	  if (mpn_cmp (refp, rp, rn))
87	    {
88	      printf ("ERROR, an = %d, bn = %d, rn = %d\n",
89		      (int) an, (int) bn, (int) rn);
90	      printf ("a: "); mpn_dump (ap, an);
91	      printf ("b: "); mpn_dump (bp, bn);
92	      printf ("r:   "); mpn_dump (rp, rn);
93	      printf ("ref: "); mpn_dump (refp, rn);
94	      abort();
95	    }
96	}
97    }
98  TMP_FREE;
99  tests_end ();
100  return 0;
101}
102