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