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