[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Guile-commits] GNU Guile branch, master, updated. release_1-9-14-145-gf
|
From: |
Andy Wingo |
|
Subject: |
[Guile-commits] GNU Guile branch, master, updated. release_1-9-14-145-gff62c16 |
|
Date: |
Sun, 30 Jan 2011 22:00:39 +0000 |
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU Guile".
http://git.savannah.gnu.org/cgit/guile.git/commit/?id=ff62c16828d41955455805dd1b427966944b7d27
The branch, master has been updated
via ff62c16828d41955455805dd1b427966944b7d27 (commit)
via a16982ca4fafe7d68907c5bf2aae4eb627a1a448 (commit)
from 074c414e297ba4095d2398fac17842c56ad24b11 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit ff62c16828d41955455805dd1b427966944b7d27
Author: Mark H Weaver <address@hidden>
Date: Sun Jan 30 08:48:28 2011 -0500
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.
commit a16982ca4fafe7d68907c5bf2aae4eb627a1a448
Author: Mark H Weaver <address@hidden>
Date: Fri Jan 28 23:58:02 2011 -0500
Add SCM_LIKELY and SCM_UNLIKELY for optimization
* libguile/numbers.c (scm_abs, scm_quotient, scm_remainder, scm_modulo):
Add SCM_LIKELY and SCM_UNLIKELY in several places for optimization.
(scm_remainder): Add comment about C99 "%" semantics.
Strip away a redundant set of braces.
-----------------------------------------------------------------------
Summary of changes:
NEWS | 29 +
doc/ref/api-data.texi | 79 ++
doc/ref/r6rs.texi | 70 ++-
libguile/numbers.c | 1292 ++++++++++++++++++++++++-
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, 1636 insertions(+), 122 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
address@hidden euclidean/
address@hidden euclidean-quotient
address@hidden euclidean-remainder
address@hidden centered/
address@hidden centered-quotient
address@hidden 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
address@hidden {Scheme Procedure} euclidean/ x y
address@hidden {Scheme Procedure} euclidean-quotient x y
address@hidden {Scheme Procedure} euclidean-remainder x y
address@hidden {C Function} scm_euclidean_quo_and_rem (x y)
address@hidden {C Function} scm_euclidean_quotient (x y)
address@hidden {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
address@hidden such that @address@hidden = @address@hidden + @var{r}} and
address@hidden <= @var{r} < abs(@var{y})}. @code{euclidean/} returns both
@var{q} and
address@hidden, and is more efficient than computing each separately. Note
+that when @address@hidden > 0}, @code{euclidean-quotient} returns
address@hidden(@var{x}/@var{y})}, otherwise it returns
address@hidden(@var{x}/@var{y})}.
+
+Note that these operators are equivalent to the R6RS operators
address@hidden, @code{mod}, and @code{div-and-mod}.
+
address@hidden
+(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
address@hidden lisp
address@hidden deffn
+
address@hidden {Scheme Procedure} centered/ x y
address@hidden {Scheme Procedure} centered-quotient x y
address@hidden {Scheme Procedure} centered-remainder x y
address@hidden {C Function} scm_centered_quo_and_rem (x y)
address@hidden {C Function} scm_centered_quotient (x y)
address@hidden {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
address@hidden such that @address@hidden = @address@hidden + @var{r}} and
address@hidden(@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 @address@hidden/@var{y}}
+rounded to the nearest integer. When @address@hidden/@var{y}} lies
+exactly half-way between two integers, the tie is broken according to
+the sign of @var{y}. If @address@hidden > 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
address@hidden, @code{mod0}, and @code{div0-and-mod0}.
+
address@hidden
+(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
address@hidden lisp
address@hidden 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.
address@hidden Copyright (C) 2010
address@hidden 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
address@hidden {Scheme Procedure} div x1 x2
address@hidden {Scheme Procedure} mod x1 x2
address@hidden {Scheme Procedure} div-and-mod x1 x2
-These procedures implement number-theoretic division.
address@hidden div
address@hidden mod
address@hidden div-and-mod
address@hidden {Scheme Procedure} div x y
address@hidden {Scheme Procedure} mod x y
address@hidden {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
address@hidden@var{x} = @address@hidden + @var{r}} and @math{0 <= @var{r} <
abs(@var{y})}.
address@hidden returns both @var{q} and @var{r}, and is more
+efficient than computing each separately. Note that when @address@hidden > 0},
address@hidden returns @math{floor(@var{x}/@var{y})}, otherwise
+it returns @math{ceiling(@var{x}/@var{y})}.
address@hidden returns two values, the respective results of
address@hidden(div x1 x2)} and @code{(mod x1 x2)}.
address@hidden
+(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
address@hidden lisp
@end deffn
address@hidden {Scheme Procedure} div0 x1 x2
address@hidden {Scheme Procedure} mod0 x1 x2
address@hidden {Scheme Procedure} div0-and-mod0 x1 x2
-These procedures are similar to @code{div}, @code{mod}, and
address@hidden, except that @code{mod0} returns values that lie
-within a half-open interval centered on zero.
address@hidden div0
address@hidden mod0
address@hidden div0-and-mod0
address@hidden {Scheme Procedure} div0 x y
address@hidden {Scheme Procedure} mod0 x y
address@hidden {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
address@hidden such that @address@hidden = @address@hidden + @var{r}} and
address@hidden(@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 @address@hidden/@var{y}} rounded to the
+nearest integer. When @address@hidden/@var{y}} lies exactly half-way
+between two integers, the tie is broken according to the sign of
address@hidden If @address@hidden > 0}, ties are rounded toward positive
+infinity, otherwise they are rounded toward negative infinity.
+This is a consequence of the requirement that
address@hidden(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.
+
address@hidden
+(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
address@hidden lisp
@end deffn
@deffn {Scheme Procedure} exact-integer-sqrt k
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 608cf7a..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)
@@ -774,18 +775,18 @@ SCM_GPROC (s_quotient, "quotient", 2, 0, 0, scm_quotient,
g_quotient);
SCM
scm_quotient (SCM x, SCM y)
{
- if (SCM_I_INUMP (x))
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
{
scm_t_inum xx = SCM_I_INUM (x);
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_quotient);
else
{
scm_t_inum z = xx / yy;
- if (SCM_FIXABLE (z))
+ if (SCM_LIKELY (SCM_FIXABLE (z)))
return SCM_I_MAKINUM (z);
else
return scm_i_inum2big (z);
@@ -809,12 +810,12 @@ scm_quotient (SCM x, SCM y)
}
else if (SCM_BIGP (x))
{
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_quotient);
- else if (yy == 1)
+ else if (SCM_UNLIKELY (yy == 1))
return x;
else
{
@@ -858,15 +859,18 @@ SCM_GPROC (s_remainder, "remainder", 2, 0, 0,
scm_remainder, g_remainder);
SCM
scm_remainder (SCM x, SCM y)
{
- if (SCM_I_INUMP (x))
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
{
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_remainder);
else
{
+ /* C99 specifies that "%" is the remainder corresponding to a
+ quotient rounded towards zero, and that's also traditional
+ for machine division, so z here should be well defined. */
scm_t_inum z = SCM_I_INUM (x) % yy;
return SCM_I_MAKINUM (z);
}
@@ -889,10 +893,10 @@ scm_remainder (SCM x, SCM y)
}
else if (SCM_BIGP (x))
{
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_remainder);
else
{
@@ -931,13 +935,13 @@ SCM_GPROC (s_modulo, "modulo", 2, 0, 0, scm_modulo,
g_modulo);
SCM
scm_modulo (SCM x, SCM y)
{
- if (SCM_I_INUMP (x))
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
{
scm_t_inum xx = SCM_I_INUM (x);
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_modulo);
else
{
@@ -1008,10 +1012,10 @@ scm_modulo (SCM x, SCM y)
}
else if (SCM_BIGP (x))
{
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum yy = SCM_I_INUM (y);
- if (yy == 0)
+ if (SCM_UNLIKELY (yy == 0))
scm_num_overflow (s_modulo);
else
{
@@ -1029,22 +1033,20 @@ scm_modulo (SCM x, SCM y)
}
else if (SCM_BIGP (y))
{
- {
- SCM result = scm_i_mkbig ();
- int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
- SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
- mpz_mod (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (pos_y));
+ SCM result = scm_i_mkbig ();
+ int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
+ SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
+ mpz_mod (SCM_I_BIG_MPZ (result),
+ SCM_I_BIG_MPZ (x),
+ SCM_I_BIG_MPZ (pos_y));
- scm_remember_upto_here_1 (x);
- if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
- mpz_add (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (y),
- SCM_I_BIG_MPZ (result));
- scm_remember_upto_here_2 (y, pos_y);
- return scm_i_normbig (result);
- }
+ scm_remember_upto_here_1 (x);
+ if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
+ mpz_add (SCM_I_BIG_MPZ (result),
+ SCM_I_BIG_MPZ (y),
+ SCM_I_BIG_MPZ (result));
+ scm_remember_upto_here_2 (y, pos_y);
+ return scm_i_normbig (result);
}
else
SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
@@ -1053,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"
+ "@address@hidden = @address@hidden + @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"
+ "@address@hidden = @address@hidden + @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 @address@hidden = @address@hidden +
@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"
+ "@address@hidden = @address@hidden + @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 @address@hidden = @address@hidden + @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 @address@hidden = @address@hidden +
@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"
@@ -5355,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))
hooks/post-receive
--
GNU Guile
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Guile-commits] GNU Guile branch, master, updated. release_1-9-14-145-gff62c16,
Andy Wingo <=