Implement centered-divide with new integer lib

* libguile/integers.c (scm_integer_centered_divide_ii)
(scm_integer_centered_divide_iz, scm_integer_centered_divide_zi)
(scm_integer_centered_divide_zz): New internal functions.
(integer_centered_divide_zz): New helper.
* libguile/integers.h: Declare internal functions.
* libguile/numbers.c (scm_centered_divide): Use the new functions.
(scm_i_bigint_centered_divide): Remove unused helper.

3 files changed, 151 insertions, 143 deletions

diff --git a/libguile/integers.c b/libguile/integers.c
index cfc71a8c8..102278e44 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -1231,3 +1231,136 @@ scm_integer_centered_remainder_zz (SCM x, SCM y)
 {
   return integer_centered_remainder_zz (scm_bignum (x), scm_bignum (y));
 }
+
+static void
+integer_centered_divide_zz (struct scm_bignum *x, struct scm_bignum *y,
+                            SCM *qp, SCM *rp)
+{
+  mpz_t q, r, min_r, zx, zy;
+  mpz_init (q);
+  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 rr 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_qr (q, r, zx, zy);
+      mpz_neg (min_r, min_r);
+      if (mpz_cmp (r, min_r) < 0)
+	{
+	  mpz_sub_ui (q, q, 1);
+	  mpz_add (r, r, zy);
+	}
+    }
+  else
+    {
+      mpz_fdiv_qr (q, r, zx, zy);
+      if (mpz_cmp (r, min_r) < 0)
+	{
+	  mpz_add_ui (q, q, 1);
+	  mpz_sub (r, r, zy);
+	}
+    }
+  scm_remember_upto_here_2 (x, y);
+  mpz_clear (min_r);
+  *qp = take_mpz (q);
+  *rp = take_mpz (r);
+}
+
+void
+scm_integer_centered_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
+{
+  if (y == 0)
+    scm_num_overflow ("centered-divide");
+
+  scm_t_inum q = x / y;
+  scm_t_inum r = x % y;
+  if (x > 0)
+    {
+      if (y > 0)
+        {
+          if (r >= (y + 1) / 2)
+            { q++; r -= y; }
+        }
+      else
+        {
+          if (r >= (1 - y) / 2)
+            { q--; r += y; }
+        }
+    }
+  else
+    {
+      if (y > 0)
+        {
+          if (r < -y / 2)
+            { q--; r += y; }
+        }
+      else
+        {
+          if (r < y / 2)
+            { q++; r -= y; }
+        }
+    }
+  *qp = long_to_scm (q);
+  *rp = SCM_I_MAKINUM (r);
+}
+
+void
+scm_integer_centered_divide_iz (scm_t_inum x, SCM y, SCM *qp, SCM *rp)
+{
+  integer_centered_divide_zz (long_to_bignum (x), scm_bignum (y), qp, rp);
+}
+
+void
+scm_integer_centered_divide_zi (SCM x, scm_t_inum y, SCM *qp, SCM *rp)
+{
+  if (y == 0)
+    scm_num_overflow ("centered-divide");
+
+  mpz_t q, zx;
+  mpz_init (q);
+  alias_bignum_to_mpz (scm_bignum (x), zx);
+  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_q_ui (q, zx, y);
+      if (r < -y / 2)
+        {
+          mpz_sub_ui (q, q, 1);
+          r += y;
+        }
+    }
+  else
+    {
+      r = - mpz_cdiv_q_ui (q, zx, -y);
+      mpz_neg (q, q);
+      if (r < y / 2)
+        {
+          mpz_add_ui (q, q, 1);
+          r -= y;
+        }
+    }
+  scm_remember_upto_here_1 (x);
+  *qp = take_mpz (q);
+  *rp = SCM_I_MAKINUM (r);
+}
+
+void
+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);
+}
+
diff --git a/libguile/integers.h b/libguile/integers.h
index 13c69e374..f941eefc6 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -96,6 +96,15 @@ 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);
 
