1/* Test mpz_cmp, mpz_mul.
2
3Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004 Free Software Foundation,
4Inc.
5
6This file is part of the GNU MP Library test suite.
7
8The GNU MP Library test suite is free software; you can redistribute it
9and/or modify it under the terms of the GNU General Public License as
10published by the Free Software Foundation; either version 3 of the License,
11or (at your option) any later version.
12
13The GNU MP Library test suite is distributed in the hope that it will be
14useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
16Public License for more details.
17
18You should have received a copy of the GNU General Public License along with
19the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
20
21#include <stdio.h>
22#include <stdlib.h>
23
24#include "gmp-impl.h"
25#include "longlong.h"
26#include "tests.h"
27
28void debug_mp (mpz_t);
29static void refmpz_mul (mpz_t, const mpz_t, const mpz_t);
30void dump_abort (int, const char *, mpz_t, mpz_t, mpz_t, mpz_t);
31
32#define FFT_MIN_BITSIZE 100000
33
34char *extra_fft;
35
36void
37one (int i, mpz_t multiplicand, mpz_t multiplier)
38{
39  mpz_t product, ref_product;
40
41  mpz_init (product);
42  mpz_init (ref_product);
43
44  /* Test plain multiplication comparing results against reference code.  */
45  mpz_mul (product, multiplier, multiplicand);
46  refmpz_mul (ref_product, multiplier, multiplicand);
47  if (mpz_cmp (product, ref_product))
48    dump_abort (i, "incorrect plain product",
49		multiplier, multiplicand, product, ref_product);
50
51  /* Test squaring, comparing results against plain multiplication  */
52  mpz_mul (product, multiplier, multiplier);
53  mpz_set (multiplicand, multiplier);
54  mpz_mul (ref_product, multiplier, multiplicand);
55  if (mpz_cmp (product, ref_product))
56    dump_abort (i, "incorrect square product",
57		multiplier, multiplier, product, ref_product);
58
59  mpz_clear (product);
60  mpz_clear (ref_product);
61}
62
63int
64main (int argc, char **argv)
65{
66  mpz_t op1, op2;
67  int i;
68  int fft_max_2exp;
69
70  gmp_randstate_ptr rands;
71  mpz_t bs;
72  unsigned long bsi, size_range, fsize_range;
73
74  tests_start ();
75  rands = RANDS;
76
77  extra_fft = getenv ("GMP_CHECK_FFT");
78  fft_max_2exp = 0;
79  if (extra_fft != NULL)
80    {
81      fft_max_2exp = atoi (extra_fft);
82      printf ("GMP_CHECK_FFT=%d (include this in bug reports)\n", fft_max_2exp);
83    }
84
85  if (fft_max_2exp <= 1)	/* compat with old use of GMP_CHECK_FFT */
86    fft_max_2exp = 22;		/* default limit, good for any machine */
87
88  mpz_init (bs);
89  mpz_init (op1);
90  mpz_init (op2);
91
92  fsize_range = 4 << 8;		/* a fraction 1/256 of size_range */
93  for (i = 0;; i++)
94    {
95      size_range = fsize_range >> 8;
96      fsize_range = fsize_range * 33 / 32;
97
98      if (size_range > fft_max_2exp)
99	break;
100
101      mpz_urandomb (bs, rands, size_range);
102      mpz_rrandomb (op1, rands, mpz_get_ui (bs));
103      if (i & 1)
104	mpz_urandomb (bs, rands, size_range);
105      mpz_rrandomb (op2, rands, mpz_get_ui (bs));
106
107      mpz_urandomb (bs, rands, 4);
108      bsi = mpz_get_ui (bs);
109      if ((bsi & 0x3) == 0)
110	mpz_neg (op1, op1);
111      if ((bsi & 0xC) == 0)
112	mpz_neg (op2, op2);
113
114      /* printf ("%d %d\n", SIZ (op1), SIZ (op2)); */
115      one (i, op2, op1);
116    }
117
118  for (i = -50; i < 0; i++)
119    {
120      mpz_urandomb (bs, rands, 32);
121      size_range = mpz_get_ui (bs) % fft_max_2exp;
122
123      mpz_urandomb (bs, rands, size_range);
124      mpz_rrandomb (op1, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
125      mpz_urandomb (bs, rands, size_range);
126      mpz_rrandomb (op2, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
127
128      /* printf ("%d: %d %d\n", i, SIZ (op1), SIZ (op2)); */
129      fflush (stdout);
130      one (-1, op2, op1);
131    }
132
133  mpz_clear (bs);
134  mpz_clear (op1);
135  mpz_clear (op2);
136
137  tests_end ();
138  exit (0);
139}
140
141static void
142refmpz_mul (mpz_t w, const mpz_t u, const mpz_t v)
143{
144  mp_size_t usize = u->_mp_size;
145  mp_size_t vsize = v->_mp_size;
146  mp_size_t wsize;
147  mp_size_t sign_product;
148  mp_ptr up, vp;
149  mp_ptr wp;
150  mp_size_t talloc;
151
152  sign_product = usize ^ vsize;
153  usize = ABS (usize);
154  vsize = ABS (vsize);
155
156  if (usize == 0 || vsize == 0)
157    {
158      SIZ (w) = 0;
159      return;
160    }
161
162  talloc = usize + vsize;
163
164  up = u->_mp_d;
165  vp = v->_mp_d;
166
167  wp = __GMP_ALLOCATE_FUNC_LIMBS (talloc);
168
169  if (usize > vsize)
170    refmpn_mul (wp, up, usize, vp, vsize);
171  else
172    refmpn_mul (wp, vp, vsize, up, usize);
173  wsize = usize + vsize;
174  wsize -= wp[wsize - 1] == 0;
175  MPZ_REALLOC (w, wsize);
176  MPN_COPY (PTR(w), wp, wsize);
177
178  SIZ(w) = sign_product < 0 ? -wsize : wsize;
179  __GMP_FREE_FUNC_LIMBS (wp, talloc);
180}
181
182void
183dump_abort (int i, const char *s,
184            mpz_t op1, mpz_t op2, mpz_t product, mpz_t ref_product)
185{
186  mp_size_t b, e;
187  fprintf (stderr, "ERROR: %s in test %d\n", s, i);
188  fprintf (stderr, "op1          = "); debug_mp (op1);
189  fprintf (stderr, "op2          = "); debug_mp (op2);
190  fprintf (stderr, "    product  = "); debug_mp (product);
191  fprintf (stderr, "ref_product  = "); debug_mp (ref_product);
192  for (b = 0; b < ABSIZ(ref_product); b++)
193    if (PTR(ref_product)[b] != PTR(product)[b])
194      break;
195  for (e = ABSIZ(ref_product) - 1; e >= 0; e--)
196    if (PTR(ref_product)[e] != PTR(product)[e])
197      break;
198  printf ("ERRORS in %ld--%ld\n", b, e);
199  abort();
200}
201
202void
203debug_mp (mpz_t x)
204{
205  size_t siz = mpz_sizeinbase (x, 16);
206
207  if (siz > 65)
208    {
209      mpz_t q;
210      mpz_init (q);
211      mpz_tdiv_q_2exp (q, x, 4 * (mpz_sizeinbase (x, 16) - 25));
212      gmp_fprintf (stderr, "%ZX...", q);
213      mpz_tdiv_r_2exp (q, x, 4 * 25);
214      gmp_fprintf (stderr, "%025ZX [%d]\n", q, (int) siz);
215      mpz_clear (q);
216    }
217  else
218    {
219      gmp_fprintf (stderr, "%ZX\n", x);
220    }
221}
222