Add two new sets of fast quotient and remainder operators

* libguile/numbers.c (scm_euclidean_quo_and_rem, scm_euclidean_quotient,
  scm_euclidean_remainder, scm_centered_quo_and_rem,
  scm_centered_quotient, scm_centered_remainder): New extensible
  procedures `euclidean/', `euclidean-quotient', `euclidean-remainder',
  `centered/', `centered-quotient', `centered-remainder'.

* libguile/numbers.h: Add function prototypes.

* module/rnrs/base.scm: Remove incorrect stub implementations of `div',
  `mod', `div-and-mod', `div0', `mod0', and `div0-and-mod0'.  Instead do
  renaming imports of `euclidean-quotient', `euclidean-remainder',
  `euclidean/', `centered-quotient', `centered-remainder', and
  `centered/', which are equivalent to the R6RS operators.

* module/rnrs/arithmetic/fixnums.scm (fxdiv, fxmod, fxdiv-and-mod,
  fxdiv0, fxmod0, fxdiv0-and-mod0): Remove redundant checks for division
  by zero and unnecessary complexity.
  (fx+/carry): Remove unneeded calls to `inexact->exact'.

* module/rnrs/arithmetic/flonums.scm (fldiv, flmod, fldiv-and-mod,
  fldiv0, flmod0, fldiv0-and-mod0): Remove redundant checks for division
  by zero and unnecessary complexity.  Remove unneeded calls to
  `inexact->exact' and `exact->inexact'

* test-suite/tests/numbers.test: (test-eqv?): New internal predicate for
  comparing numerical outputs with expected values.

  Add extensive test code for `euclidean/', `euclidean-quotient',
  `euclidean-remainder', `centered/', `centered-quotient',
  `centered-remainder'.

* test-suite/tests/r6rs-arithmetic-fixnums.test: Fix some broken test
  cases, and remove `unresolved' test markers for `fxdiv', `fxmod',
  `fxdiv-and-mod', `fxdiv0', `fxmod0', and `fxdiv0-and-mod0'.

* test-suite/tests/r6rs-arithmetic-flonums.test: Remove `unresolved'
  test markers for `fldiv', `flmod', `fldiv-and-mod', `fldiv0',
  `flmod0', and `fldiv0-and-mod0'.

* doc/ref/api-data.texi (Arithmetic): Document `euclidean/',
  `euclidean-quotient', `euclidean-remainder', `centered/',
  `centered-quotient', and `centered-remainder'.

  (Operations on Integer Values): Add cross-references to `euclidean/'
  et al, from `quotient', `remainder', and `modulo'.

* doc/ref/r6rs.texi (rnrs base): Improve documentation for `div', `mod',
  `div-and-mod', `div0', `mod0', and `div0-and-mod0'.  Add
  cross-references to `euclidean/' et al.

* NEWS: Add NEWS entry.

diff --git a/NEWS b/NEWS
index f45795e..33e49a4 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -12,6 +12,29 @@ Changes in 1.9.15 (since the 1.9.14 prerelease):
 
 ** Changes and bugfixes in numerics code
 
+*** Added two new sets of fast quotient and remainder operators
+
+Added two new sets of fast quotient and remainder operator pairs with
+different semantics than the R5RS operators.  They support not only
+integers, but all reals, including exact rationals and inexact
+floating point numbers.
+
+These procedures accept two real numbers N and D, where the divisor D
+must be non-zero.  `euclidean-quotient' returns the integer Q and
+`euclidean-remainder' returns the real R such that N = Q*D + R and
+0 <= R < |D|.  `euclidean/' returns both Q and R, and is more
+efficient than computing each separately.  Note that when D > 0,
+`euclidean-quotient' returns floor(N/D), and when D < 0 it returns
+ceiling(N/D).
+
+`centered-quotient', `centered-remainder', and `centered/' are similar
+except that the range of remainders is -abs(D/2) <= R < abs(D/2), and
+`centered-quotient' rounds N/D to the nearest integer.
+
+Note that these operators are equivalent to the R6RS integer division
+operators `div', `mod', `div-and-mod', `div0', `mod0', and
+`div0-and-mod0'.
+
 *** `eqv?' and `equal?' now compare numbers equivalently
 
 scm_equal_p `equal?' now behaves equivalently to scm_eqv_p `eqv?' for
@@ -64,6 +87,12 @@ NaNs are neither finite nor infinite.
 
 *** R6RS base library changes
 
