1 /* hgcd2_jacobi.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, 2011 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 #if GMP_NAIL_BITS > 0 39 #error Nails not supported. 40 #endif 41 42 /* FIXME: Duplicated in hgcd2.c. Should move to gmp-impl.h, and 43 possibly be renamed. */ 44 static inline mp_limb_t 45 div1 (mp_ptr rp, 46 mp_limb_t n0, 47 mp_limb_t d0) 48 { 49 mp_limb_t q = 0; 50 51 if ((mp_limb_signed_t) n0 < 0) 52 { 53 int cnt; 54 for (cnt = 1; (mp_limb_signed_t) d0 >= 0; cnt++) 55 { 56 d0 = d0 << 1; 57 } 58 59 q = 0; 60 while (cnt) 61 { 62 q <<= 1; 63 if (n0 >= d0) 64 { 65 n0 = n0 - d0; 66 q |= 1; 67 } 68 d0 = d0 >> 1; 69 cnt--; 70 } 71 } 72 else 73 { 74 int cnt; 75 for (cnt = 0; n0 >= d0; cnt++) 76 { 77 d0 = d0 << 1; 78 } 79 80 q = 0; 81 while (cnt) 82 { 83 d0 = d0 >> 1; 84 q <<= 1; 85 if (n0 >= d0) 86 { 87 n0 = n0 - d0; 88 q |= 1; 89 } 90 cnt--; 91 } 92 } 93 *rp = n0; 94 return q; 95 } 96 97 /* Two-limb division optimized for small quotients. */ 98 static inline mp_limb_t 99 div2 (mp_ptr rp, 100 mp_limb_t nh, mp_limb_t nl, 101 mp_limb_t dh, mp_limb_t dl) 102 { 103 mp_limb_t q = 0; 104 105 if ((mp_limb_signed_t) nh < 0) 106 { 107 int cnt; 108 for (cnt = 1; (mp_limb_signed_t) dh >= 0; cnt++) 109 { 110 dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1)); 111 dl = dl << 1; 112 } 113 114 while (cnt) 115 { 116 q <<= 1; 117 if (nh > dh || (nh == dh && nl >= dl)) 118 { 119 sub_ddmmss (nh, nl, nh, nl, dh, dl); 120 q |= 1; 121 } 122 dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1); 123 dh = dh >> 1; 124 cnt--; 125 } 126 } 127 else 128 { 129 int cnt; 130 for (cnt = 0; nh > dh || (nh == dh && nl >= dl); cnt++) 131 { 132 dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1)); 133 dl = dl << 1; 134 } 135 136 while (cnt) 137 { 138 dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1); 139 dh = dh >> 1; 140 q <<= 1; 141 if (nh > dh || (nh == dh && nl >= dl)) 142 { 143 sub_ddmmss (nh, nl, nh, nl, dh, dl); 144 q |= 1; 145 } 146 cnt--; 147 } 148 } 149 150 rp[0] = nl; 151 rp[1] = nh; 152 153 return q; 154 } 155 156 int 157 mpn_hgcd2_jacobi (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl, 158 struct hgcd_matrix1 *M, unsigned *bitsp) 159 { 160 mp_limb_t u00, u01, u10, u11; 161 unsigned bits = *bitsp; 162 163 if (ah < 2 || bh < 2) 164 return 0; 165 166 if (ah > bh || (ah == bh && al > bl)) 167 { 168 sub_ddmmss (ah, al, ah, al, bh, bl); 169 if (ah < 2) 170 return 0; 171 172 u00 = u01 = u11 = 1; 173 u10 = 0; 174 bits = mpn_jacobi_update (bits, 1, 1); 175 } 176 else 177 { 178 sub_ddmmss (bh, bl, bh, bl, ah, al); 179 if (bh < 2) 180 return 0; 181 182 u00 = u10 = u11 = 1; 183 u01 = 0; 184 bits = mpn_jacobi_update (bits, 0, 1); 185 } 186 187 if (ah < bh) 188 goto subtract_a; 189 190 for (;;) 191 { 192 ASSERT (ah >= bh); 193 if (ah == bh) 194 goto done; 195 196 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2))) 197 { 198 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2)); 199 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2)); 200 201 break; 202 } 203 204 /* Subtract a -= q b, and multiply M from the right by (1 q ; 0 205 1), affecting the second column of M. */ 206 ASSERT (ah > bh); 207 sub_ddmmss (ah, al, ah, al, bh, bl); 208 209 if (ah < 2) 210 goto done; 211 212 if (ah <= bh) 213 { 214 /* Use q = 1 */ 215 u01 += u00; 216 u11 += u10; 217 bits = mpn_jacobi_update (bits, 1, 1); 218 } 219 else 220 { 221 mp_limb_t r[2]; 222 mp_limb_t q = div2 (r, ah, al, bh, bl); 223 al = r[0]; ah = r[1]; 224 if (ah < 2) 225 { 226 /* A is too small, but q is correct. */ 227 u01 += q * u00; 228 u11 += q * u10; 229 bits = mpn_jacobi_update (bits, 1, q & 3); 230 goto done; 231 } 232 q++; 233 u01 += q * u00; 234 u11 += q * u10; 235 bits = mpn_jacobi_update (bits, 1, q & 3); 236 } 237 subtract_a: 238 ASSERT (bh >= ah); 239 if (ah == bh) 240 goto done; 241 242 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2))) 243 { 244 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2)); 245 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2)); 246 247 goto subtract_a1; 248 } 249 250 /* Subtract b -= q a, and multiply M from the right by (1 0 ; q 251 1), affecting the first column of M. */ 252 sub_ddmmss (bh, bl, bh, bl, ah, al); 253 254 if (bh < 2) 255 goto done; 256 257 if (bh <= ah) 258 { 259 /* Use q = 1 */ 260 u00 += u01; 261 u10 += u11; 262 bits = mpn_jacobi_update (bits, 0, 1); 263 } 264 else 265 { 266 mp_limb_t r[2]; 267 mp_limb_t q = div2 (r, bh, bl, ah, al); 268 bl = r[0]; bh = r[1]; 269 if (bh < 2) 270 { 271 /* B is too small, but q is correct. */ 272 u00 += q * u01; 273 u10 += q * u11; 274 bits = mpn_jacobi_update (bits, 0, q & 3); 275 goto done; 276 } 277 q++; 278 u00 += q * u01; 279 u10 += q * u11; 280 bits = mpn_jacobi_update (bits, 0, q & 3); 281 } 282 } 283 284 /* NOTE: Since we discard the least significant half limb, we don't 285 get a truly maximal M (corresponding to |a - b| < 286 2^{GMP_LIMB_BITS +1}). */ 287 /* Single precision loop */ 288 for (;;) 289 { 290 ASSERT (ah >= bh); 291 if (ah == bh) 292 break; 293 294 ah -= bh; 295 if (ah < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1))) 296 break; 297 298 if (ah <= bh) 299 { 300 /* Use q = 1 */ 301 u01 += u00; 302 u11 += u10; 303 bits = mpn_jacobi_update (bits, 1, 1); 304 } 305 else 306 { 307 mp_limb_t r; 308 mp_limb_t q = div1 (&r, ah, bh); 309 ah = r; 310 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1))) 311 { 312 /* A is too small, but q is correct. */ 313 u01 += q * u00; 314 u11 += q * u10; 315 bits = mpn_jacobi_update (bits, 1, q & 3); 316 break; 317 } 318 q++; 319 u01 += q * u00; 320 u11 += q * u10; 321 bits = mpn_jacobi_update (bits, 1, q & 3); 322 } 323 subtract_a1: 324 ASSERT (bh >= ah); 325 if (ah == bh) 326 break; 327 328 bh -= ah; 329 if (bh < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1))) 330 break; 331 332 if (bh <= ah) 333 { 334 /* Use q = 1 */ 335 u00 += u01; 336 u10 += u11; 337 bits = mpn_jacobi_update (bits, 0, 1); 338 } 339 else 340 { 341 mp_limb_t r; 342 mp_limb_t q = div1 (&r, bh, ah); 343 bh = r; 344 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1))) 345 { 346 /* B is too small, but q is correct. */ 347 u00 += q * u01; 348 u10 += q * u11; 349 bits = mpn_jacobi_update (bits, 0, q & 3); 350 break; 351 } 352 q++; 353 u00 += q * u01; 354 u10 += q * u11; 355 bits = mpn_jacobi_update (bits, 0, q & 3); 356 } 357 } 358 359 done: 360 M->u[0][0] = u00; M->u[0][1] = u01; 361 M->u[1][0] = u10; M->u[1][1] = u11; 362 *bitsp = bits; 363 364 return 1; 365 } 366