1 /* Reference mpz functions.
3 Copyright 1997, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
20 /* always do assertion checking */
24 #include <stdlib.h> /* for free */
31 /* Change this to "#define TRACE(x) x" for some traces. */
35 /* FIXME: Shouldn't use plain mpz functions in a reference routine. */
37 refmpz_combit (mpz_ptr r, unsigned long bit)
39 if (mpz_tstbit (r, bit))
47 refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
49 mp_size_t xsize, ysize, tsize;
53 if ((SIZ(x) < 0 && SIZ(y) >= 0)
54 || (SIZ(y) < 0 && SIZ(x) >= 0))
59 tsize = MAX (xsize, ysize);
61 xp = refmpn_malloc_limbs (tsize);
62 refmpn_zero (xp, tsize);
63 refmpn_copy (xp, PTR(x), xsize);
65 yp = refmpn_malloc_limbs (tsize);
66 refmpn_zero (yp, tsize);
67 refmpn_copy (yp, PTR(y), ysize);
70 refmpn_neg (xp, xp, tsize);
73 refmpn_neg (yp, yp, tsize);
75 ret = refmpn_hamdist (xp, yp, tsize);
83 /* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
84 #define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b))
86 /* (a/b) effect due to sign of b: mpz/mpz */
87 #define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b))
89 /* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd;
90 is (-1/b) if a<0, or +1 if a>=0 */
91 #define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0])
94 refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig)
100 if (mpz_sgn (b_orig) == 0)
101 return JACOBI_Z0 (a_orig); /* (a/0) */
103 if (mpz_sgn (a_orig) == 0)
104 return JACOBI_0Z (b_orig); /* (0/b) */
106 if (mpz_even_p (a_orig) && mpz_even_p (b_orig))
109 if (mpz_cmp_ui (b_orig, 1) == 0)
112 mpz_init_set (a, a_orig);
113 mpz_init_set (b, b_orig);
117 result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b);
122 twos = mpz_scan1 (b, 0L);
123 mpz_tdiv_q_2exp (b, b, twos);
124 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]);
129 result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]);
134 twos = mpz_scan1 (a, 0L);
135 mpz_tdiv_q_2exp (a, a, twos);
136 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
141 ASSERT (mpz_odd_p (a));
142 ASSERT (mpz_odd_p (b));
143 ASSERT (mpz_sgn (a) > 0);
144 ASSERT (mpz_sgn (b) > 0);
146 TRACE (printf ("top\n");
148 mpz_trace (" b", b));
150 if (mpz_cmp (a, b) < 0)
152 TRACE (printf ("swap\n"));
154 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]);
157 if (mpz_cmp_ui (b, 1) == 0)
161 TRACE (printf ("sub\n");
162 mpz_trace (" a", a));
163 if (mpz_sgn (a) == 0)
166 twos = mpz_scan1 (a, 0L);
167 mpz_fdiv_q_2exp (a, a, twos);
168 TRACE (printf ("twos %lu\n", twos);
169 mpz_trace (" a", a));
170 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
175 return JACOBI_BIT1_TO_PN (result_bit1);
183 /* Same as mpz_kronecker, but ignoring factors of 2 on b */
185 refmpz_jacobi (mpz_srcptr a, mpz_srcptr b)
188 mpz_init_set (b_odd, b);
189 if (mpz_sgn (b_odd) != 0)
190 mpz_fdiv_q_2exp (b_odd, b_odd, mpz_scan1 (b_odd, 0L));
191 return refmpz_kronecker (a, b_odd);
195 refmpz_legendre (mpz_srcptr a, mpz_srcptr b)
197 return refmpz_jacobi (a, b);
202 refmpz_kronecker_ui (mpz_srcptr a, unsigned long b)
206 mpz_init_set_ui (bz, b);
207 ret = refmpz_kronecker (a, bz);
213 refmpz_kronecker_si (mpz_srcptr a, long b)
217 mpz_init_set_si (bz, b);
218 ret = refmpz_kronecker (a, bz);
224 refmpz_ui_kronecker (unsigned long a, mpz_srcptr b)
228 mpz_init_set_ui (az, a);
229 ret = refmpz_kronecker (az, b);
235 refmpz_si_kronecker (long a, mpz_srcptr b)
239 mpz_init_set_si (az, a);
240 ret = refmpz_kronecker (az, b);
247 refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e)
252 mpz_init_set_ui (t, 1L);
258 for (i = 2; i <= e; i <<= 1)