From b2253c70111a59b034cb53680689ec002f7ab126 Mon Sep 17 00:00:00 2001 From: mdw Date: Sun, 8 Oct 2000 12:02:41 +0000 Subject: [PATCH] Use Euclid's algorithm rather than the binary one. --- mp-gcd.c | 165 +++++++++++---------------------------------------------------- 1 file changed, 29 insertions(+), 136 deletions(-) diff --git a/mp-gcd.c b/mp-gcd.c index 38f95d4..a17a502 100644 --- a/mp-gcd.c +++ b/mp-gcd.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mp-gcd.c,v 1.4 2000/06/17 11:34:46 mdw Exp $ + * $Id: mp-gcd.c,v 1.5 2000/10/08 12:02:41 mdw Exp $ * * Extended GCD calculation * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: mp-gcd.c,v $ + * Revision 1.5 2000/10/08 12:02:41 mdw + * Use Euclid's algorithm rather than the binary one. + * * Revision 1.4 2000/06/17 11:34:46 mdw * More hacking for the signs of the outputs. * @@ -65,10 +68,10 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) { - mp *X = MP_ONE, *Y = MP_ZERO; - mp *x = MP_ZERO, *y = MP_ONE; + mp *x = MP_ONE, *X = MP_ZERO; + mp *y = MP_ZERO, *Y = MP_ONE; mp *u, *v; - size_t shift = 0; + mp *q = MP_NEW; unsigned f = 0; enum { @@ -100,7 +103,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) /* --- Check for zeroness --- */ - if (MP_CMP(b, ==, MP_ZERO)) { + if (MP_EQ(b, MP_ZERO)) { /* --- Store %$|a|$% as the GCD --- */ @@ -123,7 +126,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) } if (xx) { if (*xx) MP_DROP(*xx); - if (MP_CMP(a, ==, MP_ZERO)) + if (MP_EQ(a, MP_ZERO)) *xx = MP_ZERO; else if (f & f_aneg) *xx = MP_MONE; @@ -148,140 +151,30 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) MP_SPLIT(a); a->f &= ~MP_NEG; MP_SPLIT(b); b->f &= ~MP_NEG; - if (((a->v[0] | b->v[0]) & 1) == 0) { - mpscan asc, bsc; - - /* --- Break off my copies --- */ - - MP_SCAN(&asc, a); - MP_SCAN(&bsc, b); - - /* --- Start scanning --- */ - - for (;;) { - if (!MP_STEP(&asc) || !MP_STEP(&bsc)) - assert(((void)"zero argument passed to mp_gcd", 0)); - if (MP_BIT(&asc) || MP_BIT(&bsc)) - break; - shift++; - } - - /* --- Shift @a@ and @b@ down --- */ - - a = mp_lsr(a, a, shift); - b = mp_lsr(b, b, shift); - } - - /* --- Set up @u@ and @v@ --- */ - u = MP_COPY(a); v = MP_COPY(b); - /* --- Start the main loop --- */ - - for (;;) { - - /* --- While @u@ is even --- */ - - { - mpscan sc, xsc, ysc; - size_t n = 0, nn = 0; - - MP_SCAN(&sc, u); - MP_SCAN(&xsc, X); MP_SCAN(&ysc, Y); - for (;;) { - MP_STEP(&sc); - MP_STEP(&xsc); MP_STEP(&ysc); - if (MP_BIT(&sc)) - break; - if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) { - if (n) { - X = mp_lsr(X, X, n); - Y = mp_lsr(Y, Y, n); - n = 0; - } - X = mp_add(X, X, b); - Y = mp_sub(Y, Y, a); - MP_SCAN(&xsc, X); - MP_SCAN(&ysc, Y); - MP_STEP(&xsc); MP_STEP(&ysc); - } - n++; nn++; - } - - if (nn) { - u = mp_lsr(u, u, nn); - if ((f & f_ext) && n) { - X = mp_lsr(X, X, n); - Y = mp_lsr(Y, Y, n); - } - } - } - - /* --- While @v@ is even --- */ - - { - mpscan sc, xsc, ysc; - size_t n = 0, nn = 0; - - MP_SCAN(&sc, v); - MP_SCAN(&xsc, x); MP_SCAN(&ysc, y); - for (;;) { - MP_STEP(&sc); - MP_STEP(&xsc); MP_STEP(&ysc); - if (MP_BIT(&sc)) - break; - if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) { - if (n) { - x = mp_lsr(x, x, n); - y = mp_lsr(y, y, n); - n = 0; - } - x = mp_add(x, x, b); - y = mp_sub(y, y, a); - MP_SCAN(&xsc, x); - MP_SCAN(&ysc, y); - MP_STEP(&xsc); MP_STEP(&ysc); - } - n++; nn++; - } - - if (nn) { - v = mp_lsr(v, v, nn); - if ((f & f_ext) && n) { - x = mp_lsr(x, x, n); - y = mp_lsr(y, y, n); - } - } - } - - /* --- End-of-loop fiddling --- */ - - if (MP_CMP(u, >=, v)) { - u = mp_sub(u, u, v); - if (f & f_ext) { - X = mp_sub(X, X, x); - Y = mp_sub(Y, Y, y); - } - } else { - v = mp_sub(v, v, u); - if (f & f_ext) { - x = mp_sub(x, x, X); - y = mp_sub(y, y, Y); - } + while (MP_LEN(v)) { + mp *t; + mp_div(&q, &u, u, v); + if (f & f_ext) { + t = mp_mul(MP_NEW, X, q); + t = mp_sub(t, x, t); + MP_DROP(x); x = X; X = t; + t = mp_mul(MP_NEW, Y, q); + t = mp_sub(t, y, t); + MP_DROP(y); y = Y; Y = t; } - - if (MP_CMP(u, ==, MP_ZERO)) - break; + t = u; u = v; v = t; } - /* --- Write the results out --- */ - + MP_DROP(q); if (!gcd) - MP_DROP(v); + MP_DROP(u); else { if (*gcd) MP_DROP(*gcd); - *gcd = mp_lsl(v, v, shift); + u->f &= ~MP_NEG; + *gcd = u; } /* --- Perform a little normalization --- * @@ -354,7 +247,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) } } - MP_DROP(u); + MP_DROP(v); MP_DROP(X); MP_DROP(Y); MP_DROP(a); MP_DROP(b); } @@ -374,7 +267,7 @@ static int gcd(dstr *v) mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW; mp_gcd(&gg, &xx, &yy, a, b); - if (MP_CMP(x, !=, xx)) { + if (!MP_EQ(x, xx)) { fputs("\n*** mp_gcd(x) failed", stderr); fputs("\na = ", stderr); mp_writefile(a, stderr, 10); fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); @@ -383,7 +276,7 @@ static int gcd(dstr *v) fputc('\n', stderr); ok = 0; } - if (MP_CMP(y, !=, yy)) { + if (!MP_EQ(y, yy)) { fputs("\n*** mp_gcd(y) failed", stderr); fputs("\na = ", stderr); mp_writefile(a, stderr, 10); fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); @@ -397,13 +290,13 @@ static int gcd(dstr *v) mp *ax = mp_mul(MP_NEW, a, xx); mp *by = mp_mul(MP_NEW, b, yy); ax = mp_add(ax, ax, by); - if (MP_CMP(ax, ==, gg)) + if (MP_EQ(ax, gg)) fputs("\n*** (Alternative result found.)\n", stderr); MP_DROP(ax); MP_DROP(by); } - if (MP_CMP(g, !=, gg)) { + if (!MP_EQ(g, gg)) { fputs("\n*** mp_gcd(gcd) failed", stderr); fputs("\na = ", stderr); mp_writefile(a, stderr, 10); fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); -- 2.11.0