gcd: rewrite for better efficiency.

* arith.c (gcd): New implementation which uses arithmetic
in the unsigned type ucnum if both operands are in that
type's range. This uses Stein's algorithm a.k.a.
binary GCD. The mpi_gcd function is used only if at least
one argument is a bignum whose value doesn't fit into
a ucnum.

* tests/016/arith.tl: gcd test cases added.

2 files changed, 104 insertions, 24 deletions

diff --git a/arith.c b/arith.c
index 1159de96..05edf13f 100644
--- a/arith.c
+++ b/arith.c
@@ -2603,33 +2603,84 @@ val square(val anum)
 
 val gcd(val anum, val bnum)
 {
-  val n;
-
-  if (!integerp(anum) || !integerp(bnum))
-    goto inval;
-
-  if (anum == zero)
-    return bnum;
-
-  if (bnum == zero)
-    return anum;
+  val self = lit("gcd");
+  ucnum ua, ub;
 
-  if (fixnump(anum))
-    anum = bignum(c_n(anum));
-
-  if (fixnump(bnum))
-    bnum = bignum(c_n(bnum));
-
-  n = make_bignum();
+  switch (TYPE_PAIR(type(anum), type(bnum))) {
+  case TYPE_PAIR(BGNUM, BGNUM):
+    if (mp_in_uintptr_range(mp(anum)) && mp_in_uintptr_range(mp(bnum))) {
+      ua = c_unum(anum, self);
+      ub = c_unum(bnum, self);
+      goto both_ucnum;
+    } else {
+      val n = make_bignum();
+      if (mp_gcd(mp(anum), mp(bnum), mp(n)) != MP_OKAY)
+        break;
+      return normalize(n);
+    }
+  case TYPE_PAIR(NUM, NUM):
+    if (anum == zero)
+      return bnum;
+    if (bnum == zero)
+      return anum;
+    ua = c_u(anum);
+    ub = c_u(bnum);
+  both_ucnum:
+    {
+      int k = 0;
+      while ((ua & 1) == 0 && (ub & 1) == 0) {
+        ua >>= 1;
+        ub >>= 1;
+        k++;
+      }
+      while ((ub & 1) == 0)
+        ub >>= 1;
+      while ((ua & 1) == 0)
+        ua >>= 1;
+      for (;;) {
+        if (ua > ub) {
+          ucnum tmp = ua;
+          ua = ub;
+          ub = tmp;
+        }
+        ub -= ua;
+        if (ub == 0)
+          break;
+        while ((ub & 1) == 0)
+          ub >>= 1;
+      }
+      return unum(ua << k);
+    }
+  case TYPE_PAIR(NUM, BGNUM):
+    {
+      val tmp = anum;
+      anum = bnum;
+      bnum = tmp;
+    }
+    /* fallthrough */
+  case TYPE_PAIR(BGNUM, NUM):
+    if (mp_in_uintptr_range(mp(anum))) {
+      ua = c_unum(anum, self);
+      ub = c_u(bnum);
+      goto both_ucnum;
+    } else {
+      mp_int bn;
+      val n = make_bignum();
 
-  if (mp_gcd(mp(anum), mp(bnum), mp(n)) != MP_OKAY)
-    goto bad;
+      mp_init(&bn);
+      mp_set_intptr(&bn, c_u(bnum));
+      if (mp_gcd(mp(anum), &bn, mp(n)) != MP_OKAY) {
+        mp_clear(&bn);
+        break;
+      }
+      mp_clear(&bn);
+      return normalize(n);
+    }
+  default:
+    uw_throwf(error_s, lit("gcd: non-integral operands ~s ~s"),
+              anum, bnum, nao);
+  }
 
-  return normalize(n);
-inval:
-  uw_throwf(error_s, lit("gcd: non-integral operands ~s ~s"),
-            anum, bnum, nao);
-bad:
   uw_throwf(error_s, lit("gcd: operation failed on ~s ~s"),
             anum, bnum, nao);
 }
diff --git a/tests/016/arith.tl b/tests/016/arith.tl
index d740835b..8f401463 100644
--- a/tests/016/arith.tl
+++ b/tests/016/arith.tl
@@ -372,3 +372,32 @@
   (some-false ((a '(1 2 3)) (b '(4 5 6))) (> a b)) t
   (some-false ((a '(1 2 3)) (b '(4 0 6))) (> a b)) t
   (some-false ((a '(1 2 3)) (b '(0 1 2))) (> a b)) nil)
+
+
+(mvtest
+  (gcd 0 0) 0
+  (gcd 0 1) 1
+  (gcd 1 0) 1
+  (gcd 100 0) 100
+  (gcd 0 100) 100
+  (gcd 0 (expt 10 60)) (expt 10 60)
+  (gcd (expt 10 60) 0) (expt 10 60))
+
+(defun power-set (s)
+  (mappend* (op comb s) (range 0 (len s))))
+
+(defun gcd-grind (primes)
+  (each-prod ((lp (cdr (power-set primes)))
+              (rp (cdr (power-set primes))))
+    (let ((ip (isec lp rp)))
+      (vtest (gcd (* . lp) (* . rp)) (* . ip)))))
+
+(each ((x 0..64)
+       (y 0..64))
+  (vtest (gcd (ash 1 x) (ash 1 y)) (ash 1 (min x y)))
+  (vtest (gcd (ash 3 x) (ash 5 y)) (ash 1 (min x y)))
+  (vtest (gcd (ash 6 x) (ash 15 y)) (ash 3 (min x y))))
+
+(gcd-grind '(2 3 5 7 11 13 17 19 23))
+
+(gcd-grind '(2 3 5 4294967291 4294967311 4294967357 4294967371))