[Guile-commits] 14/85: Implement centered-remainder with new integer lib

[Andy Wingo]

wingo pushed a commit to branch main
in repository guile.

commit 99b046d58bd02a916703c9d1f75bd3d5c6cf9ebf
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Mon Dec 13 10:51:15 2021 +0100

    Implement centered-remainder with new integer lib
    
    * libguile/integers.c (scm_integer_centered_remainder_ii)
    (scm_integer_centered_remainder_iz, scm_integer_centered_remainder_zi)
    (scm_integer_centered_remainder_zz): New internal functions.
    * libguile/integers.h: Declare internal functions.
    * libguile/numbers.c (scm_centered_remainder): Use the new functions.
    (scm_i_bigint_centered_remainder): Remove helper.
---
 libguile/integers.c | 112 +++++++++++++++++++++++++++++++++++++++++++++
 libguile/integers.h |   5 +++
 libguile/numbers.c  | 127 +++++-----------------------------------------------
 3 files changed, 127 insertions(+), 117 deletions(-)

diff --git a/libguile/integers.c b/libguile/integers.c
index e9c3a067c..cfc71a8c8 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -1119,3 +1119,115 @@ scm_integer_centered_quotient_zz (SCM x, SCM y)
 {
   return integer_centered_quotient_zz (scm_bignum (x), scm_bignum (y));
 }
