X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/da83a5610ddde70bfae47b83707f6e80bce07467..578a86d91941a0f722b87973d88e84ec2cf9a608:/mp-gcd.c

diff --git a/mp-gcd.c b/mp-gcd.c
index 19192f5..5e663a1 100644
--- a/mp-gcd.c
+++ b/mp-gcd.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mp-gcd.c,v 1.2 1999/11/22 20:49:56 mdw Exp $
+ * $Id: mp-gcd.c,v 1.7 2004/04/08 01:36:15 mdw Exp $
  *
  * Extended GCD calculation
  *
@@ -27,18 +27,6 @@
  * MA 02111-1307, USA.
  */
 
-/*----- Revision history --------------------------------------------------* 
- *
- * $Log: mp-gcd.c,v $
- * Revision 1.2  1999/11/22 20:49:56  mdw
- * Fix bug which failed to favour `x' when `y' wasn't wanted and the two
- * arguments needed swapping.
- *
- * Revision 1.1  1999/11/17 18:02:16  mdw
- * New multiprecision integer arithmetic suite.
- *
- */
-
 /*----- Header files ------------------------------------------------------*/
 
 #include "mp.h"
@@ -54,198 +42,214 @@
  *
  * Use:		Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
  *		@ax + by = gcd(a, b)@.  This is useful for computing modular
- *		inverses.  Neither @a@ nor @b@ may be zero.  Note that,
- *		unlike @mp_div@ for example, it is not possible to specify
- *		explicit destinations -- new MPs are always allocated.
+ *		inverses.
  */
 
 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;
-  int ext = xx || yy;
-  int swap = 0;
+  mp *q = MP_NEW;
+  unsigned f = 0;
 
-  /* --- Ensure that @a@ is larger than @b@ --- */
+#define f_swap 1u
+#define f_aneg 2u
+#define f_bneg 4u
+#define f_ext 8u
 
-  if (MP_CMP(a, <, b)) {
-    { mp *t = a; a = b; b = t; }
-    swap = 1;
-  }
+  /* --- Sort out some initial flags --- */
 
-  /* --- Take a reference to the arguments --- */
+  if (xx || yy)
+    f |= f_ext;
 
-  a = MP_COPY(a);
-  b = MP_COPY(b);
+  if (a->f & MP_NEG)
+    f |= f_aneg;
+  if (b->f & MP_NEG)
+    f |= f_bneg;
 
-  /* --- Make sure @a@ and @b@ are not both even --- */
+  /* --- Ensure that @a@ is larger than @b@ --- *
+   *
+   * Use absolute values here!
+   */
 
-  if (((a->v[0] | b->v[0]) & 1) == 0) {
-    mpscan asc, bsc;
+  if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
+    { mp *t = a; a = b; b = t; }
+    f |= f_swap;
+  }
 
-    /* --- Break off my copies --- */
+  /* --- Check for zeroness --- */
 
-    MP_SPLIT(a);
-    MP_SPLIT(b);
-    MP_SCAN(&asc, a);
-    MP_SCAN(&bsc, b);
+  if (MP_EQ(b, MP_ZERO)) {
 
-    /* --- Start scanning --- */
+    /* --- Store %$|a|$% as the GCD --- */
 
-    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++;
+    if (gcd) {
+      if (*gcd) MP_DROP(*gcd);
+      a = MP_COPY(a);
+      if (a->f & MP_NEG) {
+	MP_SPLIT(a);
+	a->f &= ~MP_NEG;
+	f |= f_aneg;
+      }
+      *gcd = a;
     }
 
-    /* --- Shift @a@ and @b@ down --- */
+    /* --- Store %$1$% and %$0$% in the appropriate bins --- */
 
-    a = mp_lsr(a, a, shift);
-    b = mp_lsr(b, b, shift);
+    if (f & f_ext) {
+      if (f & f_swap) {
+	mp **t = xx; xx = yy; yy = t;
+      }
+      if (xx) {
+	if (*xx) MP_DROP(*xx);
+	if (MP_EQ(a, MP_ZERO))
+	  *xx = MP_ZERO;
+	else if (f & f_aneg)
+	  *xx = MP_MONE;
+	else
+	  *xx = MP_ONE;
+      }
+      if (yy) {
+	if (*yy) MP_DROP(*yy);
+	*yy = MP_ZERO;
+      }
+    }
+    return;
   }
 
-  /* --- Set up @u@ and @v@ --- */
-
-  u = MP_COPY(a);
-  v = MP_COPY(b);
+  /* --- Take a reference to the arguments --- */
 
-  /* --- 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 (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++;
-      }
+  a = MP_COPY(a);
+  b = MP_COPY(b);
 
-      if (nn) {
-	u = mp_lsr(u, u, nn);
-	if (ext && n) {
-	  X = mp_lsr(X, X, n);
-	  Y = mp_lsr(Y, Y, n);
-	}
-      }      
-    }
+  /* --- Make sure @a@ and @b@ are not both even --- */
 
