1 /* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
2
3 Copyright 2000-2002, 2005, 2010-2012 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
8 it under the terms of either:
9
10 * the GNU Lesser General Public License as published by the Free
11 Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 or
15
16 * the GNU General Public License as published by the Free Software
17 Foundation; either version 2 of the License, or (at your option) any
18 later version.
19
20 or both in parallel, as here.
21
22 The GNU MP Library is distributed in the hope that it will be useful, but
23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25 for more details.
26
27 You should have received copies of the GNU General Public License and the
28 GNU Lesser General Public License along with the GNU MP Library. If not,
29 see https://www.gnu.org/licenses/. */
30
31 #include <stdio.h>
32 #include "gmp-impl.h"
33 #include "longlong.h"
34
35
36 /* This code does triple duty as mpz_jacobi, mpz_legendre and
37 mpz_kronecker. For ABI compatibility, the link symbol is
38 __gmpz_jacobi, not __gmpz_kronecker, even though the latter would
39 be more logical.
40
41 mpz_jacobi could assume b is odd, but the improvements from that seem
42 small compared to other operations, and anything significant should be
43 checked at run-time since we'd like odd b to go fast in mpz_kronecker
44 too.
45
46 mpz_legendre could assume b is an odd prime, but knowing this doesn't
47 present any obvious benefits. Result 0 wouldn't arise (unless "a" is a
48 multiple of b), but the checking for that takes little time compared to
49 other operations.
50
51 Enhancements:
52
53 mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
54
55 */
56
57 int
mpz_jacobi(mpz_srcptr a,mpz_srcptr b)58 mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
59 {
60 mp_srcptr asrcp, bsrcp;
61 mp_size_t asize, bsize;
62 mp_limb_t alow, blow;
63 mp_ptr ap, bp;
64 unsigned btwos;
65 int result_bit1;
66 int res;
67 TMP_DECL;
68
69 asize = SIZ(a);
70 asrcp = PTR(a);
71 alow = asrcp[0];
72
73 bsize = SIZ(b);
74 bsrcp = PTR(b);
75 blow = bsrcp[0];
76
77 /* The MPN jacobi functions require positive a and b, and b odd. So
78 we must to handle the cases of a or b zero, then signs, and then
79 the case of even b.
80 */
81
82 if (bsize == 0)
83 /* (a/0) = [ a = 1 or a = -1 ] */
84 return JACOBI_LS0 (alow, asize);
85
86 if (asize == 0)
87 /* (0/b) = [ b = 1 or b = - 1 ] */
88 return JACOBI_0LS (blow, bsize);
89
90 if ( (((alow | blow) & 1) == 0))
91 /* Common factor of 2 ==> (a/b) = 0 */
92 return 0;
93
94 if (bsize < 0)
95 {
96 /* (a/-1) = -1 if a < 0, +1 if a >= 0 */
97 result_bit1 = (asize < 0) << 1;
98 bsize = -bsize;
99 }
100 else
101 result_bit1 = 0;
102
103 JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
104
105 count_trailing_zeros (btwos, blow);
106 blow >>= btwos;
107
108 if (bsize > 1 && btwos > 0)
109 {
110 mp_limb_t b1 = bsrcp[1];
111 blow |= b1 << (GMP_NUMB_BITS - btwos);
112 if (bsize == 2 && (b1 >> btwos) == 0)
113 bsize = 1;
114 }
115
116 if (asize < 0)
117 {
118 /* (-1/b) = -1 iff b = 3 (mod 4) */
119 result_bit1 ^= JACOBI_N1B_BIT1(blow);
120 asize = -asize;
121 }
122
123 JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
124
125 /* Ensure asize >= bsize. Take advantage of the generalized
126 reciprocity law (a/b*2^n) = (b*2^n / a) * RECIP(a,b) */
127
128 if (asize < bsize)
129 {
130 MPN_SRCPTR_SWAP (asrcp, asize, bsrcp, bsize);
131 MP_LIMB_T_SWAP (alow, blow);
132
133 /* NOTE: The value of alow (old blow) is a bit subtle. For this code
134 path, we get alow as the low, always odd, limb of shifted A. Which is
135 what we need for the reciprocity update below.
136
137 However, all other uses of alow assumes that it is *not*
138 shifted. Luckily, alow matters only when either
139
140 + btwos > 0, in which case A is always odd
141
142 + asize == bsize == 1, in which case this code path is never
143 taken. */
144
145 count_trailing_zeros (btwos, blow);
146 blow >>= btwos;
147
148 if (bsize > 1 && btwos > 0)
149 {
150 mp_limb_t b1 = bsrcp[1];
151 blow |= b1 << (GMP_NUMB_BITS - btwos);
152 if (bsize == 2 && (b1 >> btwos) == 0)
153 bsize = 1;
154 }
155
156 result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
157 }
158
159 if (bsize == 1)
160 {
161 result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
162
163 if (blow == 1)
164 return JACOBI_BIT1_TO_PN (result_bit1);
165
166 if (asize > 1)
167 JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
168
169 return mpn_jacobi_base (alow, blow, result_bit1);
170 }
171
172 /* Allocation strategy: For A, we allocate a working copy only for A % B, but
173 when A is much larger than B, we have to allocate space for the large
174 quotient. We use the same area, pointed to by bp, for both the quotient
175 A/B and the working copy of B. */
176
177 TMP_MARK;
178
179 if (asize >= 2*bsize)
180 TMP_ALLOC_LIMBS_2 (ap, bsize, bp, asize - bsize + 1);
181 else
182 TMP_ALLOC_LIMBS_2 (ap, bsize, bp, bsize);
183
184 /* In the case of even B, we conceptually shift out the powers of two first,
185 and then divide A mod B. Hence, when taking those powers of two into
186 account, we must use alow *before* the division. Doing the actual division
187 first is ok, because the point is to remove multiples of B from A, and
188 multiples of 2^k B are good enough. */
189 if (asize > bsize)
190 mpn_tdiv_qr (bp, ap, 0, asrcp, asize, bsrcp, bsize);
191 else
192 MPN_COPY (ap, asrcp, bsize);
193
194 if (btwos > 0)
195 {
196 result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
197
198 ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, btwos));
199 bsize -= (ap[bsize-1] | bp[bsize-1]) == 0;
200 }
201 else
202 MPN_COPY (bp, bsrcp, bsize);
203
204 ASSERT (blow == bp[0]);
205 res = mpn_jacobi_n (ap, bp, bsize,
206 mpn_jacobi_init (ap[0], blow, (result_bit1>>1) & 1));
207
208 TMP_FREE;
209 return res;
210 }
211