+
+static SCM
+integer_centered_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
+{
+  mpz_t r, min_r, zx, zy;
+  mpz_init (r);
+  mpz_init (min_r);
+  alias_bignum_to_mpz (x, zx);
+  alias_bignum_to_mpz (y, zy);
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+
+  /* min_r will eventually become -abs(y)/2 */
+  mpz_tdiv_q_2exp (min_r, zy, 1);
+
+  /* Arrange for r to initially be non-positive, because that simplifies
+     the test to see if it is within the needed bounds. */
+  if (mpz_sgn (zy) > 0)
+    {
+      mpz_cdiv_r (r, zx, zy);
+      mpz_neg (min_r, min_r);
+      if (mpz_cmp (r, min_r) < 0)
+       mpz_add (r, r, zy);
+    }
+  else
+    {
+      mpz_fdiv_r (r, zx, zy);
+      if (mpz_cmp (r, min_r) < 0)
+       mpz_sub (r, r, zy);
+    }
+  scm_remember_upto_here_2 (x, y);
+  mpz_clear (min_r);
+  return take_mpz (r);
+}
+
+SCM
+scm_integer_centered_remainder_ii (scm_t_inum x, scm_t_inum y)
+{
+  if (y == 0)
+    scm_num_overflow ("centered-remainder");
+
+  scm_t_inum r = x % y;
+  if (x > 0)
+    {
+      if (y > 0)
+        {
+          if (r >= (y + 1) / 2)
+            r -= y;
+        }
+      else
+        {
+          if (r >= (1 - y) / 2)
+            r += y;
+        }
+    }
+  else
+    {
+      if (y > 0)
+        {
+          if (r < -y / 2)
+            r += y;
+        }
+      else
+        {
+          if (r < y / 2)
+            r -= y;
+        }
+    }
+  return SCM_I_MAKINUM (r);
+}
+
+SCM
+scm_integer_centered_remainder_iz (scm_t_inum x, SCM y)
+{
+  return integer_centered_remainder_zz (long_to_bignum (x),
+                                        scm_bignum (y));
+}
+
+SCM
+scm_integer_centered_remainder_zi (SCM x, scm_t_inum y)
+{
+  mpz_t zx;
+  alias_bignum_to_mpz (scm_bignum (x), zx);
+
+  if (y == 0)
+    scm_num_overflow ("centered-remainder");
+
+  scm_t_inum r;
+  /* Arrange for r to initially be non-positive, because that simplifies
+     the test to see if it is within the needed bounds. */
+  if (y > 0)
+    {
+      r = - mpz_cdiv_ui (zx, y);
+      if (r < -y / 2)
+        r += y;
+    }
+  else
+    {
+      r = - mpz_cdiv_ui (zx, -y);
+      if (r < y / 2)
+        r -= y;
+    }
+  scm_remember_upto_here_1 (x);
+  return SCM_I_MAKINUM (r);
+}
+
+SCM
+scm_integer_centered_remainder_zz (SCM x, SCM y)
+{
+  return integer_centered_remainder_zz (scm_bignum (x), scm_bignum (y));
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index 44dc0e376..13c69e374 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -91,6 +91,11 @@ SCM_INTERNAL SCM scm_integer_centered_quotient_iz 
(scm_t_inum x, SCM y);
 SCM_INTERNAL SCM scm_integer_centered_quotient_zi (SCM x, scm_t_inum y);
 SCM_INTERNAL SCM scm_integer_centered_quotient_zz (SCM x, SCM y);
 
+SCM_INTERNAL SCM scm_integer_centered_remainder_ii (scm_t_inum x, scm_t_inum 
y);
+SCM_INTERNAL SCM scm_integer_centered_remainder_iz (scm_t_inum x, SCM y);
+SCM_INTERNAL SCM scm_integer_centered_remainder_zi (SCM x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_centered_remainder_zz (SCM x, SCM y);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index fd3a8c6c7..1c915b75a 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -2208,7 +2208,6 @@ scm_i_exact_rational_centered_quotient (SCM x, SCM y)
 }
 
 static SCM scm_i_inexact_centered_remainder (double x, double y);
-static SCM scm_i_bigint_centered_remainder (SCM x, SCM y);
 static SCM scm_i_exact_rational_centered_remainder (SCM x, SCM y);
 
 SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
@@ -2227,54 +2226,16 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, 
"centered-remainder", 2, 0, 0,
                       "@end lisp")
 #define FUNC_NAME s_scm_centered_remainder
 {
-  if (SCM_LIKELY (SCM_I_INUMP (x)))
+  if (SCM_I_INUMP (x))
     {
-      scm_t_inum xx = SCM_I_INUM (x);
-      if (SCM_LIKELY (SCM_I_INUMP (y)))
-       {
-         scm_t_inum yy = SCM_I_INUM (y);
-         if (SCM_UNLIKELY (yy == 0))
-           scm_num_overflow (s_scm_centered_remainder);
-         else
-           {
-             scm_t_inum rr = xx % yy;
-             if (SCM_LIKELY (xx > 0))
-               {
-                 if (SCM_LIKELY (yy > 0))
-                   {
-                     if (rr >= (yy + 1) / 2)
-                       rr -= yy;
-                   }
-                 else
-                   {
-                     if (rr >= (1 - yy) / 2)
-                       rr += yy;
-                   }
-               }
-             else
-               {
-                 if (SCM_LIKELY (yy > 0))
-                   {
-                     if (rr < -yy / 2)
-                       rr += yy;
-                   }
-                 else
-                   {
-                     if (rr < yy / 2)
-                       rr -= yy;
-                   }
-               }
-             return SCM_I_MAKINUM (rr);
-           }
-       }
+      if (SCM_I_INUMP (y))
+        return scm_integer_centered_remainder_ii (SCM_I_INUM (x),
+                                                  SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-       {
-         /* Pass a denormalized bignum version of x (even though it
-            can fit in a fixnum) to scm_i_bigint_centered_remainder */
-         return scm_i_bigint_centered_remainder (scm_i_long2big (xx), y);
-       }
+        return scm_integer_centered_remainder_iz (SCM_I_INUM (x), y);
       else if (SCM_REALP (y))
-       return scm_i_inexact_centered_remainder (xx, SCM_REAL_VALUE (y));
+       return scm_i_inexact_centered_remainder (SCM_I_INUM (x),
+                                                 SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_exact_rational_centered_remainder (x, y);
       else
@@ -2283,36 +2244,10 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, 
"centered-remainder", 2, 0, 0,
     }
   else if (SCM_BIGP (x))
     {
-      if (SCM_LIKELY (SCM_I_INUMP (y)))
-       {
-         scm_t_inum yy = SCM_I_INUM (y);
-         if (SCM_UNLIKELY (yy == 0))
-           scm_num_overflow (s_scm_centered_remainder);
-         else
-           {
-             scm_t_inum rr;
-             /* Arrange for rr to initially be non-positive,
-                because that simplifies the test to see
-                if it is within the needed bounds. */
-             if (yy > 0)
-               {
-                 rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy);
-                 scm_remember_upto_here_1 (x);
-                 if (rr < -yy / 2)
-                   rr += yy;
-               }
-             else
-               {
-                 rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy);
-                 scm_remember_upto_here_1 (x);
-                 if (rr < yy / 2)
-                   rr -= yy;
-               }
-             return SCM_I_MAKINUM (rr);
-           }
-       }
+      if (SCM_I_INUMP (y))
+        return scm_integer_centered_remainder_zi (x, SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-       return scm_i_bigint_centered_remainder (x, y);
+        return scm_integer_centered_remainder_zz (x, y);
       else if (SCM_REALP (y))
        return scm_i_inexact_centered_remainder
          (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
@@ -2372,48 +2307,6 @@ scm_i_inexact_centered_remainder (double x, double y)
   return scm_i_from_double (x - q * y);
 }
 
-/* Assumes that both x and y are bigints, though
-   x might be able to fit into a fixnum. */
-static SCM
-scm_i_bigint_centered_remainder (SCM x, SCM y)
-{
-  SCM r, min_r;
-
-  /* Note that x might be small enough to fit into a
-     fixnum, so we must not let it escape into the wild */
-  r = scm_i_mkbig ();
-
-  /* min_r will eventually become -abs(y)/2 */
-  min_r = scm_i_mkbig ();
-  mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
-                  SCM_I_BIG_MPZ (y), 1);
-
-  /* Arrange for rr to initially be non-positive,
-     because that simplifies the test to see
-     if it is within the needed bounds. */
-  if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
-    {
-      mpz_cdiv_r (SCM_I_BIG_MPZ (r),
-                 SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
-      mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
-      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
-       mpz_add (SCM_I_BIG_MPZ (r),
-                SCM_I_BIG_MPZ (r),
-                SCM_I_BIG_MPZ (y));
-    }
-  else
-    {
-      mpz_fdiv_r (SCM_I_BIG_MPZ (r),
-                 SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
-      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
-       mpz_sub (SCM_I_BIG_MPZ (r),
-                SCM_I_BIG_MPZ (r),
-                SCM_I_BIG_MPZ (y));
-    }
-  scm_remember_upto_here_2 (x, y);
-  return scm_i_normbig (r);
-}
-
 static SCM
 scm_i_exact_rational_centered_remainder (SCM x, SCM y)
 {