-    /* --- 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 (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++;
-      }
+  MP_SPLIT(a); a->f &= ~MP_NEG;
+  MP_SPLIT(b); b->f &= ~MP_NEG;
 
-      if (nn) {
-	v = mp_lsr(v, v, nn);
-	if (ext && n) {
-	  x = mp_lsr(x, x, n);
-	  y = mp_lsr(y, y, n);
-	}
-      }      
-    }
-
-    /* --- End-of-loop fiddling --- */
+  u = MP_COPY(a);
+  v = MP_COPY(b);
 
-    if (MP_CMP(u, >=, v)) {
-      u = mp_sub(u, u, v);
-      if (ext) {
-	X = mp_sub(X, X, x);
-	Y = mp_sub(Y, Y, y);
-      }
-    } else {
-      v = mp_sub(v, v, u);
-      if (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 --- */
-
-  if (gcd)
-    *gcd = mp_lsl(v, v, shift);
-  else
-    MP_DROP(v);
+  MP_DROP(q);
+  if (!gcd)
+    MP_DROP(u);
+  else {
+    if (*gcd) MP_DROP(*gcd);
+    u->f &= ~MP_NEG;
+    *gcd = u;
+  }
 
   /* --- Perform a little normalization --- *
    *
    * Ensure that the coefficient returned is positive, if there is only one.
-   * If there are two, favour @y@.
+   * If there are two, favour @y@.  Of course, if the original arguments were
+   * negative then I'll need to twiddle their signs as well.
    */
 
-  if (ext) {
-    if (swap) {
+  if (f & f_ext) {
+
+    /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
+
+    if (f & f_swap) {
       mp *t = x; x = y; y = t;
       t = a; a = b; b = t;
     }
+
+    /* --- Sort out the signs --- *
+     *
+     * Note that %$ax + by = a(x - b) + b(y + a)$%.
+     *
+     * This is currently bodgy.  It needs sorting out at some time.
+     */
+
     if (yy) {
       if (y->f & MP_NEG) {
-	y = mp_add(y, y, a);
-	x = mp_sub(x, x, b);
+	do {
+	  y = mp_add(y, y, a);
+	  x = mp_sub(x, x, b);
+	} while (y->f & MP_NEG);
+      } else {
+	while (MP_CMP(y, >=, a)) {
+	  y = mp_sub(y, y, a);
+	  x = mp_add(x, x, b);
+	}
       }
-    } else if (x->f & MP_NEG)
-      x = mp_add(x, x, b);
+    } else {
+      if (x->f & MP_NEG) {
+	do
+	  x = mp_add(x, x, b);
+	while (x->f & MP_NEG);
+      } else {
+	while (MP_CMP(x, >=, b))
+	  x = mp_sub(x, x, b);
+      }
+    }
+
+    /* --- Twiddle the signs --- */
+
+    if (f & f_aneg)
+      x->f ^= MP_NEG;
+    if (f & f_bneg)
+      y->f ^= MP_NEG;
 
-    if (xx) *xx = x; else MP_DROP(x);
-    if (yy) *yy = y; else MP_DROP(y);
+    /* --- Store the results --- */
+
+    if (!xx)
+      MP_DROP(x);
+    else {
+      if (*xx) MP_DROP(*xx);
+      *xx = x;
+    }
+
+    if (!yy)
+      MP_DROP(y);
+    else {
+      if (*yy) MP_DROP(*yy);
+      *yy = y;
+    }
   }
 
-  MP_DROP(u);
+  MP_DROP(v);
   MP_DROP(X); MP_DROP(Y);
   MP_DROP(a); MP_DROP(b);
 }
 
+/* -- @mp_modinv@ --- *
+ *
+ * Arguments:	@mp *d@ = destination
+ *		@mp *x@ = argument
+ *		@mp *p@ = modulus
+ *
+ * Returns:	The inverse %$x^{-1} \bmod p$%.
+ *
+ * Use:		Computes a modular inverse.    An assertion fails if %$p$%
+ *		has no inverse.
+ */
+
+mp *mp_modinv(mp *d, mp *x, mp *p)
+{
+  mp *g = MP_NEW;
+  mp_gcd(&g, 0, &d, p, x);
+  assert(MP_EQ(g, MP_ONE));
+  mp_drop(g);
+  return (d);
+}
+
 /*----- Test rig ----------------------------------------------------------*/
 
 #ifdef TEST_RIG
@@ -259,9 +263,9 @@ static int gcd(dstr *v)
   mp *x = *(mp **)v[3].buf;
   mp *y = *(mp **)v[4].buf;
 
-  mp *gg, *xx, *yy;
+  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);
@@ -270,7 +274,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);
@@ -284,13 +288,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);
@@ -301,6 +305,7 @@ static int gcd(dstr *v)
   }
   MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
   MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }