diff options
| -rw-r--r-- | ChangeLog | 10 | ||||
| -rw-r--r-- | mpn/generic/gcd.c | 4 | ||||
| -rw-r--r-- | mpn/generic/gcd_subdiv_step.c | 84 | ||||
| -rw-r--r-- | mpn/generic/gcdext.c | 6 | ||||
| -rw-r--r-- | mpn/generic/gcdext_lehmer.c | 5 |
5 files changed, 82 insertions, 27 deletions
@@ -1,3 +1,13 @@ +2011-05-19 Niels Möller <nisse@lysator.liu.se> + + * mpn/generic/gcd_subdiv_step.c (mpn_gcd_subdiv_step): New + argument s, for use by hgcd. + + * mpn/generic/gcd.c (mpn_gcd): Pass s = 0 to mpn_gcd_subdiv_step. + * mpn/generic/gcdext.c (mpn_gcdext): Likewise. Also added an ASSERT. + * mpn/generic/gcdext_lehmer.c (mpn_gcdext_lehmer_n): Likewise. + (mpn_gcdext_hook): Added some ASSERTs. + 2011-05-17 Niels Möller <nisse@lysator.liu.se> * doc/gmp.texi (mpn_gcd, mpn_gcdext): Document input requirements: diff --git a/mpn/generic/gcd.c b/mpn/generic/gcd.c index 055ed90f8..77d85c89b 100644 --- a/mpn/generic/gcd.c +++ b/mpn/generic/gcd.c @@ -210,7 +210,7 @@ mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n) else { /* Temporary storage n */ - n = mpn_gcd_subdiv_step (up, vp, n, gcd_hook, &ctx, tp); + n = mpn_gcd_subdiv_step (up, vp, n, 0, gcd_hook, &ctx, tp); if (n == 0) goto done; } @@ -254,7 +254,7 @@ mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n) subtraction followed by one division. */ /* Temporary storage n */ - n = mpn_gcd_subdiv_step (up, vp, n, &gcd_hook, &ctx, tp); + n = mpn_gcd_subdiv_step (up, vp, n, 0, &gcd_hook, &ctx, tp); if (n == 0) goto done; } diff --git a/mpn/generic/gcd_subdiv_step.c b/mpn/generic/gcd_subdiv_step.c index 00d32db10..63514f24f 100644 --- a/mpn/generic/gcd_subdiv_step.c +++ b/mpn/generic/gcd_subdiv_step.c @@ -29,8 +29,14 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or b is small, or the difference is small. Perform one subtraction - followed by one division. Returns zero if the gcd is found. - Otherwise, compute the reduced a and b, and return the new size. */ + followed by one division. The normal case is to compute the reduced + a and b, and return the new size. + + If s == 0 (used for gcd and gcdext), returns zero if the gcd is + found. + + If s > 0, don't reduce to size <= s, and return zero if no + reduction is possible (if either a, b or |a-b| is of size <= s). */ /* The hook function is called as @@ -38,15 +44,15 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ in the following cases: - + If A == B at the start, G is the gcd, Q is NULL, d = -1. + + If A = B at the start, G is the gcd, Q is NULL, d = -1. + If one input is zero at the start, G is the gcd, Q is NULL, - d = 0 if A == G and d = 1 if B == G. + d = 0 if A = G and d = 1 if B = G. Otherwise, if d = 0 we have just subtracted a multiple of A from B, and if d = 1 we have subtracted a multiple of B from A. - + If A == B after subtraction, G is the gcd, Q is NULL. + + If A = B after subtraction, G is the gcd, Q is NULL. + If we get a zero remainder after division, G is the gcd, Q is the quotient. @@ -56,7 +62,7 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ */ mp_size_t -mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, +mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t s, gcd_subdiv_step_hook *hook, void *ctx, mp_ptr tp) { @@ -82,8 +88,10 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, MPN_CMP (c, ap, bp, an); if (UNLIKELY (c == 0)) { - /* For gcdext, return the smallest of the two cofactors. */ - hook (ctx, ap, an, NULL, 0, -1); + /* For gcdext, return the smallest of the two cofactors, so + pass d = -1. */ + if (s == 0) + hook (ctx, ap, an, NULL, 0, -1); return 0; } else if (c > 0) @@ -99,16 +107,27 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, MPN_PTR_SWAP (ap, an, bp, bn); swapped ^= 1; } - if (an == 0) - { - hook (ctx, bp, bn, NULL, 0, swapped ^ 1); - return 0; - } } + if (an <= s) + { + if (s == 0) + hook (ctx, bp, bn, NULL, 0, swapped ^ 1); + return 0; + } + ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an)); MPN_NORMALIZE (bp, bn); ASSERT (bn > 0); + if (bn <= s) + { + /* Undo subtraction. */ + mp_limb_t cy = mpn_add (bp, ap, an, bp, bn); + if (cy > 0) + bp[an] = cy; + return 0; + } + /* Arrange so that a < b */ if (an == bn) { @@ -116,7 +135,12 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, MPN_CMP (c, ap, bp, an); if (UNLIKELY (c == 0)) { - hook (ctx, bp, bn, NULL, 0, swapped); + if (s > 0) + /* Just record subtraction and return */ + hook (ctx, NULL, 0, &one, 1, swapped); + else + /* Found gcd. */ + hook (ctx, bp, bn, NULL, 0, swapped); return 0; } @@ -141,14 +165,30 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an); qn = bn - an + 1; - if (mpn_zero_p (bp, an)) - { - hook (ctx, ap, an, tp, qn, swapped); - return 0; - } - else + bn = an; + MPN_NORMALIZE (bp, bn); + + if (UNLIKELY (bn <= s)) { - hook (ctx, NULL, 0, tp, qn, swapped); - return an; + if (s == 0) + { + hook (ctx, ap, an, tp, qn, swapped); + return 0; + } + + /* Quotient is one too large, so decrement it and add back A. */ + if (bn > 0) + { + mp_limb_t cy = mpn_add (bp, ap, an, bp, bn); + if (cy) + bp[an++] = cy; + } + else + MPN_COPY (bp, ap, an); + + MPN_DECR_U (tp, qn, 1); } + + hook (ctx, NULL, 0, tp, qn, swapped); + return an; } diff --git a/mpn/generic/gcdext.c b/mpn/generic/gcdext.c index 88fc47f94..402ef9d41 100644 --- a/mpn/generic/gcdext.c +++ b/mpn/generic/gcdext.c @@ -319,7 +319,7 @@ mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep, ctx.un = 1; /* Temporary storage n */ - n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp); + n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp); if (n == 0) { TMP_FREE; @@ -374,7 +374,7 @@ mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep, ctx.un = un; /* Temporary storage n */ - n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp); + n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp); if (n == 0) { TMP_FREE; @@ -386,6 +386,8 @@ mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep, } } + ASSERT ( (ap[n-1] | bp[n-1]) > 0); + if (UNLIKELY (mpn_cmp (ap, bp, n) == 0)) { /* Must return the smallest cofactor, +u1 or -u0 */ diff --git a/mpn/generic/gcdext_lehmer.c b/mpn/generic/gcdext_lehmer.c index 904f5cc08..911c53a5b 100644 --- a/mpn/generic/gcdext_lehmer.c +++ b/mpn/generic/gcdext_lehmer.c @@ -35,6 +35,9 @@ mpn_gcdext_hook (void *p, mp_srcptr gp, mp_size_t gn, { mp_srcptr up; + ASSERT (gn > 0); + ASSERT (gp[gn-1] > 0); + MPN_COPY (ctx->gp, gp, gn); ctx->gn = gn; @@ -216,7 +219,7 @@ mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize, /* Temporary storage n for the quotient and ualloc for the new cofactor. */ - n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp); + n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp); if (n == 0) return ctx.gn; |