refmpz.c revision 1.1.1.4 
1/* Reference mpz functions.
2
3Copyright 1997, 1999-2002 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19
20/* always do assertion checking */
21#define WANT_ASSERT  1
22
23#include <stdio.h>
24#include <stdlib.h> /* for free */
25#include "gmp-impl.h"
26#include "longlong.h"
27#include "tests.h"
28
29
30/* Change this to "#define TRACE(x) x" for some traces. */
31#define TRACE(x)
32
33
34/* FIXME: Shouldn't use plain mpz functions in a reference routine. */
35void
36refmpz_combit (mpz_ptr r, unsigned long bit)
37{
38  if (mpz_tstbit (r, bit))
39    mpz_clrbit (r, bit);
40  else
41    mpz_setbit (r, bit);
42}
43
44
45unsigned long
46refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
47{
48  mp_size_t      xsize, ysize, tsize;
49  mp_ptr         xp, yp;
50  unsigned long  ret;
51
52  if ((SIZ(x) < 0 && SIZ(y) >= 0)
53      || (SIZ(y) < 0 && SIZ(x) >= 0))
54    return ULONG_MAX;
55
56  xsize = ABSIZ(x);
57  ysize = ABSIZ(y);
58  tsize = MAX (xsize, ysize);
59
60  xp = refmpn_malloc_limbs (tsize);
61  refmpn_zero (xp, tsize);
62  refmpn_copy (xp, PTR(x), xsize);
63
64  yp = refmpn_malloc_limbs (tsize);
65  refmpn_zero (yp, tsize);
66  refmpn_copy (yp, PTR(y), ysize);
67
68  if (SIZ(x) < 0)
69    refmpn_neg (xp, xp, tsize);
70
71  if (SIZ(x) < 0)
72    refmpn_neg (yp, yp, tsize);
73
74  ret = refmpn_hamdist (xp, yp, tsize);
75
76  free (xp);
77  free (yp);
78  return ret;
79}
80
81void
82refmpz_gcd (mpz_ptr g, mpz_srcptr a_orig, mpz_srcptr b_orig)
83{
84  mp_bitcnt_t a_twos, b_twos, common_twos;
85  mpz_t a;
86  mpz_t b;
87  mpz_init (a);
88  mpz_init (b);
89  mpz_abs (a, a_orig);
90  mpz_abs (b, b_orig);
91
92  if (mpz_sgn (a) == 0)
93    {
94      mpz_set (g, b);
95      return;
96    }
97  if (mpz_sgn (b) == 0)
98    {
99      mpz_set (g, a);
100      return;
101    }
102  a_twos = mpz_scan1 (a, 0);
103  mpz_tdiv_q_2exp (a, a, a_twos);
104
105  b_twos = mpz_scan1 (b, 0);
106  mpz_tdiv_q_2exp (b, b, b_twos);
107
108  common_twos = MIN(a_twos, b_twos);
109  for (;;)
110    {
111      int c;
112      mp_bitcnt_t twos;
113      c = mpz_cmp (a, b);
114      if (c == 0)
115	break;
116      if (c < 0)
117	mpz_swap (a, b);
118      mpz_sub (a, a, b);
119      twos = mpz_scan1 (a, 0);
120      mpz_tdiv_q_2exp (a, a, twos);
121    }
122  mpz_mul_2exp (g, a, common_twos);
123
124  mpz_clear (a);
125  mpz_clear (b);
126}
127
128/* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
129#define JACOBI_0Z(b)  JACOBI_0LS (PTR(b)[0], SIZ(b))
130
131/* (a/b) effect due to sign of b: mpz/mpz */
132#define JACOBI_BSGN_ZZ_BIT1(a, b)   JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b))
133
134/* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd;
135   is (-1/b) if a<0, or +1 if a>=0 */
136#define JACOBI_ASGN_ZZU_BIT1(a, b)  JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0])
137
138int
139refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig)
140{
141  unsigned long  twos;
142  mpz_t  a, b;
143  int    result_bit1 = 0;
144
145  if (mpz_sgn (b_orig) == 0)
146    return JACOBI_Z0 (a_orig);  /* (a/0) */
147
148  if (mpz_sgn (a_orig) == 0)
149    return JACOBI_0Z (b_orig);  /* (0/b) */
150
151  if (mpz_even_p (a_orig) && mpz_even_p (b_orig))
152    return 0;
153
154  if (mpz_cmp_ui (b_orig, 1) == 0)
155    return 1;
156
157  mpz_init_set (a, a_orig);
158  mpz_init_set (b, b_orig);
159
160  if (mpz_sgn (b) < 0)
161    {
162      result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b);
163      mpz_neg (b, b);
164    }
165  if (mpz_even_p (b))
166    {
167      twos = mpz_scan1 (b, 0L);
168      mpz_tdiv_q_2exp (b, b, twos);
169      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]);
170    }
171
172  if (mpz_sgn (a) < 0)
173    {
174      result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]);
175      mpz_neg (a, a);
176    }
177  if (mpz_even_p (a))
178    {
179      twos = mpz_scan1 (a, 0L);
180      mpz_tdiv_q_2exp (a, a, twos);
181      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
182    }
183
184  for (;;)
185    {
186      ASSERT (mpz_odd_p (a));
187      ASSERT (mpz_odd_p (b));
188      ASSERT (mpz_sgn (a) > 0);
189      ASSERT (mpz_sgn (b) > 0);
190
191      TRACE (printf ("top\n");
192	     mpz_trace (" a", a);
193	     mpz_trace (" b", b));
194
195      if (mpz_cmp (a, b) < 0)
196	{
197	  TRACE (printf ("swap\n"));
198	  mpz_swap (a, b);
199	  result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]);
200	}
201
202      if (mpz_cmp_ui (b, 1) == 0)
203	break;
204
205      mpz_sub (a, a, b);
206      TRACE (printf ("sub\n");
207	     mpz_trace (" a", a));
208      if (mpz_sgn (a) == 0)
209	goto zero;
210
211      twos = mpz_scan1 (a, 0L);
212      mpz_fdiv_q_2exp (a, a, twos);
213      TRACE (printf ("twos %lu\n", twos);
214	     mpz_trace (" a", a));
215      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
216    }
217
218  mpz_clear (a);
219  mpz_clear (b);
220  return JACOBI_BIT1_TO_PN (result_bit1);
221
222 zero:
223  mpz_clear (a);
224  mpz_clear (b);
225  return 0;
226}
227
228/* Same as mpz_kronecker, but ignoring factors of 2 on b */
229int
230refmpz_jacobi (mpz_srcptr a, mpz_srcptr b)
231{
232  ASSERT_ALWAYS (mpz_sgn (b) > 0);
233  ASSERT_ALWAYS (mpz_odd_p (b));
234
235  return refmpz_kronecker (a, b);
236}
237
238/* Legendre symbol via powm. p must be an odd prime. */
239int
240refmpz_legendre (mpz_srcptr a, mpz_srcptr p)
241{
242  int res;
243
244  mpz_t r;
245  mpz_t e;
246
247  ASSERT_ALWAYS (mpz_sgn (p) > 0);
248  ASSERT_ALWAYS (mpz_odd_p (p));
249
250  mpz_init (r);
251  mpz_init (e);
252
253  mpz_fdiv_r (r, a, p);
254
255  mpz_set (e, p);
256  mpz_sub_ui (e, e, 1);
257  mpz_fdiv_q_2exp (e, e, 1);
258  mpz_powm (r, r, e, p);
259
260  /* Normalize to a more or less symmetric range around zero */
261  if (mpz_cmp (r, e) > 0)
262    mpz_sub (r, r, p);
263
264  ASSERT_ALWAYS (mpz_cmpabs_ui (r, 1) <= 0);
265
266  res = mpz_sgn (r);
267
268  mpz_clear (r);
269  mpz_clear (e);
270
271  return res;
272}
273
274
275int
276refmpz_kronecker_ui (mpz_srcptr a, unsigned long b)
277{
278  mpz_t  bz;
279  int    ret;
280  mpz_init_set_ui (bz, b);
281  ret = refmpz_kronecker (a, bz);
282  mpz_clear (bz);
283  return ret;
284}
285
286int
287refmpz_kronecker_si (mpz_srcptr a, long b)
288{
289  mpz_t  bz;
290  int    ret;
291  mpz_init_set_si (bz, b);
292  ret = refmpz_kronecker (a, bz);
293  mpz_clear (bz);
294  return ret;
295}
296
297int
298refmpz_ui_kronecker (unsigned long a, mpz_srcptr b)
299{
300  mpz_t  az;
301  int    ret;
302  mpz_init_set_ui (az, a);
303  ret = refmpz_kronecker (az, b);
304  mpz_clear (az);
305  return ret;
306}
307
308int
309refmpz_si_kronecker (long a, mpz_srcptr b)
310{
311  mpz_t  az;
312  int    ret;
313  mpz_init_set_si (az, a);
314  ret = refmpz_kronecker (az, b);
315  mpz_clear (az);
316  return ret;
317}
318
319
320void
321refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e)
322{
323  mpz_t          s, t;
324  unsigned long  i;
325
326  mpz_init_set_ui (t, 1L);
327  mpz_init_set (s, b);
328
329  if ((e & 1) != 0)
330    mpz_mul (t, t, s);
331
332  for (i = 2; i <= e; i <<= 1)
333    {
334      mpz_mul (s, s, s);
335      if ((i & e) != 0)
336	mpz_mul (t, t, s);
337    }
338
339  mpz_set (w, t);
340
341  mpz_clear (s);
342  mpz_clear (t);
343}
344