[Guile-commits] 19/85: Implement gcd with new integer lib

[Andy Wingo]

wingo pushed a commit to branch main
in repository guile.

commit 1e0797db7289ad27758a1cf77d89105459475ee7
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Sun Dec 19 08:57:06 2021 +0100

    Implement gcd with new integer lib
    
    * libguile/integers.c (scm_integer_gcd_ii)
    (scm_integer_gcd_zi, scm_integer_gcd_zz): New internal functions.
    * libguile/integers.h: Declare internal functions.
    * libguile/numbers.c (scm_gcd): Use the new functions.
---
 libguile/integers.c | 75 +++++++++++++++++++++++++++++++++++++++++++++
 libguile/integers.h |  4 +++
 libguile/numbers.c  | 88 +++++------------------------------------------------
 3 files changed, 86 insertions(+), 81 deletions(-)

diff --git a/libguile/integers.c b/libguile/integers.c
index 45f769f56..cca40bc8b 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -1742,3 +1742,78 @@ scm_integer_round_divide_zz (SCM x, SCM y, SCM *qp, SCM 
*rp)
 {
   integer_round_divide_zz (scm_bignum (x), scm_bignum (y), qp, rp);
 }
+
+SCM
+scm_integer_gcd_ii (scm_t_inum x, scm_t_inum y)
+{
+  scm_t_inum u = x < 0 ? -x : x;
+  scm_t_inum v = y < 0 ? -y : y;
+  scm_t_inum result;
+  if (x == 0)
+    result = v;
+  else if (y == 0)
+    result = u;
+  else
+    {
+      int k = 0;
+      /* Determine a common factor 2^k */
+      while (((u | v) & 1) == 0)
+        {
+          k++;
+          u >>= 1;
+          v >>= 1;
+        }
+      /* Now, any factor 2^n can be eliminated */
+      if ((u & 1) == 0)
+        while ((u & 1) == 0)
+          u >>= 1;
+      else
+        while ((v & 1) == 0)
+          v >>= 1;
+      /* Both u and v are now odd.  Subtract the smaller one
+         from the larger one to produce an even number, remove
+         more factors of two, and repeat. */
+      while (u != v)
+        {
+          if (u > v)
+            {
+              u -= v;
+              while ((u & 1) == 0)
+                u >>= 1;
+            }
+          else
+            {
+              v -= u;
+              while ((v & 1) == 0)
+                v >>= 1;
+            }
+        }
+      result = u << k;
+    }
+  return ulong_to_scm (result);
+}
+
+SCM
+scm_integer_gcd_zi (SCM x, scm_t_inum y)
+{
+  scm_t_bits result;
+  if (y == 0)
+    return scm_abs (x);
+  if (y < 0)
+    y = -y;
+  result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), y);
+  scm_remember_upto_here_1 (x);
+  return ulong_to_scm (result);
+}
+
+SCM
+scm_integer_gcd_zz (SCM x, SCM y)
+{
+  mpz_t result, zx, zy;
+  mpz_init (result);
+  alias_bignum_to_mpz (scm_bignum (x), zx);
+  alias_bignum_to_mpz (scm_bignum (y), zy);
+  mpz_gcd (result, zx, zy);
+  scm_remember_upto_here_2 (x, y);
+  return take_mpz (result);
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index e1193c231..699525b3f 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -124,6 +124,10 @@ SCM_INTERNAL void scm_integer_round_divide_zi (SCM x, 
scm_t_inum y,
 SCM_INTERNAL void scm_integer_round_divide_zz (SCM x, SCM y,
                                                SCM *qp, SCM *rp);
 
+SCM_INTERNAL SCM scm_integer_gcd_ii (scm_t_inum x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_gcd_zi (SCM x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_gcd_zz (SCM x, SCM y);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 32dad4d01..f0ddc9379 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -2782,68 +2782,15 @@ SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1,
 SCM
 scm_gcd (SCM x, SCM y)
 {
-  if (SCM_UNLIKELY (SCM_UNBNDP (y)))
+  if (SCM_UNBNDP (y))
     return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x);
   
-  if (SCM_LIKELY (SCM_I_INUMP (x)))
+  if (SCM_I_INUMP (x))
     {
-      if (SCM_LIKELY (SCM_I_INUMP (y)))
-        {
-          scm_t_inum xx = SCM_I_INUM (x);
-          scm_t_inum yy = SCM_I_INUM (y);
-          scm_t_inum u = xx < 0 ? -xx : xx;
-          scm_t_inum v = yy < 0 ? -yy : yy;
-          scm_t_inum result;
-          if (SCM_UNLIKELY (xx == 0))
-           result = v;
-         else if (SCM_UNLIKELY (yy == 0))
-           result = u;
-         else
-           {
-             int k = 0;
-             /* Determine a common factor 2^k */
-             while (((u | v) & 1) == 0)
-               {
-                 k++;
-                 u >>= 1;
-                 v >>= 1;
-               }
-             /* Now, any factor 2^n can be eliminated */
-             if ((u & 1) == 0)
-               while ((u & 1) == 0)
-                 u >>= 1;
-             else
-               while ((v & 1) == 0)
-                 v >>= 1;
-             /* Both u and v are now odd.  Subtract the smaller one
-                from the larger one to produce an even number, remove
-                more factors of two, and repeat. */
-             while (u != v)
-               {
-                 if (u > v)
-                   {
-                     u -= v;
-                     while ((u & 1) == 0)
-                       u >>= 1;
-                   }
-                 else
-                   {
-                     v -= u;
-                     while ((v & 1) == 0)
-                       v >>= 1;
-                   }
-               }
-             result = u << k;
-           }
-          return (SCM_POSFIXABLE (result)
-                 ? SCM_I_MAKINUM (result)
-                 : scm_i_inum2big (result));
-        }
+      if (SCM_I_INUMP (y))
+        return scm_integer_gcd_ii (SCM_I_INUM (x), SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-        {
-          SCM_SWAP (x, y);
-          goto big_inum;
-        }
+        return scm_integer_gcd_zi (y, SCM_I_INUM (x));
       else if (SCM_REALP (y) && scm_is_integer (y))
         goto handle_inexacts;
       else
@@ -2852,30 +2799,9 @@ scm_gcd (SCM x, SCM y)
   else if (SCM_BIGP (x))
     {
       if (SCM_I_INUMP (y))
-        {
-          scm_t_bits result;
-          scm_t_inum yy;
-        big_inum:
-          yy = SCM_I_INUM (y);
-          if (yy == 0)
-            return scm_abs (x);
-          if (yy < 0)
-           yy = -yy;
-          result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), yy);
-          scm_remember_upto_here_1 (x);
-          return (SCM_POSFIXABLE (result) 
-                 ? SCM_I_MAKINUM (result)
-                 : scm_from_unsigned_integer (result));
-        }
+        return scm_integer_gcd_zi (x, SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-        {
-          SCM result = scm_i_mkbig ();
-          mpz_gcd (SCM_I_BIG_MPZ (result),
-                  SCM_I_BIG_MPZ (x),
-                  SCM_I_BIG_MPZ (y));
-          scm_remember_upto_here_2 (x, y);
-          return scm_i_normbig (result);
-        }
+        return scm_integer_gcd_zz (x, y);
       else if (SCM_REALP (y) && scm_is_integer (y))
         goto handle_inexacts;
       else