[Guile-commits] 17/69: Implement round-quotient with new integer lib

[Andy Wingo]

wingo pushed a commit to branch wip-inline-digits
in repository guile.

commit aa258402d35e0b5e2776b004e47d75467dca165b
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Mon Dec 13 11:42:17 2021 +0100

    Implement round-quotient with new integer lib
    
    * libguile/integers.c (scm_integer_round_quotient_ii)
    (scm_integer_round_quotient_iz, scm_integer_round_quotient_zi)
    (scm_integer_round_quotient_zz): New internal functions.
    (integer_round_quotient_zz): New helper.
    (long_sign, bignum_cmp_long): New helpers.
    * libguile/integers.h: Declare internal functions.
    * libguile/numbers.c (scm_round_quotient): Use the new functions.
    (scm_i_bigint_round_quotient): Remove unused helper.
---
 libguile/integers.c | 141 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 libguile/integers.h |   5 ++
 libguile/numbers.c  | 127 ++++------------------------------------------
 3 files changed, 154 insertions(+), 119 deletions(-)

diff --git a/libguile/integers.c b/libguile/integers.c
index 102278e44..ebe888d7a 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -416,7 +416,7 @@ scm_integer_floor_remainder_iz (scm_t_inum x, SCM y)
 SCM
 scm_integer_floor_remainder_zi (SCM x, scm_t_inum y)
 {
-  if (SCM_UNLIKELY (y == 0))
+  if (y == 0)
     scm_num_overflow ("floor-remainder");
   else
     {
@@ -1364,3 +1364,142 @@ scm_integer_centered_divide_zz (SCM x, SCM y, SCM *qp, 
SCM *rp)
   integer_centered_divide_zz (scm_bignum (x), scm_bignum (y), qp, rp);
 }
 
+static SCM
+integer_round_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
+{
+  mpz_t q, r, r2, zx, zy;
+  int cmp, needs_adjustment;
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+  mpz_init (q);
+  mpz_init (r);
+  mpz_init (r2);
+  alias_bignum_to_mpz (x, zx);
+  alias_bignum_to_mpz (y, zy);
+
+  mpz_fdiv_qr (q, r, zx, zy);
+  mpz_mul_2exp (r2, r, 1);  /* r2 = 2*r */
+  scm_remember_upto_here_1 (x);
+
+  cmp = mpz_cmpabs (r2, zy);
+  if (mpz_odd_p (q))
+    needs_adjustment = (cmp >= 0);
+  else
+    needs_adjustment = (cmp > 0);
+  scm_remember_upto_here_1 (y);
+
+  if (needs_adjustment)
+    mpz_add_ui (q, q, 1);
+
+  mpz_clear (r);
+  mpz_clear (r2);
+  return take_mpz (q);
+}
+
+SCM
+scm_integer_round_quotient_ii (scm_t_inum x, scm_t_inum y)
+{
+  if (y == 0)
+    scm_num_overflow ("round-quotient");
+
+  scm_t_inum q = x / y;
+  scm_t_inum r = x % y;
+  scm_t_inum ay = y;
+  scm_t_inum r2 = 2 * r;
+
+  if (y < 0)
+    {
+      ay = -ay;
+      r2 = -r2;
+    }
+
+  if (q & 1L)
+    {
+      if (r2 >= ay)
+        q++;
+      else if (r2 <= -ay)
+        q--;
+    }
+  else
+    {
+      if (r2 > ay)
+        q++;
+      else if (r2 < -ay)
+        q--;
+    }
+  return long_to_scm (q);
+}
+
+SCM
+scm_integer_round_quotient_iz (scm_t_inum x, SCM y)
+{
+  return integer_round_quotient_zz (long_to_bignum (x), scm_bignum (y));
+}
+
+SCM
+scm_integer_round_quotient_zi (SCM x, scm_t_inum y)
+{
+  if (y == 0)
+    scm_num_overflow ("round-quotient");
+  if (y == 1)
+    return x;
+
+  mpz_t q, zx;
+  mpz_init (q);
+  alias_bignum_to_mpz (scm_bignum (x), zx);
+  scm_t_inum r;
+  int needs_adjustment;
+
+  if (y > 0)
+    {
+      r = mpz_fdiv_q_ui (q, zx, y);
+      if (mpz_odd_p (q))
+        needs_adjustment = (2*r >= y);
+      else
+        needs_adjustment = (2*r > y);
+    }
+  else
+    {
+      r = - mpz_cdiv_q_ui (q, zx, -y);
+      mpz_neg (q, q);
+      if (mpz_odd_p (q))
+        needs_adjustment = (2*r <= y);
+      else
+        needs_adjustment = (2*r < y);
+    }
+  scm_remember_upto_here_1 (x);
+  if (needs_adjustment)
+    mpz_add_ui (q, q, 1);
+  return take_mpz (q);
+}
+
+SCM
+scm_integer_round_quotient_zz (SCM x, SCM y)
+{
+  SCM q, r, r2;
+  int cmp, needs_adjustment;
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+  q = scm_i_mkbig ();
+  r = scm_i_mkbig ();
+  r2 = scm_i_mkbig ();
+
+  mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+              SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+  mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1);  /* r2 = 2*r */
+  scm_remember_upto_here_2 (x, r);
+
+  cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y));
+  if (mpz_odd_p (SCM_I_BIG_MPZ (q)))
+    needs_adjustment = (cmp >= 0);
+  else
+    needs_adjustment = (cmp > 0);
+  scm_remember_upto_here_2 (r2, y);
+
+  if (needs_adjustment)
+    mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
+
+  return scm_i_normbig (q);
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index f941eefc6..dcbc32a10 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -105,6 +105,11 @@ SCM_INTERNAL void scm_integer_centered_divide_zi (SCM x, 
scm_t_inum y,
 SCM_INTERNAL void scm_integer_centered_divide_zz (SCM x, SCM y,
                                                   SCM *qp, SCM *rp);
 
+SCM_INTERNAL SCM scm_integer_round_quotient_ii (scm_t_inum x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_round_quotient_iz (scm_t_inum x, SCM y);
+SCM_INTERNAL SCM scm_integer_round_quotient_zi (SCM x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_round_quotient_zz (SCM x, SCM y);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 284d17fd4..d3b794578 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -2439,7 +2439,6 @@ scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM 
*qp, SCM *rp)
 }
 
 static SCM scm_i_inexact_round_quotient (double x, double y);
-static SCM scm_i_bigint_round_quotient (SCM x, SCM y);
 static SCM scm_i_exact_rational_round_quotient (SCM x, SCM y);
 
 SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0,
@@ -2459,55 +2458,15 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, 
"round-quotient", 2, 0, 0,
                       "@end lisp")
 #define FUNC_NAME s_scm_round_quotient
 {
-  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_round_quotient);
-         else
-           {
-             scm_t_inum qq = xx / yy;
-             scm_t_inum rr = xx % yy;
-             scm_t_inum ay = yy;
-             scm_t_inum r2 = 2 * rr;
-
-             if (SCM_LIKELY (yy < 0))
-               {
-                 ay = -ay;
-                 r2 = -r2;
-               }
-
-             if (qq & 1L)
-               {
-                 if (r2 >= ay)
-                   qq++;
-                 else if (r2 <= -ay)
-                   qq--;
-               }
-             else
-               {
-                 if (r2 > ay)
-                   qq++;
-                 else if (r2 < -ay)
-                   qq--;
-               }
-             if (SCM_LIKELY (SCM_FIXABLE (qq)))
-               return SCM_I_MAKINUM (qq);
-             else
-               return scm_i_inum2big (qq);
-           }
-       }
+      if (SCM_I_INUMP (y))
+        return scm_integer_round_quotient_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_round_quotient */
-         return scm_i_bigint_round_quotient (scm_i_long2big (xx), y);
-       }
+        return scm_integer_round_quotient_iz (SCM_I_INUM (x), y);
       else if (SCM_REALP (y))
