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.

1 files changed, 75 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);
 }