1 /* mpn_gcdext -- Extended Greatest Common Divisor. 2 3 Copyright 1996, 1998, 2000-2005, 2008, 2009, 2012 Free Software Foundation, 4 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 either: 10 11 * the GNU Lesser General Public License as published by the Free 12 Software Foundation; either version 3 of the License, or (at your 13 option) any later version. 14 15 or 16 17 * the GNU General Public License as published by the Free Software 18 Foundation; either version 2 of the License, or (at your option) any 19 later version. 20 21 or both in parallel, as here. 22 23 The GNU MP Library is distributed in the hope that it will be useful, but 24 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 25 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 26 for more details. 27 28 You should have received copies of the GNU General Public License and the 29 GNU Lesser General Public License along with the GNU MP Library. If not, 30 see https://www.gnu.org/licenses/. */ 31 32 #include "gmp-impl.h" 33 #include "longlong.h" 34 35 /* Computes (r;b) = (a; b) M. Result is of size n + M->n +/- 1, and 36 the size is returned (if inputs are non-normalized, result may be 37 non-normalized too). Temporary space needed is M->n + n. 38 */ 39 static size_t 40 hgcd_mul_matrix_vector (struct hgcd_matrix *M, 41 mp_ptr rp, mp_srcptr ap, mp_ptr bp, mp_size_t n, mp_ptr tp) 42 { 43 mp_limb_t ah, bh; 44 45 /* Compute (r,b) <-- (u00 a + u10 b, u01 a + u11 b) as 46 47 t = u00 * a 48 r = u10 * b 49 r += t; 50 51 t = u11 * b 52 b = u01 * a 53 b += t; 54 */ 55 56 if (M->n >= n) 57 { 58 mpn_mul (tp, M->p[0][0], M->n, ap, n); 59 mpn_mul (rp, M->p[1][0], M->n, bp, n); 60 } 61 else 62 { 63 mpn_mul (tp, ap, n, M->p[0][0], M->n); 64 mpn_mul (rp, bp, n, M->p[1][0], M->n); 65 } 66 67 ah = mpn_add_n (rp, rp, tp, n + M->n); 68 69 if (M->n >= n) 70 { 71 mpn_mul (tp, M->p[1][1], M->n, bp, n); 72 mpn_mul (bp, M->p[0][1], M->n, ap, n); 73 } 74 else 75 { 76 mpn_mul (tp, bp, n, M->p[1][1], M->n); 77 mpn_mul (bp, ap, n, M->p[0][1], M->n); 78 } 79 bh = mpn_add_n (bp, bp, tp, n + M->n); 80 81 n += M->n; 82 if ( (ah | bh) > 0) 83 { 84 rp[n] = ah; 85 bp[n] = bh; 86 n++; 87 } 88 else 89 { 90 /* Normalize */ 91 while ( (rp[n-1] | bp[n-1]) == 0) 92 n--; 93 } 94 95 return n; 96 } 97 98 #define COMPUTE_V_ITCH(n) (2*(n)) 99 100 /* Computes |v| = |(g - u a)| / b, where u may be positive or 101 negative, and v is of the opposite sign. max(a, b) is of size n, u and 102 v at most size n, and v must have space for n+1 limbs. */ 103 static mp_size_t 104 compute_v (mp_ptr vp, 105 mp_srcptr ap, mp_srcptr bp, mp_size_t n, 106 mp_srcptr gp, mp_size_t gn, 107 mp_srcptr up, mp_size_t usize, 108 mp_ptr tp) 109 { 110 mp_size_t size; 111 mp_size_t an; 112 mp_size_t bn; 113 mp_size_t vn; 114 115 ASSERT (n > 0); 116 ASSERT (gn > 0); 117 ASSERT (usize != 0); 118 119 size = ABS (usize); 120 ASSERT (size <= n); 121 ASSERT (up[size-1] > 0); 122 123 an = n; 124 MPN_NORMALIZE (ap, an); 125 ASSERT (gn <= an); 126 127 if (an >= size) 128 mpn_mul (tp, ap, an, up, size); 129 else 130 mpn_mul (tp, up, size, ap, an); 131 132 size += an; 133 134 if (usize > 0) 135 { 136 /* |v| = -v = (u a - g) / b */ 137 138 ASSERT_NOCARRY (mpn_sub (tp, tp, size, gp, gn)); 139 MPN_NORMALIZE (tp, size); 140 if (size == 0) 141 return 0; 142 } 143 else 144 { /* |v| = v = (g - u a) / b = (g + |u| a) / b. Since g <= a, 145 (g + |u| a) always fits in (|usize| + an) limbs. */ 146 147 ASSERT_NOCARRY (mpn_add (tp, tp, size, gp, gn)); 148 size -= (tp[size - 1] == 0); 149 } 150 151 /* Now divide t / b. There must be no remainder */ 152 bn = n; 153 MPN_NORMALIZE (bp, bn); 154 ASSERT (size >= bn); 155 156 vn = size + 1 - bn; 157 ASSERT (vn <= n + 1); 158 159 mpn_divexact (vp, tp, size, bp, bn); 160 vn -= (vp[vn-1] == 0); 161 162 return vn; 163 } 164 165 /* Temporary storage: 166 167 Initial division: Quotient of at most an - n + 1 <= an limbs. 168 169 Storage for u0 and u1: 2(n+1). 170 171 Storage for hgcd matrix M, with input ceil(n/2): 5 * ceil(n/4) 172 173 Storage for hgcd, input (n + 1)/2: 9 n/4 plus some. 174 175 When hgcd succeeds: 1 + floor(3n/2) for adjusting a and b, and 2(n+1) for the cofactors. 176 177 When hgcd fails: 2n + 1 for mpn_gcdext_subdiv_step, which is less. 178 179 For the lehmer call after the loop, Let T denote 180 GCDEXT_DC_THRESHOLD. For the gcdext_lehmer call, we need T each for 181 u, a and b, and 4T+3 scratch space. Next, for compute_v, we need T 182 for u, T+1 for v and 2T scratch space. In all, 7T + 3 is 183 sufficient for both operations. 184 185 */ 186 187 /* Optimal choice of p seems difficult. In each iteration the division 188 * of work between hgcd and the updates of u0 and u1 depends on the 189 * current size of the u. It may be desirable to use a different 190 * choice of p in each iteration. Also the input size seems to matter; 191 * choosing p = n / 3 in the first iteration seems to improve 192 * performance slightly for input size just above the threshold, but 193 * degrade performance for larger inputs. */ 194 #define CHOOSE_P_1(n) ((n) / 2) 195 #define CHOOSE_P_2(n) ((n) / 3) 196 197 mp_size_t 198 mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep, 199 mp_ptr ap, mp_size_t an, mp_ptr bp, mp_size_t n) 200 { 201 mp_size_t talloc; 202 mp_size_t scratch; 203 mp_size_t matrix_scratch; 204 mp_size_t ualloc = n + 1; 205 206 struct gcdext_ctx ctx; 207 mp_size_t un; 208 mp_ptr u0; 209 mp_ptr u1; 210 211 mp_ptr tp; 212 213 TMP_DECL; 214 215 ASSERT (an >= n); 216 ASSERT (n > 0); 217 ASSERT (bp[n-1] > 0); 218 219 TMP_MARK; 220 221 /* FIXME: Check for small sizes first, before setting up temporary 222 storage etc. */ 223 talloc = MPN_GCDEXT_LEHMER_N_ITCH(n); 224 225 /* For initial division */ 226 scratch = an - n + 1; 227 if (scratch > talloc) 228 talloc = scratch; 229 230 if (ABOVE_THRESHOLD (n, GCDEXT_DC_THRESHOLD)) 231 { 232 /* For hgcd loop. */ 233 mp_size_t hgcd_scratch; 234 mp_size_t update_scratch; 235 mp_size_t p1 = CHOOSE_P_1 (n); 236 mp_size_t p2 = CHOOSE_P_2 (n); 237 mp_size_t min_p = MIN(p1, p2); 238 mp_size_t max_p = MAX(p1, p2); 239 matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - min_p); 240 hgcd_scratch = mpn_hgcd_itch (n - min_p); 241 update_scratch = max_p + n - 1; 242 243 scratch = matrix_scratch + MAX(hgcd_scratch, update_scratch); 244 if (scratch > talloc) 245 talloc = scratch; 246 247 /* Final mpn_gcdext_lehmer_n call. Need space for u and for 248 copies of a and b. */ 249 scratch = MPN_GCDEXT_LEHMER_N_ITCH (GCDEXT_DC_THRESHOLD) 250 + 3*GCDEXT_DC_THRESHOLD; 251 252 if (scratch > talloc) 253 talloc = scratch; 254 255 /* Cofactors u0 and u1 */ 256 talloc += 2*(n+1); 257 } 258 259 tp = TMP_ALLOC_LIMBS(talloc); 260 261 if (an > n) 262 { 263 mpn_tdiv_qr (tp, ap, 0, ap, an, bp, n); 264 265 if (mpn_zero_p (ap, n)) 266 { 267 MPN_COPY (gp, bp, n); 268 *usizep = 0; 269 TMP_FREE; 270 return n; 271 } 272 } 273 274 if (BELOW_THRESHOLD (n, GCDEXT_DC_THRESHOLD)) 275 { 276 mp_size_t gn = mpn_gcdext_lehmer_n(gp, up, usizep, ap, bp, n, tp); 277 278 TMP_FREE; 279 return gn; 280 } 281 282 MPN_ZERO (tp, 2*ualloc); 283 u0 = tp; tp += ualloc; 284 u1 = tp; tp += ualloc; 285 286 ctx.gp = gp; 287 ctx.up = up; 288 ctx.usize = usizep; 289 290 { 291 /* For the first hgcd call, there are no u updates, and it makes 292 some sense to use a different choice for p. */ 293 294 /* FIXME: We could trim use of temporary storage, since u0 and u1 295 are not used yet. For the hgcd call, we could swap in the u0 296 and u1 pointers for the relevant matrix elements. */ 297 298 struct hgcd_matrix M; 299 mp_size_t p = CHOOSE_P_1 (n); 300 mp_size_t nn; 301 302 mpn_hgcd_matrix_init (&M, n - p, tp); 303 nn = mpn_hgcd (ap + p, bp + p, n - p, &M, tp + matrix_scratch); 304 if (nn > 0) 305 { 306 ASSERT (M.n <= (n - p - 1)/2); 307 ASSERT (M.n + p <= (p + n - 1) / 2); 308 309 /* Temporary storage 2 (p + M->n) <= p + n - 1 */ 310 n = mpn_hgcd_matrix_adjust (&M, p + nn, ap, bp, p, tp + matrix_scratch); 311 312 MPN_COPY (u0, M.p[1][0], M.n); 313 MPN_COPY (u1, M.p[1][1], M.n); 314 un = M.n; 315 while ( (u0[un-1] | u1[un-1] ) == 0) 316 un--; 317 } 318 else 319 { 320 /* mpn_hgcd has failed. Then either one of a or b is very 321 small, or the difference is very small. Perform one 322 subtraction followed by one division. */ 323 u1[0] = 1; 324 325 ctx.u0 = u0; 326 ctx.u1 = u1; 327 ctx.tp = tp + n; /* ualloc */ 328 ctx.un = 1; 329 330 /* Temporary storage n */ 331 n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp); 332 if (n == 0) 333 { 334 TMP_FREE; 335 return ctx.gn; 336 } 337 338 un = ctx.un; 339 ASSERT (un < ualloc); 340 } 341 } 342 343 while (ABOVE_THRESHOLD (n, GCDEXT_DC_THRESHOLD)) 344 { 345 struct hgcd_matrix M; 346 mp_size_t p = CHOOSE_P_2 (n); 347 mp_size_t nn; 348 349 mpn_hgcd_matrix_init (&M, n - p, tp); 350 nn = mpn_hgcd (ap + p, bp + p, n - p, &M, tp + matrix_scratch); 351 if (nn > 0) 352 { 353 mp_ptr t0; 354 355 t0 = tp + matrix_scratch; 356 ASSERT (M.n <= (n - p - 1)/2); 357 ASSERT (M.n + p <= (p + n - 1) / 2); 358 359 /* Temporary storage 2 (p + M->n) <= p + n - 1 */ 360 n = mpn_hgcd_matrix_adjust (&M, p + nn, ap, bp, p, t0); 361 362 /* By the same analysis as for mpn_hgcd_matrix_mul */ 363 ASSERT (M.n + un <= ualloc); 364 365 /* FIXME: This copying could be avoided by some swapping of 366 * pointers. May need more temporary storage, though. */ 367 MPN_COPY (t0, u0, un); 368 369 /* Temporary storage ualloc */ 370 un = hgcd_mul_matrix_vector (&M, u0, t0, u1, un, t0 + un); 371 372 ASSERT (un < ualloc); 373 ASSERT ( (u0[un-1] | u1[un-1]) > 0); 374 } 375 else 376 { 377 /* mpn_hgcd has failed. Then either one of a or b is very 378 small, or the difference is very small. Perform one 379 subtraction followed by one division. */ 380 ctx.u0 = u0; 381 ctx.u1 = u1; 382 ctx.tp = tp + n; /* ualloc */ 383 ctx.un = un; 384 385 /* Temporary storage n */ 386 n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp); 387 if (n == 0) 388 { 389 TMP_FREE; 390 return ctx.gn; 391 } 392 393 un = ctx.un; 394 ASSERT (un < ualloc); 395 } 396 } 397 /* We have A = ... a + ... b 398 B = u0 a + u1 b 399 400 a = u1 A + ... B 401 b = -u0 A + ... B 402 403 with bounds 404 405 |u0|, |u1| <= B / min(a, b) 406 407 We always have u1 > 0, and u0 == 0 is possible only if u1 == 1, 408 in which case the only reduction done so far is a = A - k B for 409 some k. 410 411 Compute g = u a + v b = (u u1 - v u0) A + (...) B 412 Here, u, v are bounded by 413 414 |u| <= b, 415 |v| <= a 416 */ 417 418 ASSERT ( (ap[n-1] | bp[n-1]) > 0); 419 420 if (UNLIKELY (mpn_cmp (ap, bp, n) == 0)) 421 { 422 /* Must return the smallest cofactor, +u1 or -u0 */ 423 int c; 424 425 MPN_COPY (gp, ap, n); 426 427 MPN_CMP (c, u0, u1, un); 428 /* c == 0 can happen only when A = (2k+1) G, B = 2 G. And in 429 this case we choose the cofactor + 1, corresponding to G = A 430 - k B, rather than -1, corresponding to G = - A + (k+1) B. */ 431 ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1)); 432 if (c < 0) 433 { 434 MPN_NORMALIZE (u0, un); 435 MPN_COPY (up, u0, un); 436 *usizep = -un; 437 } 438 else 439 { 440 MPN_NORMALIZE_NOT_ZERO (u1, un); 441 MPN_COPY (up, u1, un); 442 *usizep = un; 443 } 444 445 TMP_FREE; 446 return n; 447 } 448 else if (UNLIKELY (u0[0] == 0) && un == 1) 449 { 450 mp_size_t gn; 451 ASSERT (u1[0] == 1); 452 453 /* g = u a + v b = (u u1 - v u0) A + (...) B = u A + (...) B */ 454 gn = mpn_gcdext_lehmer_n (gp, up, usizep, ap, bp, n, tp); 455 456 TMP_FREE; 457 return gn; 458 } 459 else 460 { 461 mp_size_t u0n; 462 mp_size_t u1n; 463 mp_size_t lehmer_un; 464 mp_size_t lehmer_vn; 465 mp_size_t gn; 466 467 mp_ptr lehmer_up; 468 mp_ptr lehmer_vp; 469 int negate; 470 471 lehmer_up = tp; tp += n; 472 473 /* Call mpn_gcdext_lehmer_n with copies of a and b. */ 474 MPN_COPY (tp, ap, n); 475 MPN_COPY (tp + n, bp, n); 476 gn = mpn_gcdext_lehmer_n (gp, lehmer_up, &lehmer_un, tp, tp + n, n, tp + 2*n); 477 478 u0n = un; 479 MPN_NORMALIZE (u0, u0n); 480 ASSERT (u0n > 0); 481 482 if (lehmer_un == 0) 483 { 484 /* u == 0 ==> v = g / b == 1 ==> g = - u0 A + (...) B */ 485 MPN_COPY (up, u0, u0n); 486 *usizep = -u0n; 487 488 TMP_FREE; 489 return gn; 490 } 491 492 lehmer_vp = tp; 493 /* Compute v = (g - u a) / b */ 494 lehmer_vn = compute_v (lehmer_vp, 495 ap, bp, n, gp, gn, lehmer_up, lehmer_un, tp + n + 1); 496 497 if (lehmer_un > 0) 498 negate = 0; 499 else 500 { 501 lehmer_un = -lehmer_un; 502 negate = 1; 503 } 504 505 u1n = un; 506 MPN_NORMALIZE (u1, u1n); 507 ASSERT (u1n > 0); 508 509 ASSERT (lehmer_un + u1n <= ualloc); 510 ASSERT (lehmer_vn + u0n <= ualloc); 511 512 /* We may still have v == 0 */ 513 514 /* Compute u u0 */ 515 if (lehmer_un <= u1n) 516 /* Should be the common case */ 517 mpn_mul (up, u1, u1n, lehmer_up, lehmer_un); 518 else 519 mpn_mul (up, lehmer_up, lehmer_un, u1, u1n); 520 521 un = u1n + lehmer_un; 522 un -= (up[un - 1] == 0); 523 524 if (lehmer_vn > 0) 525 { 526 mp_limb_t cy; 527 528 /* Overwrites old u1 value */ 529 if (lehmer_vn <= u0n) 530 /* Should be the common case */ 531 mpn_mul (u1, u0, u0n, lehmer_vp, lehmer_vn); 532 else 533 mpn_mul (u1, lehmer_vp, lehmer_vn, u0, u0n); 534 535 u1n = u0n + lehmer_vn; 536 u1n -= (u1[u1n - 1] == 0); 537 538 if (u1n <= un) 539 { 540 cy = mpn_add (up, up, un, u1, u1n); 541 } 542 else 543 { 544 cy = mpn_add (up, u1, u1n, up, un); 545 un = u1n; 546 } 547 up[un] = cy; 548 un += (cy != 0); 549 550 ASSERT (un < ualloc); 551 } 552 *usizep = negate ? -un : un; 553 554 TMP_FREE; 555 return gn; 556 } 557 } 558
Indexes created Sat Feb 28 05:31:39 UTC 2026