-       return scm_i_inexact_round_quotient (xx, SCM_REAL_VALUE (y));
+       return scm_i_inexact_round_quotient (SCM_I_INUM (x),
+                                             SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_exact_rational_round_quotient (x, y);
       else
@@ -2516,46 +2475,10 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, 
"round-quotient", 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_round_quotient);
-         else if (SCM_UNLIKELY (yy == 1))
-           return x;
-         else
-           {
-             SCM q = scm_i_mkbig ();
-             scm_t_inum rr;
-             int needs_adjustment;
-
-             if (yy > 0)
-               {
-                 rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q),
-                                     SCM_I_BIG_MPZ (x), yy);
-                 if (mpz_odd_p (SCM_I_BIG_MPZ (q)))
-                   needs_adjustment = (2*rr >= yy);
-                 else
-                   needs_adjustment = (2*rr > yy);
-               }
-             else
-               {
-                 rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
-                                       SCM_I_BIG_MPZ (x), -yy);
-                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
-                 if (mpz_odd_p (SCM_I_BIG_MPZ (q)))
-                   needs_adjustment = (2*rr <= yy);
-                 else
-                   needs_adjustment = (2*rr < yy);
-               }
-             scm_remember_upto_here_1 (x);
-             if (needs_adjustment)
-               mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
-             return scm_i_normbig (q);
-           }
-       }
+      if (SCM_I_INUMP (y))
+        return scm_integer_round_quotient_zi (x, SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-       return scm_i_bigint_round_quotient (x, y);
+        return scm_integer_round_quotient_zz (x, y);
       else if (SCM_REALP (y))
        return scm_i_inexact_round_quotient
          (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
@@ -2601,38 +2524,6 @@ scm_i_inexact_round_quotient (double x, double y)
     return scm_i_from_double (scm_c_round (x / y));
 }
 
-/* Assumes that both x and y are bigints, though
-   x might be able to fit into a fixnum. */
-static SCM
-scm_i_bigint_round_quotient (SCM x, SCM y)
-{
-  SCM q, r, r2;
-  int cmp, needs_adjustment;
-
-  /* Note that x might be small enough to fit into a
-     fixnum, so we must not let it escape into the wild */
-  q = scm_i_mkbig ();
-  r = scm_i_mkbig ();
-  r2 = scm_i_mkbig ();
-
-  mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
-              SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
-  mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1);  /* r2 = 2*r */
-  scm_remember_upto_here_2 (x, r);
-
-  cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y));
-  if (mpz_odd_p (SCM_I_BIG_MPZ (q)))
-    needs_adjustment = (cmp >= 0);
-  else
-    needs_adjustment = (cmp > 0);
-  scm_remember_upto_here_2 (r2, y);
-
-  if (needs_adjustment)
-    mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
-
-  return scm_i_normbig (q);
-}
-
 static SCM
 scm_i_exact_rational_round_quotient (SCM x, SCM y)
 {