+**** `div', `mod', `div-and-mod', `div0', `mod0', `div0-and-mod0'
+
+Efficient versions of these R6RS division operators are now supported.
+See the NEWS entry entitled `Added two new sets of fast quotient and
+remainder operators' for more information.
+
 **** `infinite?' changes
 
 `infinite?' now returns #t for non-real complex infinities, and throws

diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 4256e18..b090782 100755 (executable)
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -897,6 +897,9 @@ sign as @var{n}.  In all cases quotient and remainder satisfy
 (remainder 13 4) @result{} 1
 (remainder -13 4) @result{} -1
 @end lisp
+
+See also @code{euclidean-quotient}, @code{euclidean-remainder} and
+related operations in @ref{Arithmetic}.
 @end deffn
 
 @c begin (texi-doc-string "guile" "modulo")
@@ -911,6 +914,9 @@ sign as @var{d}.
 (modulo 13 -4) @result{} -3
 (modulo -13 -4) @result{} -1
 @end lisp
+
+See also @code{euclidean-quotient}, @code{euclidean-remainder} and
+related operations in @ref{Arithmetic}.
 @end deffn
 
 @c begin (texi-doc-string "guile" "gcd")
@@ -1130,6 +1136,12 @@ Returns the magnitude or angle of @var{z} as a @code{double}.
 @rnindex ceiling
 @rnindex truncate
 @rnindex round
+@rnindex euclidean/
+@rnindex euclidean-quotient
+@rnindex euclidean-remainder
+@rnindex centered/
+@rnindex centered-quotient
+@rnindex centered-remainder
 
 The C arithmetic functions below always takes two arguments, while the
 Scheme functions can take an arbitrary number.  When you need to
@@ -1229,6 +1241,73 @@ respectively, but these functions take and return @code{double}
 values.
 @end deftypefn
 
+@deffn {Scheme Procedure} euclidean/ x y
+@deffnx {Scheme Procedure} euclidean-quotient x y
+@deffnx {Scheme Procedure} euclidean-remainder x y
+@deffnx {C Function} scm_euclidean_quo_and_rem (x y)
+@deffnx {C Function} scm_euclidean_quotient (x y)
+@deffnx {C Function} scm_euclidean_remainder (x y)
+These procedures accept two real numbers @var{x} and @var{y}, where the
+divisor @var{y} must be non-zero.  @code{euclidean-quotient} returns the
+integer @var{q} and @code{euclidean-remainder} returns the real number
+@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and
+@math{0 <= @var{r} < abs(@var{y})}.  @code{euclidean/} returns both @var{q} and
+@var{r}, and is more efficient than computing each separately.  Note
+that when @math{@var{y} > 0}, @code{euclidean-quotient} returns
+@math{floor(@var{x}/@var{y})}, otherwise it returns
+@math{ceiling(@var{x}/@var{y})}.
+
+Note that these operators are equivalent to the R6RS operators
+@code{div}, @code{mod}, and @code{div-and-mod}.
+
+@lisp
+(euclidean-quotient 123 10) @result{} 12
+(euclidean-remainder 123 10) @result{} 3
+(euclidean/ 123 10) @result{} 12 and 3
+(euclidean/ 123 -10) @result{} -12 and 3
+(euclidean/ -123 10) @result{} -13 and 7
+(euclidean/ -123 -10) @result{} 13 and 7
+(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8
+(euclidean/ 16/3 -10/7) @result{} -3 and 22/21
+@end lisp
+@end deffn
+
+@deffn {Scheme Procedure} centered/ x y
+@deffnx {Scheme Procedure} centered-quotient x y
+@deffnx {Scheme Procedure} centered-remainder x y
+@deffnx {C Function} scm_centered_quo_and_rem (x y)
+@deffnx {C Function} scm_centered_quotient (x y)
+@deffnx {C Function} scm_centered_remainder (x y)
+These procedures accept two real numbers @var{x} and @var{y}, where the
+divisor @var{y} must be non-zero.  @code{centered-quotient} returns the
+integer @var{q} and @code{centered-remainder} returns the real number
+@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and
+@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.  @code{centered/}
+returns both @var{q} and @var{r}, and is more efficient than computing
+each separately.
+
+Note that @code{centered-quotient} returns @math{@var{x}/@var{y}}
+rounded to the nearest integer.  When @math{@var{x}/@var{y}} lies
+exactly half-way between two integers, the tie is broken according to
+the sign of @var{y}.  If @math{@var{y} > 0}, ties are rounded toward
+positive infinity, otherwise they are rounded toward negative infinity.
+This is a consequence of the requirement that @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.
+
+Note that these operators are equivalent to the R6RS operators
+@code{div0}, @code{mod0}, and @code{div0-and-mod0}.
+
+@lisp
+(centered-quotient 123 10) @result{} 12
+(centered-remainder 123 10) @result{} 3
+(centered/ 123 10) @result{} 12 and 3
+(centered/ 123 -10) @result{} -12 and 3
+(centered/ -123 10) @result{} -12 and -3
+(centered/ -123 -10) @result{} 12 and -3
+(centered/ -123.2 -63.5) @result{} 2.0 and 3.8
+(centered/ 16/3 -10/7) @result{} -4 and -8/21
+@end lisp
+@end deffn
+
 @node Scientific
 @subsubsection Scientific Functions
 

diff --git a/doc/ref/r6rs.texi b/doc/ref/r6rs.texi
index 5fee65f..d6d9d9f 100644 (file)
--- a/doc/ref/r6rs.texi
+++ b/doc/ref/r6rs.texi
@@ -1,6 +1,6 @@
 @c -*-texinfo-*-
 @c This is part of the GNU Guile Reference Manual.
-@c Copyright (C)  2010
+@c Copyright (C)  2010, 2011
 @c   Free Software Foundation, Inc.
 @c See the file guile.texi for copying conditions.
 
@@ -464,21 +464,65 @@ grouped below by the existing manual sections to which they correspond.
 @xref{Arithmetic}, for documentation.
 @end deffn
 
-@deffn {Scheme Procedure} div x1 x2
-@deffnx {Scheme Procedure} mod x1 x2
-@deffnx {Scheme Procedure} div-and-mod x1 x2
-These procedures implement number-theoretic division.
+@rnindex div
+@rnindex mod
+@rnindex div-and-mod
+@deffn {Scheme Procedure} div x y
+@deffnx {Scheme Procedure} mod x y
+@deffnx {Scheme Procedure} div-and-mod x y
+These procedures accept two real numbers @var{x} and @var{y}, where the
+divisor @var{y} must be non-zero.  @code{div} returns the integer @var{q}
+and @code{mod} returns the real number @var{r} such that
+@math{@var{x} = @var{q}*@var{y} + @var{r}} and @math{0 <= @var{r} < abs(@var{y})}.
+@code{div-and-mod} returns both @var{q} and @var{r}, and is more
+efficient than computing each separately.  Note that when @math{@var{y} > 0},
+@code{div} returns @math{floor(@var{x}/@var{y})}, otherwise
+it returns @math{ceiling(@var{x}/@var{y})}.
 
-@code{div-and-mod} returns two values, the respective results of
-@code{(div x1 x2)} and @code{(mod x1 x2)}.
+@lisp
+(div 123 10) @result{} 12
+(mod 123 10) @result{} 3
+(div-and-mod 123 10) @result{} 12 and 3
+(div-and-mod 123 -10) @result{} -12 and 3
+(div-and-mod -123 10) @result{} -13 and 7
+(div-and-mod -123 -10) @result{} 13 and 7
+(div-and-mod -123.2 -63.5) @result{} 2.0 and 3.8
+(div-and-mod 16/3 -10/7) @result{} -3 and 22/21
+@end lisp
 @end deffn
 
-@deffn {Scheme Procedure} div0 x1 x2
-@deffnx {Scheme Procedure} mod0 x1 x2
-@deffnx {Scheme Procedure} div0-and-mod0 x1 x2
-These procedures are similar to @code{div}, @code{mod}, and 
-@code{div-and-mod}, except that @code{mod0} returns values that lie
-within a half-open interval centered on zero.
+@rnindex div0
+@rnindex mod0
+@rnindex div0-and-mod0
+@deffn {Scheme Procedure} div0 x y
+@deffnx {Scheme Procedure} mod0 x y
+@deffnx {Scheme Procedure} div0-and-mod0 x y
+These procedures accept two real numbers @var{x} and @var{y}, where the
+divisor @var{y} must be non-zero.  @code{div0} returns the
+integer @var{q} and @code{rem0} returns the real number
+@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and
+@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.  @code{div0-and-mod0}
+returns both @var{q} and @var{r}, and is more efficient than computing
+each separately.
+
+Note that @code{div0} returns @math{@var{x}/@var{y}} rounded to the
+nearest integer.  When @math{@var{x}/@var{y}} lies exactly half-way
+between two integers, the tie is broken according to the sign of
+@var{y}.  If @math{@var{y} > 0}, ties are rounded toward positive
+infinity, otherwise they are rounded toward negative infinity.
+This is a consequence of the requirement that
+@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.
+
+@lisp
+(div0 123 10) @result{} 12
+(mod0 123 10) @result{} 3
+(div0-and-mod0 123 10) @result{} 12 and 3
+(div0-and-mod0 123 -10) @result{} -12 and 3
+(div0-and-mod0 -123 10) @result{} -12 and -3
+(div0-and-mod0 -123 -10) @result{} 12 and -3
+(div0-and-mod0 -123.2 -63.5) @result{} 2.0 and 3.8
+(div0-and-mod0 16/3 -10/7) @result{} -4 and -8/21
+@end lisp
 @end deffn
 
 @deffn {Scheme Procedure} exact-integer-sqrt k

diff --git a/libguile/numbers.c b/libguile/numbers.c
index e15a6ac..4515dc9 100644 (file)
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -105,6 +105,7 @@ typedef scm_t_signed_bits scm_t_inum;
 
 
 static SCM flo0;
+static SCM exactly_one_half;
 
 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
 
@@ -1054,6 +1055,1230 @@ scm_modulo (SCM x, SCM y)
     SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo);
 }
 
+static SCM scm_i_inexact_euclidean_quotient (double x, double y);
+static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} such that\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "where @math{0 <= @var{r} < abs(@var{y})}.\n"
+                      "@lisp\n"
+                      "(euclidean-quotient 123 10) @result{} 12\n"
+                      "(euclidean-quotient 123 -10) @result{} -12\n"
+                      "(euclidean-quotient -123 10) @result{} -13\n"
+                      "(euclidean-quotient -123 -10) @result{} 13\n"
+                      "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
+                      "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_quotient
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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_euclidean_quotient);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             if (xx < qq * yy)
+               {
+                 if (yy > 0)
+                   qq--;
+                 else
+                   qq++;
+               }
+             return SCM_I_MAKINUM (qq);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         if (SCM_I_INUM (x) >= 0)
+           return SCM_INUM0;
+         else
+           return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y)));
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  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_euclidean_quotient);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             if (yy > 0)
+               mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy);
+             else
+               {
+                 mpz_fdiv_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));
+               }
+             scm_remember_upto_here_1 (x);
+             return scm_i_normbig (q);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM q = scm_i_mkbig ();
+         if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           mpz_fdiv_q (SCM_I_BIG_MPZ (q),
+                       SCM_I_BIG_MPZ (x),
+                       SCM_I_BIG_MPZ (y));
+         else
+           mpz_cdiv_q (SCM_I_BIG_MPZ (q),
+                       SCM_I_BIG_MPZ (x),
+                       SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_i_normbig (q);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_quotient
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
+                       s_scm_euclidean_quotient);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_quotient (double x, double y)
+{
+  if (SCM_LIKELY (y > 0))
+    return scm_from_double (floor (x / y));
+  else if (SCM_LIKELY (y < 0))
+    return scm_from_double (ceil (x / y));
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_quotient);  /* or return a NaN? */
+  else
+    return scm_nan ();
+}
+
+/* Compute exact euclidean_quotient the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_quotient (SCM x, SCM y)
+{
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
+                       s_scm_euclidean_quotient);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                       s_scm_euclidean_quotient);
+  else if (scm_is_true (scm_positive_p (y)))
+    return scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    return scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_quotient);
+}
+
+static SCM scm_i_inexact_euclidean_remainder (double x, double y);
+static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the real number @var{r} such that\n"
+                      "@math{0 <= @var{r} < abs(@var{y})} and\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "for some integer @var{q}.\n"
+                      "@lisp\n"
+                      "(euclidean-remainder 123 10) @result{} 3\n"
+                      "(euclidean-remainder 123 -10) @result{} 3\n"
+                      "(euclidean-remainder -123 10) @result{} 7\n"
+                      "(euclidean-remainder -123 -10) @result{} 7\n"
+                      "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
+                      "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_remainder
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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_euclidean_remainder);
+         else
+           {
+             scm_t_inum rr = SCM_I_INUM (x) % yy;
+             if (rr >= 0)
+               return SCM_I_MAKINUM (rr);
+             else if (yy > 0)
+               return SCM_I_MAKINUM (rr + yy);
+             else
+               return SCM_I_MAKINUM (rr - yy);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         scm_t_inum xx = SCM_I_INUM (x);
+         if (xx >= 0)
+           return x;
+         else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             return scm_i_normbig (r);
+           }
+         else
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r));
+             return scm_i_normbig (r);
+           }
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  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_euclidean_remainder);
+         else
+           {
+             scm_t_inum rr;
+             if (yy < 0)
+               yy = -yy;
+             rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy);
+             scm_remember_upto_here_1 (x);
+             return SCM_I_MAKINUM (rr);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM r = scm_i_mkbig ();
+         mpz_mod (SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x),
+                  SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_i_normbig (r);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_remainder
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
+                       s_scm_euclidean_remainder);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_remainder (double x, double y)
+{
+  double q;
+
+  /* Although it would be more efficient to use fmod here, we can't
+     because it would in some cases produce results inconsistent with
+     scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
+     even close).  In particular, when x is very close to a multiple of
+     y, then r might be either 0.0 or abs(y)-epsilon, but those two
+     cases must correspond to different choices of q.  If r = 0.0 then q
+     must be x/y, and if r = abs(y) then q must be (x-r)/y.  If quotient
+     chooses one and remainder chooses the other, it would be bad.  This
+     problem was observed with x = 130.0 and y = 10/7. */
+  if (SCM_LIKELY (y > 0))
+    q = floor (x / y);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x / y);
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_remainder);  /* or return a NaN? */
+  else
+    return scm_nan ();
+  return scm_from_double (x - q * y);
+}
+
+/* Compute exact euclidean_remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_remainder (SCM x, SCM y)
+{
+  SCM q;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
+                       s_scm_euclidean_remainder);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                       s_scm_euclidean_remainder);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_remainder);
+  return scm_difference (x, scm_product (y, q));
+}
+
+
+static SCM scm_i_inexact_euclidean_quo_and_rem (double x, double y);
+static SCM scm_i_slow_exact_euclidean_quo_and_rem (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} and the real number @var{r}\n"
+                      "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "and @math{0 <= @var{r} < abs(@var{y})}.\n"
+                      "@lisp\n"
+                      "(euclidean/ 123 10) @result{} 12 and 3\n"
+                      "(euclidean/ 123 -10) @result{} -12 and 3\n"
+                      "(euclidean/ -123 10) @result{} -13 and 7\n"
+                      "(euclidean/ -123 -10) @result{} 13 and 7\n"
+                      "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+                      "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_quo_and_rem
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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_euclidean_quo_and_rem);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * yy;
+             if (rr < 0)
+               {
+                 if (yy > 0)
+                   { rr += yy; qq--; }
+                 else
+                   { rr -= yy; qq++; }
+               }
+             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         scm_t_inum xx = SCM_I_INUM (x);
+         if (xx >= 0)
+           return scm_values (scm_list_2 (SCM_INUM0, x));
+         else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             return scm_values
+               (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r)));
+           }
+         else
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r));
+             return scm_values (scm_list_2 (SCM_INUM1, scm_i_normbig (r)));
+           }
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quo_and_rem
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_euclidean_quo_and_rem);
+    }
+  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_euclidean_quo_and_rem);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             SCM r = scm_i_mkbig ();
+             if (yy > 0)
+               mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                               SCM_I_BIG_MPZ (x), yy);
+             else
+               {
+                 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                                 SCM_I_BIG_MPZ (x), -yy);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+               }
+             scm_remember_upto_here_1 (x);
+             return scm_values (scm_list_2 (scm_i_normbig (q),
+                                            scm_i_normbig (r)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM q = scm_i_mkbig ();
+         SCM r = scm_i_mkbig ();
+         if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                        SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+         else
+           mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                        SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_values (scm_list_2 (scm_i_normbig (q),
+                                        scm_i_normbig (r)));
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quo_and_rem
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_euclidean_quo_and_rem);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_quo_and_rem
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_euclidean_quo_and_rem);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quo_and_rem
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG1,
+                       s_scm_euclidean_quo_and_rem);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_quo_and_rem (double x, double y)
+{
+  double q, r;
+
+  if (SCM_LIKELY (y > 0))
+    q = floor (x / y);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x / y);
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_quo_and_rem);  /* or return a NaN? */
+  else
+    q = guile_NaN;
+  r = x - q * y;
+  return scm_values (scm_list_2 (scm_from_double (q),
+                                scm_from_double (r)));
+}
+
+/* Compute exact euclidean quotient and remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_quo_and_rem (SCM x, SCM y)
+{
+  SCM q, r;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG1,
+                       s_scm_euclidean_quo_and_rem);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
+                       s_scm_euclidean_quo_and_rem);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_quo_and_rem);
+  r = scm_difference (x, scm_product (q, y));
+  return scm_values (scm_list_2 (q, r));
+}
+
+static SCM scm_i_inexact_centered_quotient (double x, double y);
+static SCM scm_i_bigint_centered_quotient (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_quotient (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} such that\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
+                      "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+                      "@lisp\n"
+                      "(centered-quotient 123 10) @result{} 12\n"
+                      "(centered-quotient 123 -10) @result{} -12\n"
+                      "(centered-quotient -123 10) @result{} -12\n"
+                      "(centered-quotient -123 -10) @result{} 12\n"
+                      "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
+                      "(centered-quotient 16/3 -10/7) @result{} -4\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_quotient
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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_quotient);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * yy;
+             if (SCM_LIKELY (xx > 0))
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr >= (yy + 1) / 2)
+                       qq++;
+                   }
+                 else
+                   {
+                     if (rr >= (1 - yy) / 2)
+                       qq--;
+                   }
+               }
+             else
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr < -yy / 2)
+                       qq--;
+                   }
+                 else
+                   {
+                     if (rr < yy / 2)
+                       qq++;
+                   }
+               }
+             return SCM_I_MAKINUM (qq);
+           }
+       }
+      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_quotient */
+         return scm_i_bigint_centered_quotient
+           (scm_i_long2big (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                           s_scm_centered_quotient);
+    }
+  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_quotient);
+         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);
+               }
+             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);
+               }
+             return scm_i_normbig (q);
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_quotient (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                           s_scm_centered_quotient);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_centered_quotient
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                           s_scm_centered_quotient);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_quotient (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
+                       s_scm_centered_quotient);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_quotient (double x, double y)
+{
+  if (SCM_LIKELY (y > 0))
+    return scm_from_double (floor (x/y + 0.5));
+  else if (SCM_LIKELY (y < 0))
+    return scm_from_double (ceil (x/y - 0.5));
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_quotient);  /* or return a NaN? */
+  else
+    return scm_nan ();
+}
+
+/* Assumes that both x and y are bigints, though
+   x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_quotient (SCM x, SCM y)
+{
+  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));
+      scm_remember_upto_here_2 (x, 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);
+    }
+  else
+    {
+      mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      scm_remember_upto_here_2 (x, 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);
+    }
+  scm_remember_upto_here_2 (r, min_r);
+  return scm_i_normbig (q);
+}
+
+/* Compute exact centered quotient the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_quotient (SCM x, SCM y)
+{
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
+                       s_scm_centered_quotient);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                       s_scm_centered_quotient);
+  else if (scm_is_true (scm_positive_p (y)))
+    return scm_floor (scm_sum (scm_divide (x, y),
+                              exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    return scm_ceiling (scm_difference (scm_divide (x, y),
+                                       exactly_one_half));
+  else
+    scm_num_overflow (s_scm_centered_quotient);
+}
+
+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_slow_exact_centered_remainder (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the real number @var{r} such that\n"
+                      "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
+                      "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "for some integer @var{q}.\n"
+                      "@lisp\n"
+                      "(centered-remainder 123 10) @result{} 3\n"
+                      "(centered-remainder 123 -10) @result{} 3\n"
+                      "(centered-remainder -123 10) @result{} -3\n"
+                      "(centered-remainder -123 -10) @result{} -3\n"
+                      "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
+                      "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_remainder
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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 xx = SCM_I_INUM (x);
+             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);
+           }
+       }
+      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 (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                           s_scm_centered_remainder);
+    }
+  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);
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_remainder (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                           s_scm_centered_remainder);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_centered_remainder
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                           s_scm_centered_remainder);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_remainder (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
+                       s_scm_centered_remainder);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_remainder (double x, double y)
+{
+  double q;
+
+  /* Although it would be more efficient to use fmod here, we can't
+     because it would in some cases produce results inconsistent with
+     scm_i_inexact_centered_quotient, such that x != r + q * y (not even
+     close).  In particular, when x-y/2 is very close to a multiple of
+     y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
+     two cases must correspond to different choices of q.  If quotient
+     chooses one and remainder chooses the other, it would be bad. */
+  if (SCM_LIKELY (y > 0))
+    q = floor (x/y + 0.5);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x/y - 0.5);
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_remainder);  /* or return a NaN? */
+  else
+    return scm_nan ();
+  return scm_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);
+}
+
+/* Compute exact centered_remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_remainder (SCM x, SCM y)
+{
+  SCM q;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
+                       s_scm_centered_remainder);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                       s_scm_centered_remainder);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_sum (scm_divide (x, y), exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_difference (scm_divide (x, y), exactly_one_half));
+  else
+    scm_num_overflow (s_scm_centered_remainder);
+  return scm_difference (x, scm_product (y, q));
+}
+
+
+static SCM scm_i_inexact_centered_quo_and_rem (double x, double y);
+static SCM scm_i_bigint_centered_quo_and_rem (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_quo_and_rem (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} and the real number @var{r}\n"
+                      "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+                      "@lisp\n"
+                      "(centered/ 123 10) @result{} 12 and 3\n"
+                      "(centered/ 123 -10) @result{} -12 and 3\n"
+                      "(centered/ -123 10) @result{} -12 and -3\n"
+                      "(centered/ -123 -10) @result{} 12 and -3\n"
+                      "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+                      "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_quo_and_rem
+{
+  if (SCM_LIKELY (SCM_I_INUMP (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_quo_and_rem);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * 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; }
+                   }
+               }
+             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      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_quo_and_rem */
+         return scm_i_bigint_centered_quo_and_rem
+           (scm_i_long2big (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quo_and_rem
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quo_and_rem (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_centered_quo_and_rem);
+    }
+  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_quo_and_rem);
+         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;
+                   }
+               }
+             return scm_values (scm_list_2 (scm_i_normbig (q),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_quo_and_rem (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quo_and_rem
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quo_and_rem (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_centered_quo_and_rem);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_centered_quo_and_rem
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+     else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
+                           s_scm_centered_quo_and_rem);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_quo_and_rem
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_quo_and_rem (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG1,
+                       s_scm_centered_quo_and_rem);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_quo_and_rem (double x, double y)
+{
+  double q, r;
+
+  if (SCM_LIKELY (y > 0))
+    q = floor (x/y + 0.5);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x/y - 0.5);
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_quo_and_rem);  /* or return a NaN? */
+  else
+    q = guile_NaN;
+  r = x - q * y;
+  return scm_values (scm_list_2 (scm_from_double (q),
+                                scm_from_double (r)));
+}
+
+/* Assumes that both x and y are bigints, though
+   x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_quo_and_rem (SCM x, SCM y)
+{
+  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);
+  return scm_values (scm_list_2 (scm_i_normbig (q),
+                                scm_i_normbig (r)));
+}
+
+/* Compute exact centered quotient and remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_quo_and_rem (SCM x, SCM y)
+{
+  SCM q, r;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG1,
+                       s_scm_centered_quo_and_rem);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
+                       s_scm_centered_quo_and_rem);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_sum (scm_divide (x, y),
+                           exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_difference (scm_divide (x, y),
+                                    exactly_one_half));
+  else
+    scm_num_overflow (s_scm_centered_quo_and_rem);
+  r = scm_difference (x, scm_product (q, y));
+  return scm_values (scm_list_2 (q, r));
+}
+
+
 SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1,
                        (SCM x, SCM y, SCM rest),
                        "Return the greatest common divisor of all parameter values.\n"
@@ -5356,8 +6581,6 @@ SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-static SCM exactly_one_half;
-
 SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
            (SCM x),
            "Round the number @var{x} towards the nearest integer. "

diff --git a/libguile/numbers.h b/libguile/numbers.h
index 740dc80..76d2972 100644 (file)
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -177,6 +177,12 @@ SCM_API SCM scm_abs (SCM x);
 SCM_API SCM scm_quotient (SCM x, SCM y);
 SCM_API SCM scm_remainder (SCM x, SCM y);
 SCM_API SCM scm_modulo (SCM x, SCM y);
+SCM_API SCM scm_euclidean_quo_and_rem (SCM x, SCM y);
+SCM_API SCM scm_euclidean_quotient (SCM x, SCM y);
+SCM_API SCM scm_euclidean_remainder (SCM x, SCM y);
+SCM_API SCM scm_centered_quo_and_rem (SCM x, SCM y);
+SCM_API SCM scm_centered_quotient (SCM x, SCM y);
+SCM_API SCM scm_centered_remainder (SCM x, SCM y);
 SCM_API SCM scm_gcd (SCM x, SCM y);
 SCM_API SCM scm_lcm (SCM n1, SCM n2);
 SCM_API SCM scm_logand (SCM n1, SCM n2);

diff --git a/module/rnrs/arithmetic/fixnums.scm b/module/rnrs/arithmetic/fixnums.scm
index c1f3571..befbe9d 100644 (file)
--- a/module/rnrs/arithmetic/fixnums.scm
+++ b/module/rnrs/arithmetic/fixnums.scm
@@ -1,6 +1,6 @@
 ;;; fixnums.scm --- The R6RS fixnums arithmetic library
 
-;;      Copyright (C) 2010 Free Software Foundation, Inc.
+;;      Copyright (C) 2010, 2011 Free Software Foundation, Inc.
 ;;
 ;; This library is free software; you can redistribute it and/or
 ;; modify it under the terms of the GNU Lesser General Public
@@ -175,40 +175,33 @@
 
   (define (fxdiv fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
-    (let ((r (div fx1 fx2))) r))
+    (div fx1 fx2))
 
   (define (fxmod fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
-    (let ((r (mod fx1 fx2))) r))
+    (mod fx1 fx2))
 
   (define (fxdiv-and-mod fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
     (div-and-mod fx1 fx2))
 
   (define (fxdiv0 fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
-    (let ((r (div0 fx1 fx2))) r))
+    (div0 fx1 fx2))
   
   (define (fxmod0 fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
-    (let ((r (mod0 fx1 fx2))) r))    
+    (mod0 fx1 fx2))
 
   (define (fxdiv0-and-mod0 fx1 fx2)
     (assert-fixnum fx1 fx2)
-    (if (zero? fx2) (raise (make-assertion-violation)))
-    (call-with-values (lambda () (div0-and-mod0 fx1 fx2))
-      (lambda (q r) (values q r))))
+    (div0-and-mod0 fx1 fx2))
 
   (define (fx+/carry fx1 fx2 fx3)
     (assert-fixnum fx1 fx2 fx3)
     (let* ((s (+ fx1 fx2 fx3))
-          (s0 (mod0 s (inexact->exact (expt 2 (fixnum-width)))))
-          (s1 (div0 s (inexact->exact (expt 2 (fixnum-width))))))
+          (s0 (mod0 s (expt 2 (fixnum-width))))
+          (s1 (div0 s (expt 2 (fixnum-width)))))
       (values s0 s1)))
 
   (define (fx-/carry fx1 fx2 fx3)

diff --git a/module/rnrs/arithmetic/flonums.scm b/module/rnrs/arithmetic/flonums.scm
index 4fadbd0..b65c294 100644 (file)
--- a/module/rnrs/arithmetic/flonums.scm
+++ b/module/rnrs/arithmetic/flonums.scm
@@ -1,6 +1,6 @@
 ;;; flonums.scm --- The R6RS flonums arithmetic library
 
-;;      Copyright (C) 2010 Free Software Foundation, Inc.
+;;      Copyright (C) 2010, 2011 Free Software Foundation, Inc.
 ;;
 ;; This library is free software; you can redistribute it and/or
 ;; modify it under the terms of the GNU Lesser General Public
@@ -127,40 +127,27 @@
 
   (define (fldiv-and-mod fl1 fl2)
     (assert-iflonum fl1 fl2)
-    (if (zero? fl2) (raise (make-assertion-violation)))
-    (let ((fx1 (inexact->exact fl1))
-         (fx2 (inexact->exact fl2)))
-      (call-with-values (lambda () (div-and-mod fx1 fx2))
-       (lambda (div mod) (values (exact->inexact div)
-                                 (exact->inexact mod))))))
+    (div-and-mod fl1 fl2))
 
   (define (fldiv fl1 fl2)
     (assert-iflonum fl1 fl2)
-    (if (zero? fl2) (raise (make-assertion-violation)))
-    (let ((fx1 (inexact->exact fl1))
-         (fx2 (inexact->exact fl2)))
-      (exact->inexact (quotient fx1 fx2))))
+    (div fl1 fl2))
 
   (define (flmod fl1 fl2)
     (assert-iflonum fl1 fl2)
-    (if (zero? fl2) (raise (make-assertion-violation)))
-    (let ((fx1 (inexact->exact fl1))
-         (fx2 (inexact->exact fl2)))
-      (exact->inexact (modulo fx1 fx2))))
+    (mod fl1 fl2))
 
   (define (fldiv0-and-mod0 fl1 fl2)
     (assert-iflonum fl1 fl2)
-    (if (zero? fl2) (raise (make-assertion-violation)))
-    (let* ((fx1 (inexact->exact fl1))
-          (fx2 (inexact->exact fl2)))
-      (call-with-values (lambda () (div0-and-mod0 fx1 fx2))
-       (lambda (q r) (values (real->flonum q) (real->flonum r))))))
+    (div0-and-mod0 fl1 fl2))
 
   (define (fldiv0 fl1 fl2)
-    (call-with-values (lambda () (fldiv0-and-mod0 fl1 fl2)) (lambda (q r) q)))
+    (assert-iflonum fl1 fl2)
+    (div0 fl1 fl2))
 
   (define (flmod0 fl1 fl2)
-    (call-with-values (lambda () (fldiv0-and-mod0 fl1 fl2)) (lambda (q r) r)))
+    (assert-iflonum fl1 fl2)
+    (mod0 fl1 fl2))
 
   (define (flnumerator fl) 
     (assert-flonum fl) 

diff --git a/module/rnrs/base.scm b/module/rnrs/base.scm
index 04a7e23..c81ded1 100644 (file)
--- a/module/rnrs/base.scm
+++ b/module/rnrs/base.scm
@@ -74,8 +74,12 @@
 
          syntax-rules identifier-syntax)
   (import (rename (except (guile) error raise)
-                  (quotient div) 
-                  (modulo mod)
+                  (euclidean-quotient div)
+                  (euclidean-remainder mod)
+                  (euclidean/ div-and-mod)
+                  (centered-quotient div0)
+                  (centered-remainder mod0)
+                  (centered/ div0-and-mod0)
                   (inf? infinite?)
                   (exact->inexact inexact)
                   (inexact->exact exact))
@@ -119,21 +123,6 @@
  (define (vector-map proc . vecs)
    (list->vector (apply map (cons proc (map vector->list vecs)))))
 
- (define (div-and-mod x y) (let ((q (div x y)) (r (mod x y))) (values q r)))
-
- (define (div0 x y)
-   (call-with-values (lambda () (div0-and-mod0 x y)) (lambda (q r) q)))
-
- (define (mod0 x y)
-   (call-with-values (lambda () (div0-and-mod0 x y)) (lambda (q r) r)))
-
- (define (div0-and-mod0 x y)
-   (call-with-values (lambda () (div-and-mod x y))
-     (lambda (q r)
-       (cond ((< r (abs (/ y 2))) (values q r))
-            ((negative? y) (values (- q 1) (+ r y)))
-            (else (values (+ q 1) (+ r y)))))))
-
  (define raise
    (@ (rnrs exceptions) raise))
  (define condition

diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 36e3128..9cf9202 100644 (file)
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -17,7 +17,8 @@
 
 (define-module (test-suite test-numbers)
   #:use-module (test-suite lib)
-  #:use-module (ice-9 documentation))
+  #:use-module (ice-9 documentation)
+  #:use-module (srfi srfi-11))  ; let-values
 
 ;;;
 ;;; miscellaneous
@@ -92,6 +93,35 @@
        (negative? obj)
        (inf? obj)))
 
+;;
+;; Tolerance used by test-eqv? for inexact numbers.
+;;
+(define test-epsilon 1e-10)
+
+;;
+;; Like eqv?, except that inexact finite numbers need only be within
+;; test-epsilon (1e-10) to be considered equal.  An exception is made
+;; for zeroes, however.  If X is zero, then it is tested using eqv?
+;; without any allowance for imprecision.  In particular, 0.0 is
+;; considered distinct from -0.0.  For non-real complex numbers,
+;; each component is tested according to these rules.  The intent
+;; is that the known-correct value will be the first parameter.
+;;
+(define (test-eqv? x y)
+  (cond ((real? x)
+        (and (real? y) (test-real-eqv? x y)))
+       ((complex? x)
+        (and (not (real? y))
+             (test-real-eqv? (real-part x) (real-part y))
+             (test-real-eqv? (imag-part x) (imag-part y))))
+       (else (eqv? x y))))
+
+;; Auxiliary predicate used by test-eqv?
+(define (test-real-eqv? x y)
+  (cond ((or (exact? x) (zero? x) (nan? x) (inf? x))
+        (eqv? x y))
+       (else (and (inexact? y) (> test-epsilon (abs (- x y)))))))
+
 (define const-e    2.7182818284590452354)
 (define const-e^2  7.3890560989306502274)
 (define const-1/e  0.3678794411714423215)
@@ -3480,3 +3510,144 @@
   (pass-if "-100i swings back to 45deg down"
     (eqv-loosely? +7.071-7.071i (sqrt -100.0i))))
 
+;;;
+;;; euclidean/
+;;; euclidean-quotient
+;;; euclidean-remainder
+;;; centered/
+;;; centered-quotient
+;;; centered-remainder
+;;;
+
+(with-test-prefix "Number-theoretic division"
+
+  ;; Tests that (lo <= x < hi),
+  ;; but allowing for imprecision
+  ;; if x is inexact.
+  (define (test-within-range? lo hi x)
+    (if (exact? x)
+        (and (<= lo x) (< x hi))
+        (let ((lo (- lo test-epsilon))
+              (hi (+ hi test-epsilon)))
+          (<= lo x hi))))
+
+  (define (safe-euclidean-quotient x y)
+    (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
+          ((zero? y) (throw 'divide-by-zero))
+          ((nan?  y) (nan))
+          ((positive? y) (floor   (/ x y)))
+          ((negative? y) (ceiling (/ x y)))
+          (else (throw 'unknown-problem))))
+
+  (define (safe-euclidean-remainder x y)
+    (- x (* y (safe-euclidean-quotient x y))))
+
+  (define (safe-euclidean/ x y)
+    (let ((q (safe-euclidean-quotient x y))
+          (r (safe-euclidean-remainder x y)))
+      (if (not (and (eq? (exact? q) (exact? r))
+                    (eq? (exact? q) (and (exact? x) (exact? y)))
+                    (test-real-eqv? r (- x (* q y)))
+                    (or (and (integer? q)
+                             (test-within-range? 0 (abs y) r))
+                        (not (finite? x))
+                        (not (finite? y)))))
+          (throw 'safe-euclidean/-is-broken (list x y q r))
+          (values q r))))
+
+  (define (safe-centered-quotient x y)
+    (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
+          ((zero? y) (throw 'divide-by-zero))
+          ((nan?  y) (nan))
+          ((positive? y) (floor   (+  1/2 (/ x y))))
+          ((negative? y) (ceiling (+ -1/2 (/ x y))))
+          (else (throw 'unknown-problem))))
+
+  (define (safe-centered-remainder x y)
+    (- x (* y (safe-centered-quotient x y))))
+
+  (define (safe-centered/ x y)
+    (let ((q (safe-centered-quotient x y))
+          (r (safe-centered-remainder x y)))
+      (if (not (and (eq? (exact? q) (exact? r))
+                    (eq? (exact? q) (and (exact? x) (exact? y)))
+                    (test-real-eqv? r (- x (* q y)))
+                    (or (and (integer? q)
+                             (test-within-range? (* -1/2 (abs y))
+                                                 (* +1/2 (abs y))
+                                                 r))
+                        (not (finite? x))
+                        (not (finite? y)))))
+          (throw 'safe-centered/-is-broken (list x y q r))
+          (values q r))))
+
+  (define test-numerators
+    (append
+     (list  123  125  127  130  3  5  10  123.2  125.0
+            -123 -125 -127 -130 -3 -5 -10 -123.2 -125.0
+            127.2  130.0  123/7  125/7  127/7  130/7
+            -127.2 -130.0 -123/7 -125/7 -127/7 -130/7
+            0 +0.0 -0.0 +inf.0 -inf.0 +nan.0
+            most-negative-fixnum (1+ most-positive-fixnum)
+            (1- most-negative-fixnum))
+     (apply append
+            (map (lambda (x) (list (* x (+ 1 most-positive-fixnum))
+                                   (* x (+ 2 most-positive-fixnum))))
+                 '( 123  125  127  130  3  5  10
+                   -123 -125 -127 -130 -3 -5 -10)))))
+
+  (define test-denominators
+    (list   10  5  10/7  127/2  10.0  63.5
+            -10 -5 -10/7 -127/2 -10.0 -63.5
+            +inf.0 -inf.0 +nan.0 most-negative-fixnum
+            (+ 1 most-positive-fixnum) (+ -1 most-negative-fixnum)
+            (+ 2 most-positive-fixnum) (+ -2 most-negative-fixnum)))
+
+  (define (do-tests-1 op-name real-op safe-op)
+    (for-each (lambda (d)
+                (for-each (lambda (n)
+                            (run-test (list op-name n d) #t
+                                      (lambda ()
+                                        (test-eqv? (real-op n d)
+                                                   (safe-op n d)))))
+                          test-numerators))
+              test-denominators))
+
+  (define (do-tests-2 op-name real-op safe-op)
+    (for-each (lambda (d)
+                (for-each (lambda (n)
+                            (run-test (list op-name n d) #t
+                                      (lambda ()
+                                        (let-values
+                                            (((q  r)  (safe-op n d))
+                                             ((q1 r1) (real-op n d)))
+                                          (and (test-eqv? q q1)
+                                               (test-eqv? r r1))))))
+                          test-numerators))
+              test-denominators))
+
+  (with-test-prefix "euclidean-quotient"
+    (do-tests-1 'euclidean-quotient
+                euclidean-quotient
+                safe-euclidean-quotient))
+  (with-test-prefix "euclidean-remainder"
+    (do-tests-1 'euclidean-remainder
+                euclidean-remainder
+                safe-euclidean-remainder))
+  (with-test-prefix "euclidean/"
+    (do-tests-2 'euclidean/
+                euclidean/
+                safe-euclidean/))
+
+  (with-test-prefix "centered-quotient"
+    (do-tests-1 'centered-quotient
+                centered-quotient
+                safe-centered-quotient))
+  (with-test-prefix "centered-remainder"
+    (do-tests-1 'centered-remainder
+                centered-remainder
+                safe-centered-remainder))
+  (with-test-prefix "centered/"
+    (do-tests-2 'centered/
+                centered/
+                safe-centered/)))

diff --git a/test-suite/tests/r6rs-arithmetic-fixnums.test b/test-suite/tests/r6rs-arithmetic-fixnums.test
index fed72eb..d39d544 100644 (file)
--- a/test-suite/tests/r6rs-arithmetic-fixnums.test
+++ b/test-suite/tests/r6rs-arithmetic-fixnums.test
@@ -1,6 +1,6 @@
 ;;; arithmetic-fixnums.test --- Test suite for R6RS (rnrs arithmetic bitwise)
 
-;;      Copyright (C) 2010 Free Software Foundation, Inc.
+;;      Copyright (C) 2010, 2011 Free Software Foundation, Inc.
 ;;
 ;; This library is free software; you can redistribute it and/or
 ;; modify it under the terms of the GNU Lesser General Public
@@ -121,32 +121,25 @@
   (pass-if "simple"
     (call-with-values (lambda () (fxdiv-and-mod 123 10))
       (lambda (d m) 
-       (or (and (fx=? d 12) (fx=? m 3))
-           (throw 'unresolved))))))
+       (and (fx=? d 12) (fx=? m 3))))))
 
-(with-test-prefix "fxdiv"
-  (pass-if "simple" (or (fx=? (fxdiv -123 10) -13) (throw 'unresolved))))
-
-(with-test-prefix "fxmod"
-  (pass-if "simple" (or (fx=? (fxmod -123 10) 7) (throw 'unresolved))))
+(with-test-prefix "fxdiv" (pass-if "simple" (fx=? (fxdiv -123 10) -13)))
+(with-test-prefix "fxmod" (pass-if "simple" (fx=? (fxmod -123 10) 7)))
 
 (with-test-prefix "fxdiv0-and-mod0"
   (pass-if "simple"
     (call-with-values (lambda () (fxdiv0-and-mod0 -123 10))
       (lambda (d m)
-       (or (and (fx=? d 12) (fx=? m -3))
-           (throw 'unresolved))))))
-
-(with-test-prefix "fxdiv0"
-  (pass-if "simple" (or (fx=? (fxdiv0 -123 10) 12) (throw 'unresolved))))
+       (and (fx=? d -12) (fx=? m -3))))))
 
-(with-test-prefix "fxmod0"
-  (pass-if "simple" (or (fx=? (fxmod0 -123 10) -3) (throw 'unresolved))))
+(with-test-prefix "fxdiv0" (pass-if "simple" (fx=? (fxdiv0 -123 10) -12)))
+(with-test-prefix "fxmod0" (pass-if "simple" (fx=? (fxmod0 -123 10) -3)))
 
 
 ;; Without working div and mod implementations and without any example results
 ;; from the spec, I have no idea what the results of these functions should
 ;; be.  -juliang
+;; UPDATE: div and mod implementations are now working properly  -mhw
 
 (with-test-prefix "fx+/carry" (pass-if "simple" (throw 'unresolved)))
 

diff --git a/test-suite/tests/r6rs-arithmetic-flonums.test b/test-suite/tests/r6rs-arithmetic-flonums.test
index 873447b..af9dbbf 100644 (file)
--- a/test-suite/tests/r6rs-arithmetic-flonums.test
+++ b/test-suite/tests/r6rs-arithmetic-flonums.test
@@ -1,6 +1,6 @@
 ;;; arithmetic-flonums.test --- Test suite for R6RS (rnrs arithmetic flonums)
 
-;;      Copyright (C) 2010 Free Software Foundation, Inc.
+;;      Copyright (C) 2010, 2011 Free Software Foundation, Inc.
 ;;
 ;; This library is free software; you can redistribute it and/or
 ;; modify it under the terms of the GNU Lesser General Public
@@ -195,14 +195,13 @@
   (pass-if "simple"
     (call-with-values (lambda () (fldiv0-and-mod0 -123.0 10.0))
       (lambda (div mod) 
-       (or (and (fl=? div -12.0) (fl=? mod -3.0))
-           (throw 'unresolved))))))
+       (and (fl=? div -12.0) (fl=? mod -3.0))))))
 
 (with-test-prefix "fldiv0" 
-  (pass-if "simple" (or (fl=? (fldiv0 -123.0 10.0) -12.0) (throw 'unresolved))))
+  (pass-if "simple" (fl=? (fldiv0 -123.0 10.0) -12.0)))
 
 (with-test-prefix "flmod0" 
-  (pass-if "simple" (or (fl=? (flmod0 -123.0 10.0) -3.0) (throw 'unresolved))))
+  (pass-if "simple" (fl=? (flmod0 -123.0 10.0) -3.0)))
 
 (with-test-prefix "flnumerator"
   (pass-if "simple" (fl=? (flnumerator 0.5) 1.0))