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