1 /* Reference mpz functions.
2
3 Copyright 1997, 1999-2002 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the 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. */
35 void
36 refmpz_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
45 unsigned long
46 refmpz_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
81 void
82 refmpz_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
138 int
139 refmpz_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 */
229 int
230 refmpz_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. */
239 int
240 refmpz_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
275 int
276 refmpz_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
286 int
287 refmpz_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
297 int
298 refmpz_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
308 int
309 refmpz_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
320 void
321 refmpz_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 }