Line | Count | Source (jump to first uncovered line) |
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 |
58 | | mpz_jacobi (mpz_srcptr a, mpz_srcptr b) |
59 | 720 | { |
60 | 720 | mp_srcptr asrcp, bsrcp; |
61 | 720 | mp_size_t asize, bsize; |
62 | 720 | mp_limb_t alow, blow; |
63 | 720 | mp_ptr ap, bp; |
64 | 720 | unsigned btwos; |
65 | 720 | int result_bit1; |
66 | 720 | int res; |
67 | 720 | TMP_DECL; |
68 | | |
69 | 720 | asize = SIZ(a); |
70 | 720 | asrcp = PTR(a); |
71 | 720 | alow = asrcp[0]; |
72 | | |
73 | 720 | bsize = SIZ(b); |
74 | 720 | bsrcp = PTR(b); |
75 | 720 | 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 | 720 | if (bsize == 0) |
83 | | /* (a/0) = [ a = 1 or a = -1 ] */ |
84 | 9 | return JACOBI_LS0 (alow, asize); |
85 | | |
86 | 711 | if (asize == 0) |
87 | | /* (0/b) = [ b = 1 or b = - 1 ] */ |
88 | 5 | return JACOBI_0LS (blow, bsize); |
89 | | |
90 | 706 | if ( (((alow | blow) & 1) == 0)) |
91 | | /* Common factor of 2 ==> (a/b) = 0 */ |
92 | 8 | return 0; |
93 | | |
94 | 698 | if (bsize < 0) |
95 | 0 | { |
96 | | /* (a/-1) = -1 if a < 0, +1 if a >= 0 */ |
97 | 0 | result_bit1 = (asize < 0) << 1; |
98 | 0 | bsize = -bsize; |
99 | 0 | } |
100 | 698 | else |
101 | 698 | result_bit1 = 0; |
102 | | |
103 | 698 | JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow); |
104 | | |
105 | 698 | count_trailing_zeros (btwos, blow); |
106 | 698 | blow >>= btwos; |
107 | | |
108 | 698 | if (bsize > 1 && btwos > 0) |
109 | 111 | { |
110 | 111 | mp_limb_t b1 = bsrcp[1]; |
111 | 111 | blow |= b1 << (GMP_NUMB_BITS - btwos); |
112 | 111 | if (bsize == 2 && (b1 >> btwos) == 0) |
113 | 0 | bsize = 1; |
114 | 111 | } |
115 | | |
116 | 698 | if (asize < 0) |
117 | 0 | { |
118 | | /* (-1/b) = -1 iff b = 3 (mod 4) */ |
119 | 0 | result_bit1 ^= JACOBI_N1B_BIT1(blow); |
120 | 0 | asize = -asize; |
121 | 0 | } |
122 | | |
123 | 698 | 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 | 698 | if (asize < bsize) |
129 | 237 | { |
130 | 237 | MPN_SRCPTR_SWAP (asrcp, asize, bsrcp, bsize); |
131 | 237 | 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 | 237 | count_trailing_zeros (btwos, blow); |
146 | 237 | blow >>= btwos; |
147 | | |
148 | 237 | if (bsize > 1 && btwos > 0) |
149 | 124 | { |
150 | 124 | mp_limb_t b1 = bsrcp[1]; |
151 | 124 | blow |= b1 << (GMP_NUMB_BITS - btwos); |
152 | 124 | if (bsize == 2 && (b1 >> btwos) == 0) |
153 | 2 | bsize = 1; |
154 | 124 | } |
155 | | |
156 | 237 | result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow); |
157 | 237 | } |
158 | | |
159 | 698 | if (bsize == 1) |
160 | 56 | { |
161 | 56 | result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow); |
162 | | |
163 | 56 | if (blow == 1) |
164 | 4 | return JACOBI_BIT1_TO_PN (result_bit1); |
165 | | |
166 | 52 | if (asize > 1) |
167 | 15 | JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow); |
168 | | |
169 | 52 | return mpn_jacobi_base (alow, blow, result_bit1); |
170 | 52 | } |
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 | 642 | TMP_MARK; |
178 | | |
179 | 642 | if (asize >= 2*bsize) |
180 | 66 | TMP_ALLOC_LIMBS_2 (ap, bsize, bp, asize - bsize + 1); |
181 | 576 | else |
182 | 576 | 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 | 642 | if (asize > bsize) |
190 | 455 | mpn_tdiv_qr (bp, ap, 0, asrcp, asize, bsrcp, bsize); |
191 | 187 | else |
192 | 187 | MPN_COPY (ap, asrcp, bsize); |
193 | | |
194 | 642 | if (btwos > 0) |
195 | 213 | { |
196 | 213 | result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow); |
197 | | |
198 | 213 | ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, btwos)); |
199 | 213 | bsize -= (ap[bsize-1] | bp[bsize-1]) == 0; |
200 | 213 | } |
201 | 429 | else |
202 | 429 | MPN_COPY (bp, bsrcp, bsize); |
203 | | |
204 | 642 | ASSERT (blow == bp[0]); |
205 | 642 | res = mpn_jacobi_n (ap, bp, bsize, |
206 | 642 | mpn_jacobi_init (ap[0], blow, (result_bit1>>1) & 1)); |
207 | | |
208 | 642 | TMP_FREE; |
209 | 642 | return res; |
210 | 642 | } |