* Re: Simplified code for fast R6RS division operators
[not found] <87bp2y14zu.fsf@yeeloong.netris.org>
@ 2011-01-30 21:36 ` Andy Wingo
0 siblings, 0 replies; only message in thread
From: Andy Wingo @ 2011-01-30 21:36 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 732 bytes --]
On Sun 30 Jan 2011 22:00, Mark H Weaver <mhw@netris.org> writes:
> I'm sending this to you directly instead of the list,
> because I've probably been posting too much :)
The list is still probably the right thing :) It's not too much traffic
for people that care about development _of_ guile. People who just care
about development _with_ guile have guile-users :)
> Regarding my code for the fast R6RS division operators: I realized in
> the last couple hours that I could cut out some code and simplify other
> bits. So here's a new slightly simpler version of that patch in case
> you haven't yet started reviewed the code.
Thanks for the note, I'll use the new patch (attached, for the list's
benefit :).
Cheers,
Andy
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: r6rs division patch --]
[-- Type: text/x-patch, Size: 69171 bytes --]
From 235e4643f63c6d0ec105d8eb606c04c5c3168f94 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 30 Jan 2011 08:48:28 -0500
Subject: [PATCH] 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.
---
NEWS | 29 +
doc/ref/api-data.texi | 79 ++
doc/ref/r6rs.texi | 70 ++-
libguile/numbers.c | 1227 ++++++++++++++++++++++++-
libguile/numbers.h | 6 +
module/rnrs/arithmetic/fixnums.scm | 23 +-
module/rnrs/arithmetic/flonums.scm | 31 +-
module/rnrs/base.scm | 23 +-
test-suite/tests/numbers.test | 173 ++++-
test-suite/tests/r6rs-arithmetic-fixnums.test | 23 +-
test-suite/tests/r6rs-arithmetic-flonums.test | 9 +-
11 files changed, 1603 insertions(+), 90 deletions(-)
diff --git a/NEWS b/NEWS
index f45795e..33e49a4 100644
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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))
--
1.5.6.5
[-- Attachment #3: Type: text/plain, Size: 26 bytes --]
--
http://wingolog.org/
^ permalink raw reply related [flat|nested] only message in thread
only message in thread, other threads:[~2011-01-30 21:36 UTC | newest]
Thread overview: (only message) (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
[not found] <87bp2y14zu.fsf@yeeloong.netris.org>
2011-01-30 21:36 ` Simplified code for fast R6RS division operators Andy Wingo
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).