1 /* jacobi_2.c 2 3 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY 4 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 5 GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. 6 7 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2008, 2010 Free Software 8 Foundation, Inc. 9 10 This file is part of the GNU MP Library. 11 12 The GNU MP Library is free software; you can redistribute it and/or modify 13 it under the terms of the GNU Lesser General Public License as published by 14 the Free Software Foundation; either version 3 of the License, or (at your 15 option) any later version. 16 17 The GNU MP Library is distributed in the hope that it will be useful, but 18 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 19 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 20 License for more details. 21 22 You should have received a copy of the GNU Lesser General Public License 23 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 24 25 #include "gmp.h" 26 #include "gmp-impl.h" 27 #include "longlong.h" 28 29 #ifndef JACOBI_2_METHOD 30 #define JACOBI_2_METHOD 2 31 #endif 32 33 /* Computes (a / b) where b is odd, and a and b are otherwise arbitrary 34 two-limb numbers. */ 35 #if JACOBI_2_METHOD == 1 36 int 37 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit) 38 { 39 mp_limb_t ah, al, bh, bl; 40 int c; 41 42 al = ap[0]; 43 ah = ap[1]; 44 bl = bp[0]; 45 bh = bp[1]; 46 47 ASSERT (bl & 1); 48 49 bl = ((bh << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK) | (bl >> 1); 50 bh >>= 1; 51 52 if ( (bh | bl) == 0) 53 return 1 - 2*(bit & 1); 54 55 if ( (ah | al) == 0) 56 return 0; 57 58 if (al == 0) 59 { 60 al = ah; 61 ah = 0; 62 bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1)); 63 } 64 count_trailing_zeros (c, al); 65 bit ^= c & (bl ^ (bl >> 1)); 66 67 c++; 68 if (UNLIKELY (c == GMP_NUMB_BITS)) 69 { 70 al = ah; 71 ah = 0; 72 } 73 else 74 { 75 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c); 76 ah >>= c; 77 } 78 while ( (ah | bh) > 0) 79 { 80 mp_limb_t th, tl; 81 mp_limb_t bgta; 82 83 sub_ddmmss (th, tl, ah, al, bh, bl); 84 if ( (tl | th) == 0) 85 return 0; 86 87 bgta = LIMB_HIGHBIT_TO_MASK (th); 88 89 /* If b > a, invoke reciprocity */ 90 bit ^= (bgta & al & bl); 91 92 /* b <-- min (a, b) */ 93 add_ssaaaa (bh, bl, bh, bl, th & bgta, tl & bgta); 94 95 if ( (bh | bl) == 0) 96 return 1 - 2*(bit & 1); 97 98 /* a <-- |a - b| */ 99 al = (bgta ^ tl) - bgta; 100 ah = (bgta ^ th); 101 102 if (UNLIKELY (al == 0)) 103 { 104 /* If b > a, al == 0 implies that we have a carry to 105 propagate. */ 106 al = ah - bgta; 107 ah = 0; 108 bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1)); 109 } 110 count_trailing_zeros (c, al); 111 c++; 112 bit ^= c & (bl ^ (bl >> 1)); 113 114 if (UNLIKELY (c == GMP_NUMB_BITS)) 115 { 116 al = ah; 117 ah = 0; 118 } 119 else 120 { 121 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c); 122 ah >>= c; 123 } 124 } 125 126 ASSERT (bl > 0); 127 128 while ( (al | bl) & GMP_LIMB_HIGHBIT) 129 { 130 /* Need an extra comparison to get the mask. */ 131 mp_limb_t t = al - bl; 132 mp_limb_t bgta = - (bl > al); 133 134 if (t == 0) 135 return 0; 136 137 /* If b > a, invoke reciprocity */ 138 bit ^= (bgta & al & bl); 139 140 /* b <-- min (a, b) */ 141 bl += (bgta & t); 142 143 /* a <-- |a - b| */ 144 al = (t ^ bgta) - bgta; 145 146 /* Number of trailing zeros is the same no matter if we look at 147 * t or a, but using t gives more parallelism. */ 148 count_trailing_zeros (c, t); 149 c ++; 150 /* (2/b) = -1 if b = 3 or 5 mod 8 */ 151 bit ^= c & (bl ^ (bl >> 1)); 152 153 if (UNLIKELY (c == GMP_NUMB_BITS)) 154 return 1 - 2*(bit & 1); 155 156 al >>= c; 157 } 158 159 /* Here we have a little impedance mismatch. Better to inline it? */ 160 return mpn_jacobi_base (2*al+1, 2*bl+1, bit << 1); 161 } 162 #elif JACOBI_2_METHOD == 2 163 int 164 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit) 165 { 166 mp_limb_t ah, al, bh, bl; 167 int c; 168 169 al = ap[0]; 170 ah = ap[1]; 171 bl = bp[0]; 172 bh = bp[1]; 173 174 ASSERT (bl & 1); 175 176 /* Use bit 1. */ 177 bit <<= 1; 178 179 if (bh == 0 && bl == 1) 180 /* (a|1) = 1 */ 181 return 1 - (bit & 2); 182 183 if (al == 0) 184 { 185 if (ah == 0) 186 /* (0|b) = 0, b > 1 */ 187 return 0; 188 189 count_trailing_zeros (c, ah); 190 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1)); 191 192 al = bl; 193 bl = ah >> c; 194 195 if (bl == 1) 196 /* (1|b) = 1 */ 197 return 1 - (bit & 2); 198 199 ah = bh; 200 201 bit ^= al & bl; 202 203 goto b_reduced; 204 } 205 if ( (al & 1) == 0) 206 { 207 count_trailing_zeros (c, al); 208 209 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c); 210 ah >>= c; 211 bit ^= (c << 1) & (bl ^ (bl >> 1)); 212 } 213 if (ah == 0) 214 { 215 if (bh > 0) 216 { 217 bit ^= al & bl; 218 MP_LIMB_T_SWAP (al, bl); 219 ah = bh; 220 goto b_reduced; 221 } 222 goto ab_reduced; 223 } 224 225 while (bh > 0) 226 { 227 /* Compute (a|b) */ 228 while (ah > bh) 229 { 230 sub_ddmmss (ah, al, ah, al, bh, bl); 231 if (al == 0) 232 { 233 count_trailing_zeros (c, ah); 234 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1)); 235 236 al = bl; 237 bl = ah >> c; 238 ah = bh; 239 240 bit ^= al & bl; 241 goto b_reduced; 242 } 243 count_trailing_zeros (c, al); 244 bit ^= (c << 1) & (bl ^ (bl >> 1)); 245 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c); 246 ah >>= c; 247 } 248 if (ah == bh) 249 goto cancel_hi; 250 251 if (ah == 0) 252 { 253 bit ^= al & bl; 254 MP_LIMB_T_SWAP (al, bl); 255 ah = bh; 256 break; 257 } 258 259 bit ^= al & bl; 260 261 /* Compute (b|a) */ 262 while (bh > ah) 263 { 264 sub_ddmmss (bh, bl, bh, bl, ah, al); 265 if (bl == 0) 266 { 267 count_trailing_zeros (c, bh); 268 bit ^= ((GMP_NUMB_BITS + c) << 1) & (al ^ (al >> 1)); 269 270 bl = bh >> c; 271 bit ^= al & bl; 272 goto b_reduced; 273 } 274 count_trailing_zeros (c, bl); 275 bit ^= (c << 1) & (al ^ (al >> 1)); 276 bl = ((bh << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (bl >> c); 277 bh >>= c; 278 } 279 bit ^= al & bl; 280 281 /* Compute (a|b) */ 282 if (ah == bh) 283 { 284 cancel_hi: 285 if (al < bl) 286 { 287 MP_LIMB_T_SWAP (al, bl); 288 bit ^= al & bl; 289 } 290 al -= bl; 291 if (al == 0) 292 return 0; 293 294 count_trailing_zeros (c, al); 295 bit ^= (c << 1) & (bl ^ (bl >> 1)); 296 al >>= c; 297 298 if (al == 1) 299 return 1 - (bit & 2); 300 301 MP_LIMB_T_SWAP (al, bl); 302 bit ^= al & bl; 303 break; 304 } 305 } 306 307 b_reduced: 308 /* Compute (a|b), with b a single limb. */ 309 ASSERT (bl & 1); 310 311 if (bl == 1) 312 /* (a|1) = 1 */ 313 return 1 - (bit & 2); 314 315 while (ah > 0) 316 { 317 ah -= (al < bl); 318 al -= bl; 319 if (al == 0) 320 { 321 if (ah == 0) 322 return 0; 323 count_trailing_zeros (c, ah); 324 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1)); 325 al = ah >> c; 326 goto ab_reduced; 327 } 328 count_trailing_zeros (c, al); 329 330 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c); 331 ah >>= c; 332 bit ^= (c << 1) & (bl ^ (bl >> 1)); 333 } 334 ab_reduced: 335 ASSERT (bl & 1); 336 ASSERT (bl > 1); 337 338 return mpn_jacobi_base (al, bl, bit); 339 } 340 #else 341 #error Unsupported value for JACOBI_2_METHOD 342 #endif 343