1 /* mpn_gcdext -- Extended Greatest Common Divisor. 2 3 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008 Free Software 4 Foundation, Inc. 5 6 This file is part of the GNU MP Library. 7 8 The GNU MP Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 20 21 /* Default to binary gcdext_1, since it is best on most current machines. 22 We should teach tuneup to choose the right gcdext_1. */ 23 #define GCDEXT_1_USE_BINARY 1 24 25 #include "gmp.h" 26 #include "gmp-impl.h" 27 #include "longlong.h" 28 29 #ifndef NULL 30 # define NULL ((void *) 0) 31 #endif 32 33 /* FIXME: Takes two single-word limbs. It could be extended to a 34 * function that accepts a bignum for the first input, and only 35 * returns the first co-factor. */ 36 37 /* Returns g, u and v such that g = u A - v B. There are three 38 different cases for the result: 39 40 g = u A - v B, 0 < u < b, 0 < v < a 41 g = A u = 1, v = 0 42 g = B u = B, v = A - 1 43 44 We always return with 0 < u <= b, 0 <= v < a. 45 */ 46 #if GCDEXT_1_USE_BINARY 47 48 static mp_limb_t 49 gcdext_1_odd (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) 50 { 51 mp_limb_t u0; 52 mp_limb_t v0; 53 mp_limb_t v1; 54 mp_limb_t u1; 55 56 mp_limb_t B = b; 57 mp_limb_t A = a; 58 59 /* Through out this function maintain 60 61 a = u0 A - v0 B 62 b = u1 A - v1 B 63 64 where A and B are odd. */ 65 66 u0 = 1; v0 = 0; 67 u1 = b; v1 = a-1; 68 69 if (A == 1) 70 { 71 *up = u0; *vp = v0; 72 return 1; 73 } 74 else if (B == 1) 75 { 76 *up = u1; *vp = v1; 77 return 1; 78 } 79 80 while (a != b) 81 { 82 mp_limb_t mask; 83 84 ASSERT (a % 2 == 1); 85 ASSERT (b % 2 == 1); 86 87 ASSERT (0 < u0); ASSERT (u0 <= B); 88 ASSERT (0 < u1); ASSERT (u1 <= B); 89 90 ASSERT (0 <= v0); ASSERT (v0 < A); 91 ASSERT (0 <= v1); ASSERT (v1 < A); 92 93 if (a > b) 94 { 95 MP_LIMB_T_SWAP (a, b); 96 MP_LIMB_T_SWAP (u0, u1); 97 MP_LIMB_T_SWAP (v0, v1); 98 } 99 100 ASSERT (a < b); 101 102 /* Makes b even */ 103 b -= a; 104 105 mask = - (mp_limb_t) (u1 < u0); 106 u1 += B & mask; 107 v1 += A & mask; 108 u1 -= u0; 109 v1 -= v0; 110 111 ASSERT (b % 2 == 0); 112 113 do 114 { 115 /* As b = u1 A + v1 B is even, while A and B are odd, 116 either both or none of u1, v1 is even */ 117 118 ASSERT (u1 % 2 == v1 % 2); 119 120 mask = -(u1 & 1); 121 u1 = u1 / 2 + ((B / 2) & mask) - mask; 122 v1 = v1 / 2 + ((A / 2) & mask) - mask; 123 124 b /= 2; 125 } 126 while (b % 2 == 0); 127 } 128 129 /* Now g = a = b */ 130 ASSERT (a == b); 131 ASSERT (u1 <= B); 132 ASSERT (v1 < A); 133 134 ASSERT (A % a == 0); 135 ASSERT (B % a == 0); 136 ASSERT (u0 % (B/a) == u1 % (B/a)); 137 ASSERT (v0 % (A/a) == v1 % (A/a)); 138 139 *up = u0; *vp = v0; 140 141 return a; 142 } 143 144 mp_limb_t 145 mpn_gcdext_1 (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) 146 { 147 unsigned shift = 0; 148 mp_limb_t g; 149 mp_limb_t u; 150 mp_limb_t v; 151 152 /* We use unsigned values in the range 0, ... B - 1. As the values 153 are uniquely determined only modulo B, we can add B at will, to 154 get numbers in range or flip the least significant bit. */ 155 /* Deal with powers of two */ 156 while ((a | b) % 2 == 0) 157 { 158 a /= 2; b /= 2; shift++; 159 } 160 161 if (b % 2 == 0) 162 { 163 unsigned k = 0; 164 165 do { 166 b /= 2; k++; 167 } while (b % 2 == 0); 168 169 g = gcdext_1_odd (&u, &v, a, b); 170 171 while (k--) 172 { 173 /* We have g = u a + v b, and need to construct 174 g = u'a + v'(2b). 175 176 If v is even, we can just set u' = u, v' = v/2 177 If v is odd, we can set v' = (v + a)/2, u' = u + b 178 */ 179 180 if (v % 2 == 0) 181 v /= 2; 182 else 183 { 184 u = u + b; 185 v = v/2 + a/2 + 1; 186 } 187 b *= 2; 188 } 189 } 190 else if (a % 2 == 0) 191 { 192 unsigned k = 0; 193 194 do { 195 a /= 2; k++; 196 } while (a % 2 == 0); 197 198 g = gcdext_1_odd (&u, &v, a, b); 199 200 while (k--) 201 { 202 /* We have g = u a + v b, and need to construct 203 g = u'(2a) + v'b. 204 205 If u is even, we can just set u' = u/2, v' = v. 206 If u is odd, we can set u' = (u + b)/2 207 */ 208 209 if (u % 2 == 0) 210 u /= 2; 211 else 212 { 213 u = u/2 + b/2 + 1; 214 v = v + a; 215 } 216 a *= 2; 217 } 218 } 219 else 220 /* Ok, both are odd */ 221 g = gcdext_1_odd (&u, &v, a, b); 222 223 *up = u; 224 *vp = v; 225 226 return g << shift; 227 } 228 229 #else /* ! GCDEXT_1_USE_BINARY */ 230 static mp_limb_t 231 gcdext_1_u (mp_limb_t *up, mp_limb_t a, mp_limb_t b) 232 { 233 /* Maintain 234 235 a = u0 A mod B 236 b = - u1 A mod B 237 */ 238 mp_limb_t u0 = 1; 239 mp_limb_t u1 = 0; 240 mp_limb_t B = b; 241 242 ASSERT (a >= b); 243 ASSERT (b > 0); 244 245 for (;;) 246 { 247 mp_limb_t q; 248 249 q = a / b; 250 a -= q * b; 251 252 if (a == 0) 253 { 254 *up = B - u1; 255 return b; 256 } 257 u0 += q * u1; 258 259 q = b / a; 260 b -= q * a; 261 262 if (b == 0) 263 { 264 *up = u0; 265 return a; 266 } 267 u1 += q * u0; 268 } 269 } 270 271 mp_limb_t 272 mpn_gcdext_1 (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) 273 { 274 /* Maintain 275 276 a = u0 A - v0 B 277 b = - u1 A + v1 B = (B - u1) A - (A - v1) B 278 */ 279 mp_limb_t u0 = 1; 280 mp_limb_t v0 = 0; 281 mp_limb_t u1 = 0; 282 mp_limb_t v1 = 1; 283 284 mp_limb_t A = a; 285 mp_limb_t B = b; 286 287 ASSERT (a >= b); 288 ASSERT (b > 0); 289 290 for (;;) 291 { 292 mp_limb_t q; 293 294 q = a / b; 295 a -= q * b; 296 297 if (a == 0) 298 { 299 *up = B - u1; 300 *vp = A - v1; 301 return b; 302 } 303 u0 += q * u1; 304 v0 += q * v1; 305 306 q = b / a; 307 b -= q * a; 308 309 if (b == 0) 310 { 311 *up = u0; 312 *vp = v0; 313 return a; 314 } 315 u1 += q * u0; 316 v1 += q * v0; 317 } 318 } 319 #endif /* ! GCDEXT_1_USE_BINARY */ 320