1 /* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
2
3 Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library.
6
7 The GNU MP Library is free software; you can redistribute it and/or modify it
8 under the terms of the GNU Lesser General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
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 or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
15 for more details.
16
17 You should have received a copy of the GNU Lesser General Public License along
18 with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
19
20 #include <stdio.h>
21 #include "gmp.h"
22 #include "gmp-impl.h"
23 #include "longlong.h"
24
25
26 /* Change this to "#define TRACE(x) x" for some traces. */
27 #define TRACE(x)
28
29
30 #define MPN_RSHIFT_OR_COPY(dst,src,size,shift) \
31 do { \
32 if ((shift) != 0) \
33 { \
34 ASSERT_NOCARRY (mpn_rshift (dst, src, size, shift)); \
35 (size) -= ((dst)[(size)-1] == 0); \
36 } \
37 else \
38 MPN_COPY (dst, src, size); \
39 } while (0)
40
41
42 /* This code does triple duty as mpz_jacobi, mpz_legendre and mpz_kronecker.
43
44 mpz_jacobi could assume b is odd, but the improvements from that seem
45 small compared to other operations, and anything significant should be
46 checked at run-time since we'd like odd b to go fast in mpz_kronecker
47 too.
48
49 mpz_legendre could assume b is an odd prime, but knowing this doesn't
50 present any obvious benefits. Result 0 wouldn't arise (unless "a" is a
51 multiple of b), but the checking for that takes little time compared to
52 other operations.
53
54 The main loop is just a simple binary GCD with the jacobi symbol result
55 tracked during the reduction.
56
57 The special cases for a or b fitting in one limb let mod_1 or modexact_1
58 get used, without any copying, and end up just as efficient as the mixed
59 precision mpz_kronecker_ui etc.
60
61 When tdiv_qr is called it's not necessary to make "a" odd or make a
62 working copy of it, but tdiv_qr is going to be pretty slow so it's not
63 worth bothering trying to save anything for that case.
64
65 Enhancements:
66
67 mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
68
69 Some sort of multi-step algorithm should be used. The current subtract
70 and shift for every bit is very inefficient. Lehmer (per current gcdext)
71 would need some low bits included in its calculation to apply the sign
72 change for reciprocity. Binary Lehmer keeps low bits to strip twos
73 anyway, so might be better suited. Maybe the accelerated GCD style k-ary
74 reduction would work, if sign changes due to the extra factors it
75 introduces can be accounted for (or maybe they can be ignored). */
76
77
78 int
mpz_jacobi(mpz_srcptr a,mpz_srcptr b)79 mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
80 {
81 mp_srcptr asrcp, bsrcp;
82 mp_size_t asize, bsize;
83 mp_ptr ap, bp;
84 mp_limb_t alow, blow, ahigh, bhigh, asecond, bsecond;
85 unsigned atwos, btwos;
86 int result_bit1;
87 TMP_DECL;
88
89 TRACE (printf ("start asize=%d bsize=%d\n", SIZ(a), SIZ(b));
90 mpz_trace (" a", a);
91 mpz_trace (" b", b));
92
93 asize = SIZ(a);
94 asrcp = PTR(a);
95 alow = asrcp[0];
96
97 bsize = SIZ(b);
98 if (bsize == 0)
99 return JACOBI_LS0 (alow, asize); /* (a/0) */
100
101 bsrcp = PTR(b);
102 blow = bsrcp[0];
103
104 if (asize == 0)
105 return JACOBI_0LS (blow, bsize); /* (0/b) */
106
107 /* (even/even)=0 */
108 if (((alow | blow) & 1) == 0)
109 return 0;
110
111 /* account for effect of sign of b, then ignore it */
112 result_bit1 = JACOBI_BSGN_SS_BIT1 (asize, bsize);
113 bsize = ABS (bsize);
114
115 /* low zero limbs on b can be discarded */
116 JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
117
118 count_trailing_zeros (btwos, blow);
119 TRACE (printf ("b twos %u\n", btwos));
120
121 /* establish shifted blow */
122 blow >>= btwos;
123 if (bsize > 1)
124 {
125 bsecond = bsrcp[1];
126 if (btwos != 0)
127 blow |= (bsecond << (GMP_NUMB_BITS - btwos)) & GMP_NUMB_MASK;
128 }
129
130 /* account for effect of sign of a, then ignore it */
131 result_bit1 ^= JACOBI_ASGN_SU_BIT1 (asize, blow);
132 asize = ABS (asize);
133
134 if (bsize == 1 || (bsize == 2 && (bsecond >> btwos) == 0))
135 {
136 /* special case one limb b, use modexact and no copying */
137
138 /* (a/2)=(2/a) with a odd, and if b is even then a is odd here */
139 result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
140
141 if (blow == 1) /* (a/1)=1 always */
142 return JACOBI_BIT1_TO_PN (result_bit1);
143
144 JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
145 TRACE (printf ("base (%lu/%lu) with %d\n",
146 alow, blow, JACOBI_BIT1_TO_PN (result_bit1)));
147 return mpn_jacobi_base (alow, blow, result_bit1);
148 }
149
150 /* Discard low zero limbs of a. Usually there won't be anything to
151 strip, hence not bothering with it for the bsize==1 case. */
152 JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
153
154 count_trailing_zeros (atwos, alow);
155 TRACE (printf ("a twos %u\n", atwos));
156 result_bit1 ^= JACOBI_TWOS_U_BIT1 (atwos, blow);
157
158 /* establish shifted alow */
159 alow >>= atwos;
160 if (asize > 1)
161 {
162 asecond = asrcp[1];
163 if (atwos != 0)
164 alow |= (asecond << (GMP_NUMB_BITS - atwos)) & GMP_NUMB_MASK;
165 }
166
167 /* (a/2)=(2/a) with a odd */
168 result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
169
170 if (asize == 1 || (asize == 2 && (asecond >> atwos) == 0))
171 {
172 /* another special case with modexact and no copying */
173
174 if (alow == 1) /* (1/b)=1 always */
175 return JACOBI_BIT1_TO_PN (result_bit1);
176
177 /* b still has its twos, so cancel out their effect */
178 result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
179
180 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow); /* now (b/a) */
181 JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, blow, bsrcp, bsize, alow);
182 TRACE (printf ("base (%lu/%lu) with %d\n",
183 blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
184 return mpn_jacobi_base (blow, alow, result_bit1);
185 }
186
187
188 TMP_MARK;
189 TMP_ALLOC_LIMBS_2 (ap, asize, bp, bsize);
190
191 MPN_RSHIFT_OR_COPY (ap, asrcp, asize, atwos);
192 ASSERT (alow == ap[0]);
193 TRACE (mpn_trace ("stripped a", ap, asize));
194
195 MPN_RSHIFT_OR_COPY (bp, bsrcp, bsize, btwos);
196 ASSERT (blow == bp[0]);
197 TRACE (mpn_trace ("stripped b", bp, bsize));
198
199 /* swap if necessary to make a longer than b */
200 if (asize < bsize)
201 {
202 TRACE (printf ("swap\n"));
203 MPN_PTR_SWAP (ap,asize, bp,bsize);
204 MP_LIMB_T_SWAP (alow, blow);
205 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
206 }
207
208 /* If a is bigger than b then reduce to a mod b.
209 Division is much faster than chipping away at "a" bit-by-bit. */
210 if (asize > bsize)
211 {
212 mp_ptr rp, qp;
213
214 TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize));
215
216 TMP_ALLOC_LIMBS_2 (rp, bsize, qp, asize-bsize+1);
217 mpn_tdiv_qr (qp, rp, (mp_size_t) 0, ap, asize, bp, bsize);
218 ap = rp;
219 asize = bsize;
220 MPN_NORMALIZE (ap, asize);
221
222 TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize);
223 mpn_trace (" a", ap, asize);
224 mpn_trace (" b", bp, bsize));
225
226 if (asize == 0) /* (0/b)=0 for b!=1 */
227 goto zero;
228
229 alow = ap[0];
230 goto strip_a;
231 }
232
233 for (;;)
234 {
235 ASSERT (asize >= 1); /* a,b non-empty */
236 ASSERT (bsize >= 1);
237 ASSERT (ap[asize-1] != 0); /* a,b normalized (and hence non-zero) */
238 ASSERT (bp[bsize-1] != 0);
239 ASSERT (alow == ap[0]); /* low limb copies should be correct */
240 ASSERT (blow == bp[0]);
241 ASSERT (alow & 1); /* a,b odd */
242 ASSERT (blow & 1);
243
244 TRACE (printf ("top asize=%ld bsize=%ld\n", asize, bsize);
245 mpn_trace (" a", ap, asize);
246 mpn_trace (" b", bp, bsize));
247
248 /* swap if necessary to make a>=b, applying reciprocity
249 high limbs are almost always enough to tell which is bigger */
250 if (asize < bsize
251 || (asize == bsize
252 && ((ahigh=ap[asize-1]) < (bhigh=bp[asize-1])
253 || (ahigh == bhigh
254 && mpn_cmp (ap, bp, asize-1) < 0))))
255 {
256 TRACE (printf ("swap\n"));
257 MPN_PTR_SWAP (ap,asize, bp,bsize);
258 MP_LIMB_T_SWAP (alow, blow);
259 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
260 }
261
262 if (asize == 1)
263 break;
264
265 /* a = a-b */
266 ASSERT (asize >= bsize);
267 ASSERT_NOCARRY (mpn_sub (ap, ap, asize, bp, bsize));
268 MPN_NORMALIZE (ap, asize);
269 alow = ap[0];
270
271 /* (0/b)=0 for b!=1. b!=1 when a==0 because otherwise would have had
272 a==1 which is asize==1 and would have exited above. */
273 if (asize == 0)
274 goto zero;
275
276 strip_a:
277 /* low zero limbs on a can be discarded */
278 JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, ap, asize, alow);
279
280 if ((alow & 1) == 0)
281 {
282 /* factors of 2 from a */
283 unsigned twos;
284 count_trailing_zeros (twos, alow);
285 TRACE (printf ("twos %u\n", twos));
286 result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, blow);
287 ASSERT_NOCARRY (mpn_rshift (ap, ap, asize, twos));
288 asize -= (ap[asize-1] == 0);
289 alow = ap[0];
290 }
291 }
292
293 ASSERT (asize == 1 && bsize == 1); /* just alow and blow left */
294 TMP_FREE;
295
296 /* (1/b)=1 always (in this case have b==1 because a>=b) */
297 if (alow == 1)
298 return JACOBI_BIT1_TO_PN (result_bit1);
299
300 /* swap with reciprocity and do (b/a) */
301 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
302 TRACE (printf ("base (%lu/%lu) with %d\n",
303 blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
304 return mpn_jacobi_base (blow, alow, result_bit1);
305
306 zero:
307 TMP_FREE;
308 return 0;
309 }
310