+SCM_INTERNAL void scm_integer_centered_divide_ii (scm_t_inum x, scm_t_inum y,
+                                                  SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_centered_divide_iz (scm_t_inum x, SCM y,
+                                                  SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_centered_divide_zi (SCM x, scm_t_inum y,
+                                                  SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_centered_divide_zz (SCM x, SCM y,
+                                                  SCM *qp, SCM *rp);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 1c915b75a..284d17fd4 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -2320,7 +2320,6 @@ scm_i_exact_rational_centered_remainder (SCM x, SCM y)
 
 static void scm_i_inexact_centered_divide (double x, double y,
 					   SCM *qp, SCM *rp);
-static void scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp);
 static void scm_i_exact_rational_centered_divide (SCM x, SCM y,
 						  SCM *qp, SCM *rp);
 
@@ -2352,57 +2351,15 @@ SCM_PRIMITIVE_GENERIC (scm_i_centered_divide, "centered/", 2, 0, 0,
 void
 scm_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
 {
-  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_divide);
-	  else
-	    {
-	      scm_t_inum qq = xx / yy;
-	      scm_t_inum rr = xx % yy;
-	      if (SCM_LIKELY (xx > 0))
-		{
-		  if (SCM_LIKELY (yy > 0))
-		    {
-		      if (rr >= (yy + 1) / 2)
-			{ qq++; rr -= yy; }
-		    }
-		  else
-		    {
-		      if (rr >= (1 - yy) / 2)
-			{ qq--; rr += yy; }
-		    }
-		}
-	      else
-		{
-		  if (SCM_LIKELY (yy > 0))
-		    {
-		      if (rr < -yy / 2)
-			{ qq--; rr += yy; }
-		    }
-		  else
-		    {
-		      if (rr < yy / 2)
-			{ qq++; rr -= yy; }
-		    }
-		}
-	      if (SCM_LIKELY (SCM_FIXABLE (qq)))
-		*qp = SCM_I_MAKINUM (qq);
-	      else
-		*qp = scm_i_inum2big (qq);
-	      *rp = SCM_I_MAKINUM (rr);
-	    }
-	}
+      if (SCM_I_INUMP (y))
+        scm_integer_centered_divide_ii (SCM_I_INUM (x), SCM_I_INUM (y), qp, rp);
       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_divide */
-        scm_i_bigint_centered_divide (scm_i_long2big (xx), y, qp, rp);
+        scm_integer_centered_divide_iz (SCM_I_INUM (x), y, qp, rp);
       else if (SCM_REALP (y))
-	scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y), qp, rp);
+	scm_i_inexact_centered_divide (SCM_I_INUM (x), SCM_REAL_VALUE (y),
+                                       qp, rp);
       else if (SCM_FRACTIONP (y))
 	scm_i_exact_rational_centered_divide (x, y, qp, rp);
       else
@@ -2411,49 +2368,10 @@ scm_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
     }
   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_divide);
-	  else
-	    {
-	      SCM q = scm_i_mkbig ();
-	      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_q_ui (SCM_I_BIG_MPZ (q),
-					SCM_I_BIG_MPZ (x), yy);
-		  scm_remember_upto_here_1 (x);
-		  if (rr < -yy / 2)
-		    {
-		      mpz_sub_ui (SCM_I_BIG_MPZ (q),
-				  SCM_I_BIG_MPZ (q), 1);
-		      rr += yy;
-		    }
-		}
-	      else
-		{
-		  rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
-					SCM_I_BIG_MPZ (x), -yy);
-		  scm_remember_upto_here_1 (x);
-		  mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
-		  if (rr < yy / 2)
-		    {
-		      mpz_add_ui (SCM_I_BIG_MPZ (q),
-				  SCM_I_BIG_MPZ (q), 1);
-		      rr -= yy;
-		    }
-		}
-	      *qp = scm_i_normbig (q);
-	      *rp = SCM_I_MAKINUM (rr);
-	    }
-	}
+      if (SCM_I_INUMP (y))
+        scm_integer_centered_divide_zi (x, SCM_I_INUM (y), qp, rp);
       else if (SCM_BIGP (y))
-	scm_i_bigint_centered_divide (x, y, qp, rp);
+        scm_integer_centered_divide_zz (x, y, qp, rp);
       else if (SCM_REALP (y))
 	scm_i_inexact_centered_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y),
                                        qp, rp);
@@ -2507,58 +2425,6 @@ scm_i_inexact_centered_divide (double x, double y, SCM *qp, SCM *rp)
   *rp = scm_i_from_double (r);
 }
 
-/* Assumes that both x and y are bigints, though
-   x might be able to fit into a fixnum. */
-static void
-scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
-{
-  SCM q, 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 */
-  q = scm_i_mkbig ();
-  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_qr (SCM_I_BIG_MPZ (q), 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_sub_ui (SCM_I_BIG_MPZ (q),
-		      SCM_I_BIG_MPZ (q), 1);
-	  mpz_add (SCM_I_BIG_MPZ (r),
-		   SCM_I_BIG_MPZ (r),
-		   SCM_I_BIG_MPZ (y));
-	}
-    }
-  else
-    {
-      mpz_fdiv_qr (SCM_I_BIG_MPZ (q), 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_add_ui (SCM_I_BIG_MPZ (q),
-		      SCM_I_BIG_MPZ (q), 1);
-	  mpz_sub (SCM_I_BIG_MPZ (r),
-		   SCM_I_BIG_MPZ (r),
-		   SCM_I_BIG_MPZ (y));
-	}
-    }
-  scm_remember_upto_here_2 (x, y);
-  *qp = scm_i_normbig (q);
-  *rp = scm_i_normbig (r);
-}
-
 static void
 scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
 {