* [PATCH] Fast R6RS div/mod; improved extensibility of numerics
@ 2011-01-30 16:27 Mark H Weaver
2011-01-30 22:24 ` Andy Wingo
0 siblings, 1 reply; 12+ messages in thread
From: Mark H Weaver @ 2011-01-30 16:27 UTC (permalink / raw)
To: guile-devel
[-- Attachment #1: Type: text/plain, Size: 1728 bytes --]
Here are another four patches of numerics improvements.
The second patch adds fast C implementations of the R6RS
number-theoretic division operations: `div-and-mod', `div', `mod',
`div0-and-mod0', `div0' and `mod0'. Note that these operations will
accept any real numbers (except division by zero), not just integers.
The first three are equivalent to `euclidean/' et al from Taylor
Campbell's proposal, and those are the names I gave them in guile core.
The other three I named `centered/', `centered-quotient', and
`centered-remainder', because the range of their remainder operation is
centered on 0 in a twos-complement sort of way. They are almost the
same as Taylor's `round/' et al, but in my opinion they have superior
properties, as I explained in a recent commentary on this list.
The last two patches improve the extensibility of numeric operations,
and as a side effect add documentation strings to many of them.
The last patch might be slightly controversial. Although it does not
make `integer-expt' a generic function, nonetheless it can now
exponentiate _any_ scheme object that can be multiplied using `*'.
raising to a negative exponent, the reciprocal operation must also be
supported. As a side effect, (expt <obj> 0) will return 1 and
(expt <obj> 1) will return <obj> for any scheme object whatsoever,
without operating on it at all. Raising to any other power will however
throw an exception if the requisite methods have not been defined.
I've attached a simple demo that I wrote in about 10 minutes, which
defines a bare-bones matrix class and uses `expt' to do fast matrix
exponentiation. The demo uses this to define a very fast and exact
fibonacci procedure (fib n).
Best,
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Add SCM_LIKELY and SCM_UNLIKELY for optimization --]
[-- Type: text/x-diff, Size: 4706 bytes --]
From 9f169efcd0fda01983f91e9b599528ab7252607a Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Fri, 28 Jan 2011 23:58:02 -0500
Subject: [PATCH] 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.
---
libguile/numbers.c | 65 ++++++++++++++++++++++++++-------------------------
1 files changed, 33 insertions(+), 32 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 608cf7a..e15a6ac 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -774,18 +774,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 +809,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 +858,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 +892,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 +934,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 +1011,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 +1032,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);
--
1.5.6.5
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: Add two new sets of fast quotient and remainder operators --]
[-- Type: text/x-diff, Size: 69938 bytes --]
From bd09baef686122598f423b594c6dc82e92a3db3f 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 | 1245 ++++++++++++++++++++++++-
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, 1621 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..875cd87 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,1248 @@ 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)))
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+ s_scm_euclidean_quotient);
+ else
+ return scm_i_inexact_euclidean_quotient
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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 should we 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 == SCM_MOST_NEGATIVE_FIXNUM) &&
+ (0 == mpz_cmp_ui (SCM_I_BIG_MPZ (y),
+ - SCM_MOST_NEGATIVE_FIXNUM)))
+ {
+ /* Special case: x == fixnum-min && y == abs (fixnum-min) */
+ scm_remember_upto_here_1 (y);
+ return SCM_INUM0;
+ }
+ else 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)))
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+ s_scm_euclidean_remainder);
+ else
+ return scm_i_inexact_euclidean_remainder
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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)
+{
+ 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)))
+ return scm_difference
+ (x, scm_product (y, scm_floor (scm_divide (x, y))));
+ else if (scm_is_true (scm_negative_p (y)))
+ return scm_difference
+ (x, scm_product (y, scm_ceiling (scm_divide (x, y))));
+ else
+ scm_num_overflow (s_scm_euclidean_remainder);
+}
+
+
+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 ((xx == SCM_MOST_NEGATIVE_FIXNUM) &&
+ (0 == mpz_cmp_ui (SCM_I_BIG_MPZ (y),
+ - SCM_MOST_NEGATIVE_FIXNUM)))
+ {
+ /* Special case: x == fixnum-min && y == abs (fixnum-min) */
+ scm_remember_upto_here_1 (y);
+ return scm_values
+ (scm_list_2 (SCM_I_MAKINUM (-1), SCM_INUM0));
+ }
+ 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)))
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
+ s_scm_euclidean_quo_and_rem);
+ else
+ return scm_i_inexact_euclidean_quo_and_rem
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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)))
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+ s_scm_centered_quotient);
+ else
+ return scm_i_inexact_centered_quotient
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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)))
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+ s_scm_centered_remainder);
+ else
+ return scm_i_inexact_centered_remainder
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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 with different choices of q. If
+ centered_quotient chooses one way and centered_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)
+{
+ 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)))
+ return scm_difference
+ (x, scm_product (y, scm_floor (scm_sum (scm_divide (x, y),
+ exactly_one_half))));
+ else if (scm_is_true (scm_negative_p (y)))
+ return scm_difference
+ (x, scm_product (y, scm_ceiling (scm_difference (scm_divide (x, y),
+ exactly_one_half))));
+ else
+ scm_num_overflow (s_scm_centered_remainder);
+}
+
+
+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)))
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
+ s_scm_centered_quo_and_rem);
+ else
+ return scm_i_inexact_centered_quo_and_rem
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ }
+ 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 +6599,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
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: Improve extensibility of core numeric procedures --]
[-- Type: text/x-diff, Size: 33895 bytes --]
From 6e9ee1bdbd3a449d7d864d52a42cbab630d57ca1 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 30 Jan 2011 09:52:51 -0500
Subject: [PATCH] Improve extensibility of core numeric procedures
* libguile/numbers.c (scm_quotient, scm_remainder, scm_modulo,
scm_zero_p, scm_positive_p, scm_negative_p, scm_real_part,
scm_imag_part, scm_numerator, scm_denominator, scm_magnitude,
scm_angle, scm_exact_to_inexact): Change from SCM_GPROC to
SCM_PRIMITIVE_GENERIC. As a side effect, all of these procedures now
have documentation strings.
(scm_exact_p, scm_inexact_p, scm_odd_p, scm_even_p, scm_finite_p,
scm_inf_p, scm_nan_p, scm_expt, scm_inexact_to_exact, scm_log,
scm_log10, scm_exp, scm_sqrt): Change from SCM_DEFINE to
SCM_PRIMITIVE_GENERIC, and make sure the code allows these functions
to be extended in practice.
(scm_real_part, scm_imag_part, scm_numerator, scm_denominator,
scm_inexact_to_exact): Simplify type dispatch code.
(scm_sqrt): Rename formal argument from x to z, since complex numbers
are supported.
(scm_abs): Fix empty FUNC_NAME.
* libguile/numbers.h (scm_finite_p): Add missing prototype.
(scm_inf_p, scm_nan_p): Rename formal parameter from n to x, since
the domain is the real numbers.
* test-suite/tests/numbers.test: Test for documentation strings. Change
from `expect-fail' to `pass-if' for several of these, and add tests
for others. Also add other tests for `real-part' and `imag-part',
which previously had none.
---
libguile/numbers.c | 416 ++++++++++++++++++++---------------------
libguile/numbers.h | 5 +-
test-suite/tests/numbers.test | 69 +++++--
3 files changed, 257 insertions(+), 233 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 875cd87..955404d 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -498,8 +498,8 @@ scm_i_fraction2double (SCM z)
SCM_FRACTION_DENOMINATOR (z)));
}
-SCM_DEFINE (scm_exact_p, "exact?", 1, 0, 0,
- (SCM x),
+SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0,
+ (SCM x),
"Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
"otherwise.")
#define FUNC_NAME s_scm_exact_p
@@ -509,12 +509,12 @@ SCM_DEFINE (scm_exact_p, "exact?", 1, 0, 0,
else if (SCM_NUMBERP (x))
return SCM_BOOL_T;
else
- SCM_WRONG_TYPE_ARG (1, x);
+ SCM_WTA_DISPATCH_1 (g_scm_exact_p, x, 1, s_scm_exact_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_inexact_p, "inexact?", 1, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_inexact_p, "inexact?", 1, 0, 0,
(SCM x),
"Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
"else.")
@@ -525,12 +525,12 @@ SCM_DEFINE (scm_inexact_p, "inexact?", 1, 0, 0,
else if (SCM_NUMBERP (x))
return SCM_BOOL_F;
else
- SCM_WRONG_TYPE_ARG (1, x);
+ SCM_WTA_DISPATCH_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0,
(SCM n),
"Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
"otherwise.")
@@ -547,25 +547,24 @@ SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
scm_remember_upto_here_1 (n);
return scm_from_bool (odd_p);
}
- else if (scm_is_true (scm_inf_p (n)))
- SCM_WRONG_TYPE_ARG (1, n);
else if (SCM_REALP (n))
{
- double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
- if (rem == 1.0)
- return SCM_BOOL_T;
- else if (rem == 0.0)
- return SCM_BOOL_F;
- else
- SCM_WRONG_TYPE_ARG (1, n);
+ double val = SCM_REAL_VALUE (n);
+ if (DOUBLE_IS_FINITE (val))
+ {
+ double rem = fabs (fmod (val, 2.0));
+ if (rem == 1.0)
+ return SCM_BOOL_T;
+ else if (rem == 0.0)
+ return SCM_BOOL_F;
+ }
}
- else
- SCM_WRONG_TYPE_ARG (1, n);
+ SCM_WTA_DISPATCH_1 (g_scm_odd_p, n, 1, s_scm_odd_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_even_p, "even?", 1, 0, 0,
(SCM n),
"Return @code{#t} if @var{n} is an even number, @code{#f}\n"
"otherwise.")
@@ -582,25 +581,24 @@ SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
scm_remember_upto_here_1 (n);
return scm_from_bool (even_p);
}
- else if (scm_is_true (scm_inf_p (n)))
- SCM_WRONG_TYPE_ARG (1, n);
else if (SCM_REALP (n))
{
- double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
- if (rem == 1.0)
- return SCM_BOOL_F;
- else if (rem == 0.0)
- return SCM_BOOL_T;
- else
- SCM_WRONG_TYPE_ARG (1, n);
+ double val = SCM_REAL_VALUE (n);
+ if (DOUBLE_IS_FINITE (val))
+ {
+ double rem = fabs (fmod (val, 2.0));
+ if (rem == 1.0)
+ return SCM_BOOL_F;
+ else if (rem == 0.0)
+ return SCM_BOOL_T;
+ }
}
- else
- SCM_WRONG_TYPE_ARG (1, n);
+ SCM_WTA_DISPATCH_1 (g_scm_even_p, n, 1, s_scm_even_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_finite_p, "finite?", 1, 0, 0,
- (SCM x),
+SCM_PRIMITIVE_GENERIC (scm_finite_p, "finite?", 1, 0, 0,
+ (SCM x),
"Return @code{#t} if the real number @var{x} is neither\n"
"infinite nor a NaN, @code{#f} otherwise.")
#define FUNC_NAME s_scm_finite_p
@@ -610,14 +608,14 @@ SCM_DEFINE (scm_finite_p, "finite?", 1, 0, 0,
else if (scm_is_real (x))
return SCM_BOOL_T;
else
- SCM_WRONG_TYPE_ARG (1, x);
+ SCM_WTA_DISPATCH_1 (g_scm_finite_p, x, 1, s_scm_finite_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0,
- (SCM x),
- "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
- "@samp{-inf.0}. Otherwise return @code{#f}.")
+SCM_PRIMITIVE_GENERIC (scm_inf_p, "inf?", 1, 0, 0,
+ (SCM x),
+ "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
+ "@samp{-inf.0}. Otherwise return @code{#f}.")
#define FUNC_NAME s_scm_inf_p
{
if (SCM_REALP (x))
@@ -625,12 +623,12 @@ SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0,
else if (scm_is_real (x))
return SCM_BOOL_F;
else
- SCM_WRONG_TYPE_ARG (1, x);
+ SCM_WTA_DISPATCH_1 (g_scm_inf_p, x, 1, s_scm_inf_p);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0,
- (SCM x),
+SCM_PRIMITIVE_GENERIC (scm_nan_p, "nan?", 1, 0, 0,
+ (SCM x),
"Return @code{#t} if the real number @var{x} is a NaN,\n"
"or @code{#f} otherwise.")
#define FUNC_NAME s_scm_nan_p
@@ -640,7 +638,7 @@ SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0,
else if (scm_is_real (x))
return SCM_BOOL_F;
else
- SCM_WRONG_TYPE_ARG (1, x);
+ SCM_WTA_DISPATCH_1 (g_scm_nan_p, x, 1, s_scm_nan_p);
}
#undef FUNC_NAME
@@ -727,7 +725,7 @@ SCM_DEFINE (scm_nan, "nan", 0, 0, 0,
SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
(SCM x),
"Return the absolute value of @var{x}.")
-#define FUNC_NAME
+#define FUNC_NAME s_scm_abs
{
if (SCM_I_INUMP (x))
{
@@ -769,11 +767,10 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
#undef FUNC_NAME
-SCM_GPROC (s_quotient, "quotient", 2, 0, 0, scm_quotient, g_quotient);
-/* "Return the quotient of the numbers @var{x} and @var{y}."
- */
-SCM
-scm_quotient (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_quotient, "quotient", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the quotient of the numbers @var{x} and @var{y}.")
+#define FUNC_NAME s_scm_quotient
{
if (SCM_LIKELY (SCM_I_INUMP (x)))
{
@@ -782,7 +779,7 @@ scm_quotient (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_quotient);
+ scm_num_overflow (s_scm_quotient);
else
{
scm_t_inum z = xx / yy;
@@ -806,7 +803,7 @@ scm_quotient (SCM x, SCM y)
return SCM_INUM0;
}
else
- SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
}
else if (SCM_BIGP (x))
{
@@ -814,7 +811,7 @@ scm_quotient (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_quotient);
+ scm_num_overflow (s_scm_quotient);
else if (SCM_UNLIKELY (yy == 1))
return x;
else
@@ -843,21 +840,21 @@ scm_quotient (SCM x, SCM y)
return scm_i_normbig (result);
}
else
- SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
}
else
- SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG1, s_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient);
}
+#undef FUNC_NAME
-SCM_GPROC (s_remainder, "remainder", 2, 0, 0, scm_remainder, g_remainder);
-/* "Return the remainder of the numbers @var{x} and @var{y}.\n"
- * "@lisp\n"
- * "(remainder 13 4) @result{} 1\n"
- * "(remainder -13 4) @result{} -1\n"
- * "@end lisp"
- */
-SCM
-scm_remainder (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_remainder, "remainder", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the remainder of the numbers @var{x} and @var{y}.\n"
+ "@lisp\n"
+ "(remainder 13 4) @result{} 1\n"
+ "(remainder -13 4) @result{} -1\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_remainder
{
if (SCM_LIKELY (SCM_I_INUMP (x)))
{
@@ -865,7 +862,7 @@ scm_remainder (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_remainder);
+ scm_num_overflow (s_scm_remainder);
else
{
/* C99 specifies that "%" is the remainder corresponding to a
@@ -889,7 +886,7 @@ scm_remainder (SCM x, SCM y)
return x;
}
else
- SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
}
else if (SCM_BIGP (x))
{
@@ -897,7 +894,7 @@ scm_remainder (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_remainder);
+ scm_num_overflow (s_scm_remainder);
else
{
SCM result = scm_i_mkbig ();
@@ -918,22 +915,22 @@ scm_remainder (SCM x, SCM y)
return scm_i_normbig (result);
}
else
- SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
}
else
- SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG1, s_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder);
}
+#undef FUNC_NAME
-SCM_GPROC (s_modulo, "modulo", 2, 0, 0, scm_modulo, g_modulo);
-/* "Return the modulo of the numbers @var{x} and @var{y}.\n"
- * "@lisp\n"
- * "(modulo 13 4) @result{} 1\n"
- * "(modulo -13 4) @result{} 3\n"
- * "@end lisp"
- */
-SCM
-scm_modulo (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the modulo of the numbers @var{x} and @var{y}.\n"
+ "@lisp\n"
+ "(modulo 13 4) @result{} 1\n"
+ "(modulo -13 4) @result{} 3\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_modulo
{
if (SCM_LIKELY (SCM_I_INUMP (x)))
{
@@ -942,7 +939,7 @@ scm_modulo (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_modulo);
+ scm_num_overflow (s_scm_modulo);
else
{
/* C99 specifies that "%" is the remainder corresponding to a
@@ -1008,7 +1005,7 @@ scm_modulo (SCM x, SCM y)
}
}
else
- SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
+ SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
}
else if (SCM_BIGP (x))
{
@@ -1016,7 +1013,7 @@ scm_modulo (SCM x, SCM y)
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_modulo);
+ scm_num_overflow (s_scm_modulo);
else
{
SCM result = scm_i_mkbig ();
@@ -1049,11 +1046,12 @@ scm_modulo (SCM x, SCM y)
return scm_i_normbig (result);
}
else
- SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
+ SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
}
else
- SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo);
+ SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG1, s_scm_modulo);
}
+#undef FUNC_NAME
static SCM scm_i_inexact_euclidean_quotient (double x, double y);
static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y);
@@ -3054,8 +3052,9 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
"Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
"exact integer, @var{n} can be any number.\n"
"\n"
- "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
- "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
+ "Negative @var{k} is supported, and results in\n"
+ "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
+ "@math{@var{n}^0} is 1, as usual, and that\n"
"includes @math{0^0} is 1.\n"
"\n"
"@lisp\n"
@@ -5038,12 +5037,11 @@ scm_geq_p (SCM x, SCM y)
#undef FUNC_NAME
-SCM_GPROC (s_zero_p, "zero?", 1, 0, 0, scm_zero_p, g_zero_p);
-/* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
- * "zero."
- */
-SCM
-scm_zero_p (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_zero_p, "zero?", 1, 0, 0,
+ (SCM z),
+ "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
+ "zero.")
+#define FUNC_NAME s_scm_zero_p
{
if (SCM_I_INUMP (z))
return scm_from_bool (scm_is_eq (z, SCM_INUM0));
@@ -5057,16 +5055,16 @@ scm_zero_p (SCM z)
else if (SCM_FRACTIONP (z))
return SCM_BOOL_F;
else
- SCM_WTA_DISPATCH_1 (g_zero_p, z, SCM_ARG1, s_zero_p);
+ SCM_WTA_DISPATCH_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p);
}
+#undef FUNC_NAME
-SCM_GPROC (s_positive_p, "positive?", 1, 0, 0, scm_positive_p, g_positive_p);
-/* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
- * "zero."
- */
-SCM
-scm_positive_p (SCM x)
+SCM_PRIMITIVE_GENERIC (scm_positive_p, "positive?", 1, 0, 0,
+ (SCM x),
+ "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
+ "zero.")
+#define FUNC_NAME s_scm_positive_p
{
if (SCM_I_INUMP (x))
return scm_from_bool (SCM_I_INUM (x) > 0);
@@ -5081,16 +5079,16 @@ scm_positive_p (SCM x)
else if (SCM_FRACTIONP (x))
return scm_positive_p (SCM_FRACTION_NUMERATOR (x));
else
- SCM_WTA_DISPATCH_1 (g_positive_p, x, SCM_ARG1, s_positive_p);
+ SCM_WTA_DISPATCH_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p);
}
+#undef FUNC_NAME
-SCM_GPROC (s_negative_p, "negative?", 1, 0, 0, scm_negative_p, g_negative_p);
-/* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
- * "zero."
- */
-SCM
-scm_negative_p (SCM x)
+SCM_PRIMITIVE_GENERIC (scm_negative_p, "negative?", 1, 0, 0,
+ (SCM x),
+ "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
+ "zero.")
+#define FUNC_NAME s_scm_negative_p
{
if (SCM_I_INUMP (x))
return scm_from_bool (SCM_I_INUM (x) < 0);
@@ -5105,8 +5103,9 @@ scm_negative_p (SCM x)
else if (SCM_FRACTIONP (x))
return scm_negative_p (SCM_FRACTION_NUMERATOR (x));
else
- SCM_WTA_DISPATCH_1 (g_negative_p, x, SCM_ARG1, s_negative_p);
+ SCM_WTA_DISPATCH_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p);
}
+#undef FUNC_NAME
/* scm_min and scm_max return an inexact when either argument is inexact, as
@@ -6695,9 +6694,9 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
Written by Jerry D. Hedden, (C) FSF.
See the file `COPYING' for terms applying to this program. */
-SCM_DEFINE (scm_expt, "expt", 2, 0, 0,
- (SCM x, SCM y),
- "Return @var{x} raised to the power of @var{y}.")
+SCM_PRIMITIVE_GENERIC (scm_expt, "expt", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return @var{x} raised to the power of @var{y}.")
#define FUNC_NAME s_scm_expt
{
if (scm_is_integer (y))
@@ -6727,8 +6726,12 @@ SCM_DEFINE (scm_expt, "expt", 2, 0, 0,
{
return scm_from_double (pow (scm_to_double (x), scm_to_double (y)));
}
- else
+ else if (scm_is_complex (x) && scm_is_complex (y))
return scm_exp (scm_product (scm_log (x), y));
+ else if (scm_is_complex (x))
+ SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG2, s_scm_expt);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt);
}
#undef FUNC_NAME
@@ -7054,90 +7057,76 @@ SCM_DEFINE (scm_make_polar, "make-polar", 2, 0, 0,
#undef FUNC_NAME
-SCM_GPROC (s_real_part, "real-part", 1, 0, 0, scm_real_part, g_real_part);
-/* "Return the real part of the number @var{z}."
- */
-SCM
-scm_real_part (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_real_part, "real-part", 1, 0, 0,
+ (SCM z),
+ "Return the real part of the number @var{z}.")
+#define FUNC_NAME s_scm_real_part
{
- if (SCM_I_INUMP (z))
- return z;
- else if (SCM_BIGP (z))
- return z;
- else if (SCM_REALP (z))
- return z;
- else if (SCM_COMPLEXP (z))
+ if (SCM_COMPLEXP (z))
return scm_from_double (SCM_COMPLEX_REAL (z));
- else if (SCM_FRACTIONP (z))
+ else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z))
return z;
else
- SCM_WTA_DISPATCH_1 (g_real_part, z, SCM_ARG1, s_real_part);
+ SCM_WTA_DISPATCH_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part);
}
+#undef FUNC_NAME
-SCM_GPROC (s_imag_part, "imag-part", 1, 0, 0, scm_imag_part, g_imag_part);
-/* "Return the imaginary part of the number @var{z}."
- */
-SCM
-scm_imag_part (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_imag_part, "imag-part", 1, 0, 0,
+ (SCM z),
+ "Return the imaginary part of the number @var{z}.")
+#define FUNC_NAME s_scm_imag_part
{
- if (SCM_I_INUMP (z))
- return SCM_INUM0;
- else if (SCM_BIGP (z))
- return SCM_INUM0;
+ if (SCM_COMPLEXP (z))
+ return scm_from_double (SCM_COMPLEX_IMAG (z));
else if (SCM_REALP (z))
return flo0;
- else if (SCM_COMPLEXP (z))
- return scm_from_double (SCM_COMPLEX_IMAG (z));
- else if (SCM_FRACTIONP (z))
+ else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
return SCM_INUM0;
else
- SCM_WTA_DISPATCH_1 (g_imag_part, z, SCM_ARG1, s_imag_part);
+ SCM_WTA_DISPATCH_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part);
}
+#undef FUNC_NAME
-SCM_GPROC (s_numerator, "numerator", 1, 0, 0, scm_numerator, g_numerator);
-/* "Return the numerator of the number @var{z}."
- */
-SCM
-scm_numerator (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_numerator, "numerator", 1, 0, 0,
+ (SCM z),
+ "Return the numerator of the number @var{z}.")
+#define FUNC_NAME s_scm_numerator
{
- if (SCM_I_INUMP (z))
- return z;
- else if (SCM_BIGP (z))
+ if (SCM_I_INUMP (z) || SCM_BIGP (z))
return z;
else if (SCM_FRACTIONP (z))
return SCM_FRACTION_NUMERATOR (z);
else if (SCM_REALP (z))
return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
else
- SCM_WTA_DISPATCH_1 (g_numerator, z, SCM_ARG1, s_numerator);
+ SCM_WTA_DISPATCH_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator);
}
+#undef FUNC_NAME
-SCM_GPROC (s_denominator, "denominator", 1, 0, 0, scm_denominator, g_denominator);
-/* "Return the denominator of the number @var{z}."
- */
-SCM
-scm_denominator (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_denominator, "denominator", 1, 0, 0,
+ (SCM z),
+ "Return the denominator of the number @var{z}.")
+#define FUNC_NAME s_scm_denominator
{
- if (SCM_I_INUMP (z))
- return SCM_INUM1;
- else if (SCM_BIGP (z))
+ if (SCM_I_INUMP (z) || SCM_BIGP (z))
return SCM_INUM1;
else if (SCM_FRACTIONP (z))
return SCM_FRACTION_DENOMINATOR (z);
else if (SCM_REALP (z))
return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
else
- SCM_WTA_DISPATCH_1 (g_denominator, z, SCM_ARG1, s_denominator);
+ SCM_WTA_DISPATCH_1 (g_scm_denominator, z, SCM_ARG1, s_scm_denominator);
}
+#undef FUNC_NAME
-SCM_GPROC (s_magnitude, "magnitude", 1, 0, 0, scm_magnitude, g_magnitude);
-/* "Return the magnitude of the number @var{z}. This is the same as\n"
- * "@code{abs} for real arguments, but also allows complex numbers."
- */
-SCM
-scm_magnitude (SCM z)
+
+SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0,
+ (SCM z),
+ "Return the magnitude of the number @var{z}. This is the same as\n"
+ "@code{abs} for real arguments, but also allows complex numbers.")
+#define FUNC_NAME s_scm_magnitude
{
if (SCM_I_INUMP (z))
{
@@ -7170,15 +7159,15 @@ scm_magnitude (SCM z)
SCM_FRACTION_DENOMINATOR (z));
}
else
- SCM_WTA_DISPATCH_1 (g_magnitude, z, SCM_ARG1, s_magnitude);
+ SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude);
}
+#undef FUNC_NAME
-SCM_GPROC (s_angle, "angle", 1, 0, 0, scm_angle, g_angle);
-/* "Return the angle of the complex number @var{z}."
- */
-SCM
-scm_angle (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0,
+ (SCM z),
+ "Return the angle of the complex number @var{z}.")
+#define FUNC_NAME s_scm_angle
{
/* atan(0,-1) is pi and it'd be possible to have that as a constant like
flo0 to save allocating a new flonum with scm_from_double each time.
@@ -7216,15 +7205,15 @@ scm_angle (SCM z)
else return scm_from_double (atan2 (0.0, -1.0));
}
else
- SCM_WTA_DISPATCH_1 (g_angle, z, SCM_ARG1, s_angle);
+ SCM_WTA_DISPATCH_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
}
+#undef FUNC_NAME
-SCM_GPROC (s_exact_to_inexact, "exact->inexact", 1, 0, 0, scm_exact_to_inexact, g_exact_to_inexact);
-/* Convert the number @var{x} to its inexact representation.\n"
- */
-SCM
-scm_exact_to_inexact (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact, "exact->inexact", 1, 0, 0,
+ (SCM z),
+ "Convert the number @var{z} to its inexact representation.\n")
+#define FUNC_NAME s_scm_exact_to_inexact
{
if (SCM_I_INUMP (z))
return scm_from_double ((double) SCM_I_INUM (z));
@@ -7235,22 +7224,21 @@ scm_exact_to_inexact (SCM z)
else if (SCM_INEXACTP (z))
return z;
else
- SCM_WTA_DISPATCH_1 (g_exact_to_inexact, z, 1, s_exact_to_inexact);
+ SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact, z, 1, s_scm_exact_to_inexact);
}
+#undef FUNC_NAME
-SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
- (SCM z),
- "Return an exact number that is numerically closest to @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
+ (SCM z),
+ "Return an exact number that is numerically closest to @var{z}.")
#define FUNC_NAME s_scm_inexact_to_exact
{
- if (SCM_I_INUMP (z))
- return z;
- else if (SCM_BIGP (z))
+ if (SCM_I_INUMP (z) || SCM_BIGP (z))
return z;
else if (SCM_REALP (z))
{
- if (isinf (SCM_REAL_VALUE (z)) || isnan (SCM_REAL_VALUE (z)))
+ if (!DOUBLE_IS_FINITE (SCM_REAL_VALUE (z)))
SCM_OUT_OF_RANGE (1, z);
else
{
@@ -7272,7 +7260,7 @@ SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
else if (SCM_FRACTIONP (z))
return z;
else
- SCM_WRONG_TYPE_ARG (1, z);
+ SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact, z, 1, s_scm_inexact_to_exact);
}
#undef FUNC_NAME
@@ -7712,9 +7700,9 @@ scm_is_number (SCM z)
real-only case, and because we have to test SCM_COMPLEXP anyway so may as
well use it to go straight to the applicable C func. */
-SCM_DEFINE (scm_log, "log", 1, 0, 0,
- (SCM z),
- "Return the natural logarithm of @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0,
+ (SCM z),
+ "Return the natural logarithm of @var{z}.")
#define FUNC_NAME s_scm_log
{
if (SCM_COMPLEXP (z))
@@ -7728,7 +7716,7 @@ SCM_DEFINE (scm_log, "log", 1, 0, 0,
atan2 (im, re));
#endif
}
- else
+ else if (SCM_NUMBERP (z))
{
/* ENHANCE-ME: When z is a bignum the logarithm will fit a double
although the value itself overflows. */
@@ -7739,13 +7727,15 @@ SCM_DEFINE (scm_log, "log", 1, 0, 0,
else
return scm_c_make_rectangular (l, M_PI);
}
+ else
+ SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
- (SCM z),
- "Return the base 10 logarithm of @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0,
+ (SCM z),
+ "Return the base 10 logarithm of @var{z}.")
#define FUNC_NAME s_scm_log10
{
if (SCM_COMPLEXP (z))
@@ -7763,7 +7753,7 @@ SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
M_LOG10E * atan2 (im, re));
#endif
}
- else
+ else if (SCM_NUMBERP (z))
{
/* ENHANCE-ME: When z is a bignum the logarithm will fit a double
although the value itself overflows. */
@@ -7774,14 +7764,16 @@ SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
else
return scm_c_make_rectangular (l, M_LOG10E * M_PI);
}
+ else
+ SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_exp, "exp", 1, 0, 0,
- (SCM z),
- "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
- "base of natural logarithms (2.71828@dots{}).")
+SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0,
+ (SCM z),
+ "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
+ "base of natural logarithms (2.71828@dots{}).")
#define FUNC_NAME s_scm_exp
{
if (SCM_COMPLEXP (z))
@@ -7793,51 +7785,55 @@ SCM_DEFINE (scm_exp, "exp", 1, 0, 0,
SCM_COMPLEX_IMAG (z));
#endif
}
- else
+ else if (SCM_NUMBERP (z))
{
/* When z is a negative bignum the conversion to double overflows,
giving -infinity, but that's ok, the exp is still 0.0. */
return scm_from_double (exp (scm_to_double (z)));
}
+ else
+ SCM_WTA_DISPATCH_1 (g_scm_exp, z, 1, s_scm_exp);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
- (SCM x),
- "Return the square root of @var{z}. Of the two possible roots\n"
- "(positive and negative), the one with the a positive real part\n"
- "is returned, or if that's zero then a positive imaginary part.\n"
- "Thus,\n"
- "\n"
- "@example\n"
- "(sqrt 9.0) @result{} 3.0\n"
- "(sqrt -9.0) @result{} 0.0+3.0i\n"
- "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
- "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
- "@end example")
+SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
+ (SCM z),
+ "Return the square root of @var{z}. Of the two possible roots\n"
+ "(positive and negative), the one with the a positive real part\n"
+ "is returned, or if that's zero then a positive imaginary part.\n"
+ "Thus,\n"
+ "\n"
+ "@example\n"
+ "(sqrt 9.0) @result{} 3.0\n"
+ "(sqrt -9.0) @result{} 0.0+3.0i\n"
+ "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
+ "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
+ "@end example")
#define FUNC_NAME s_scm_sqrt
{
- if (SCM_COMPLEXP (x))
+ if (SCM_COMPLEXP (z))
{
#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
&& defined SCM_COMPLEX_VALUE
- return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x)));
+ return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z)));
#else
- double re = SCM_COMPLEX_REAL (x);
- double im = SCM_COMPLEX_IMAG (x);
+ double re = SCM_COMPLEX_REAL (z);
+ double im = SCM_COMPLEX_IMAG (z);
return scm_c_make_polar (sqrt (hypot (re, im)),
0.5 * atan2 (im, re));
#endif
}
- else
+ else if (SCM_NUMBERP (z))
{
- double xx = scm_to_double (x);
+ double xx = scm_to_double (z);
if (xx < 0)
return scm_c_make_rectangular (0.0, sqrt (-xx));
else
return scm_from_double (sqrt (xx));
}
+ else
+ SCM_WTA_DISPATCH_1 (g_scm_sqrt, z, 1, s_scm_sqrt);
}
#undef FUNC_NAME
diff --git a/libguile/numbers.h b/libguile/numbers.h
index 76d2972..2cf3fd7 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -169,8 +169,9 @@ typedef struct scm_t_complex
SCM_API SCM scm_exact_p (SCM x);
SCM_API SCM scm_odd_p (SCM n);
SCM_API SCM scm_even_p (SCM n);
-SCM_API SCM scm_inf_p (SCM n);
-SCM_API SCM scm_nan_p (SCM n);
+SCM_API SCM scm_finite_p (SCM x);
+SCM_API SCM scm_inf_p (SCM x);
+SCM_API SCM scm_nan_p (SCM x);
SCM_API SCM scm_inf (void);
SCM_API SCM scm_nan (void);
SCM_API SCM scm_abs (SCM x);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 9cf9202..01bccda 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -281,8 +281,7 @@
;;;
(with-test-prefix "exp"
- (pass-if "documented?"
- (documented? exp))
+ (pass-if (documented? exp))
(pass-if-exception "no args" exception:wrong-num-args
(exp))
@@ -426,9 +425,7 @@
;;;
(with-test-prefix "quotient"
-
- (expect-fail "documented?"
- (documented? quotient))
+ (pass-if (documented? quotient))
(with-test-prefix "0 / n"
@@ -642,9 +639,7 @@
;;;
(with-test-prefix "remainder"
-
- (expect-fail "documented?"
- (documented? remainder))
+ (pass-if (documented? remainder))
(with-test-prefix "0 / n"
@@ -837,9 +832,7 @@
;;;
(with-test-prefix "modulo"
-
- (expect-fail "documented?"
- (documented? modulo))
+ (pass-if (documented? modulo))
(with-test-prefix "0 % n"
@@ -2354,7 +2347,7 @@
;;;
(with-test-prefix "zero?"
- (expect-fail (documented? zero?))
+ (pass-if (documented? zero?))
(pass-if (zero? 0))
(pass-if (not (zero? 7)))
(pass-if (not (zero? -7)))
@@ -2368,7 +2361,7 @@
;;;
(with-test-prefix "positive?"
- (expect-fail (documented? positive?))
+ (pass-if (documented? positive?))
(pass-if (positive? 1))
(pass-if (positive? (+ fixnum-max 1)))
(pass-if (positive? 1.3))
@@ -2382,7 +2375,7 @@
;;;
(with-test-prefix "negative?"
- (expect-fail (documented? negative?))
+ (pass-if (documented? negative?))
(pass-if (not (negative? 1)))
(pass-if (not (negative? (+ fixnum-max 1))))
(pass-if (not (negative? 1.3)))
@@ -3118,6 +3111,7 @@
;;;
(with-test-prefix "expt"
+ (pass-if (documented? expt))
(pass-if-exception "non-numeric base" exception:wrong-type-arg
(expt #t 0))
(pass-if (eqv? 1 (expt 0 0)))
@@ -3199,15 +3193,32 @@
;;; real-part
;;;
+(with-test-prefix "real-part"
+ (pass-if (documented? real-part))
+ (pass-if (eqv? 5.0 (real-part 5.0)))
+ (pass-if (eqv? 0.0 (real-part +5.0i)))
+ (pass-if (eqv? 5 (real-part 5)))
+ (pass-if (eqv? 1/5 (real-part 1/5)))
+ (pass-if (eqv? (1+ fixnum-max) (real-part (1+ fixnum-max)))))
+
;;;
;;; imag-part
;;;
+(with-test-prefix "imag-part"
+ (pass-if (documented? imag-part))
+ (pass-if (eqv? 0.0 (imag-part 5.0)))
+ (pass-if (eqv? 5.0 (imag-part +5.0i)))
+ (pass-if (eqv? 0 (imag-part 5)))
+ (pass-if (eqv? 0 (imag-part 1/5)))
+ (pass-if (eqv? 0 (imag-part (1+ fixnum-max)))))
+
;;;
;;; magnitude
;;;
(with-test-prefix "magnitude"
+ (pass-if (documented? magnitude))
(pass-if (= 0 (magnitude 0)))
(pass-if (= 1 (magnitude 1)))
(pass-if (= 1 (magnitude -1)))
@@ -3227,6 +3238,8 @@
(define (almost= x y)
(> 0.01 (magnitude (- x y))))
+ (pass-if (documented? angle))
+
(pass-if "inum +ve" (= 0 (angle 1)))
(pass-if "inum -ve" (almost= pi (angle -1)))
@@ -3241,7 +3254,8 @@
;;;
(with-test-prefix "inexact->exact"
-
+ (pass-if (documented? inexact->exact))
+
(pass-if-exception "+inf" exception:out-of-range
(inexact->exact +inf.0))
@@ -3263,6 +3277,7 @@
;;;
(with-test-prefix "integer-expt"
+ (pass-if (documented? integer-expt))
(pass-if-exception "non-numeric base" exception:wrong-type-arg
(integer-expt #t 0))
@@ -3294,6 +3309,7 @@
;;;
(with-test-prefix "integer-length"
+ (pass-if (documented? integer-length))
(with-test-prefix "-2^i, ...11100..00"
(do ((n -1 (ash n 1))
@@ -3321,8 +3337,7 @@
;;;
(with-test-prefix "log"
- (pass-if "documented?"
- (documented? log))
+ (pass-if (documented? log))
(pass-if-exception "no args" exception:wrong-num-args
(log))
@@ -3349,8 +3364,7 @@
;;;
(with-test-prefix "log10"
- (pass-if "documented?"
- (documented? log10))
+ (pass-if (documented? log10))
(pass-if-exception "no args" exception:wrong-num-args
(log10))
@@ -3377,6 +3391,8 @@
;;;
(with-test-prefix "logbit?"
+ (pass-if (documented? logbit?))
+
(pass-if (eq? #f (logbit? 0 0)))
(pass-if (eq? #f (logbit? 1 0)))
(pass-if (eq? #f (logbit? 31 0)))
@@ -3412,6 +3428,7 @@
;;;
(with-test-prefix "logcount"
+ (pass-if (documented? logcount))
(with-test-prefix "-2^i, meaning ...11100..00"
(do ((n -1 (ash n 1))
@@ -3439,6 +3456,8 @@
;;;
(with-test-prefix "logior"
+ (pass-if (documented? logior))
+
(pass-if (eqv? -1 (logior (ash -1 1) 1)))
;; check that bignum or bignum+inum args will reduce to an inum
@@ -3468,6 +3487,8 @@
;;;
(with-test-prefix "lognot"
+ (pass-if (documented? lognot))
+
(pass-if (= -1 (lognot 0)))
(pass-if (= 0 (lognot -1)))
(pass-if (= -2 (lognot 1)))
@@ -3483,8 +3504,7 @@
;;;
(with-test-prefix "sqrt"
- (pass-if "documented?"
- (documented? sqrt))
+ (pass-if (documented? sqrt))
(pass-if-exception "no args" exception:wrong-num-args
(sqrt))
@@ -3626,6 +3646,13 @@
test-numerators))
test-denominators))
+ (pass-if (documented? euclidean/))
+ (pass-if (documented? euclidean-quotient))
+ (pass-if (documented? euclidean-remainder))
+ (pass-if (documented? centered/))
+ (pass-if (documented? centered-quotient))
+ (pass-if (documented? centered-remainder))
+
(with-test-prefix "euclidean-quotient"
(do-tests-1 'euclidean-quotient
euclidean-quotient
--
1.5.6.5
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #5: Improve extensibility of `expt' and `integer-expt' --]
[-- Type: text/x-diff, Size: 6282 bytes --]
From 37dd087091b3d7305a0546d5d35659249967c8cc Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 30 Jan 2011 10:51:36 -0500
Subject: [PATCH] Improve extensibility of `expt' and `integer-expt'
* libguile/numbers.c (scm_integer_expt): No longer require that the
first argument be a number, in order to improve extensibility. This
allows us to efficiently raise arbitrary objects to an integer power
as long as we can multiply those objects. For example, this allows us
to efficiently exponentiate matrices if we define only multiplication
methods for matrices. Note also that scm_expt calls this procedure
whenever the exponent is an integer, regardless of the type of the
first argument. Also rearrange the order in which we test special
cases.
* test-suite/tests/numbers.test (expt, integer-expt): Comment out tests
that required `(expt #t 0)' and `(integer-expt #t 0)' to throw
exceptions. Add tests for (expt #t 2) and `(integer-expt #t 2)
instead.
* NEWS: Add NEWS entry
---
NEWS | 11 +++++++++++
libguile/numbers.c | 28 +++++++++++++++++-----------
test-suite/tests/numbers.test | 32 ++++++++++++++++++++++++++++++--
3 files changed, 58 insertions(+), 13 deletions(-)
diff --git a/NEWS b/NEWS
index 33e49a4..b78e1d1 100644
--- a/NEWS
+++ b/NEWS
@@ -58,6 +58,17 @@ integer-expt. This is more correct, and conforming to R6RS, but seems
to be incompatible with R5RS, which would return 0 for all non-zero
values of N.
+*** `expt' and `integer-expt' are more generic, less strict
+
+When raising to an exact non-negative integer exponent, `expt' and
+`integer-expt' are now able to exponentiate any object that can be
+multiplied using `*'. They can also raise an object to an exact
+negative integer power if its reciprocal can be taken using `/'.
+In order to allow this, the type of the first argument is no longer
+checked when raising to an exact integer power. If the exponent is 0
+or 1, the first parameter is not manipulated at all, and need not
+even support multiplication.
+
*** Infinities are no longer integers, nor rationals
scm_integer_p `integer?' and scm_rational_p `rational?' now return #f
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 955404d..0dc0ed1 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -3070,23 +3070,29 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
int i2_is_big = 0;
SCM acc = SCM_I_MAKINUM (1L);
- SCM_VALIDATE_NUMBER (SCM_ARG1, n);
- if (!SCM_I_INUMP (k) && !SCM_BIGP (k))
+ /* Specifically refrain from checking the type of the first argument.
+ This allows us to exponentiate any object that can be multiplied.
+ If we must raise to a negative power, we must also be able to
+ take its reciprocal. */
+ if (!SCM_LIKELY (SCM_I_INUMP (k)) && !SCM_LIKELY (SCM_BIGP (k)))
SCM_WRONG_TYPE_ARG (2, k);
- if (scm_is_true (scm_zero_p (n)))
- {
- if (scm_is_true (scm_zero_p (k))) /* 0^0 == 1 per R5RS */
- return acc; /* return exact 1, regardless of n */
- else if (scm_is_true (scm_positive_p (k)))
+ if (SCM_UNLIKELY (scm_is_eq (k, SCM_INUM0)))
+ return SCM_INUM1; /* n^(exact0) is exact 1, regardless of n */
+ else if (SCM_UNLIKELY (scm_is_eq (n, SCM_I_MAKINUM (-1L))))
+ return scm_is_false (scm_even_p (k)) ? n : SCM_INUM1;
+ /* The next check is necessary only because R6RS specifies different
+ behavior for 0^(-k) than for (/ 0). If n is not a scheme number,
+ we simply skip this case and move on. */
+ else if (SCM_NUMBERP (n) && scm_is_true (scm_zero_p (n)))
+ {
+ /* k cannot be 0 at this point, because we
+ have already checked for that case above */
+ if (scm_is_true (scm_positive_p (k)))
return n;
else /* return NaN for (0 ^ k) for negative k per R6RS */
return scm_nan ();
}
- else if (scm_is_eq (n, acc))
- return acc;
- else if (scm_is_eq (n, SCM_I_MAKINUM (-1L)))
- return scm_is_false (scm_even_p (k)) ? n : acc;
if (SCM_I_INUMP (k))
i2 = SCM_I_INUM (k);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 01bccda..8fa7605 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -3112,8 +3112,23 @@
(with-test-prefix "expt"
(pass-if (documented? expt))
+
+ ;;
+ ;; integer-expt no longer requires its first argument to be a scheme
+ ;; number, for the sake of extensibility, and expt calls integer-expt
+ ;; for integer powers. To raise to a positive power, all that is
+ ;; required is that it can be multiplied using `*'. For negative
+ ;; powers we must also be able to find the reciprocal. If we try to
+ ;; raise #t to any power other than 0 or 1 it will throw an exception.
+ ;; However, when raising to the 0 or 1 power, the first argument is
+ ;; not manipulated at all.
+ ;;
+ ;; (pass-if-exception "non-numeric base" exception:wrong-type-arg
+ ;; (expt #t 0))
+ ;;
(pass-if-exception "non-numeric base" exception:wrong-type-arg
- (expt #t 0))
+ (expt #t 2))
+
(pass-if (eqv? 1 (expt 0 0)))
(pass-if (eqv? 1 (expt 0.0 0)))
(pass-if (eqv? 1.0 (expt 0 0.0)))
@@ -3279,8 +3294,21 @@
(with-test-prefix "integer-expt"
(pass-if (documented? integer-expt))
+ ;;
+ ;; integer-expt no longer requires its first argument to be a scheme
+ ;; number, for the sake of extensibility. To raise to a positive
+ ;; power, all that is required is that it can be multiplied using `*'.
+ ;; For negative powers we must also be able to find the reciprocal.
+ ;; If we try to raise #t to any power other than 0 or 1 it will throw
+ ;; an exception, because `*' will fail. However, when raising to the
+ ;; 0 or 1 power, the first argument is not manipulated at all.
+ ;;
+ ;; (pass-if-exception "non-numeric base" exception:wrong-type-arg
+ ;; (integer-expt #t 0))
+ ;;
(pass-if-exception "non-numeric base" exception:wrong-type-arg
- (integer-expt #t 0))
+ (integer-expt #t 2))
+
(pass-if-exception "2^+inf" exception:wrong-type-arg
(integer-expt 2 +inf.0))
(pass-if-exception "2^-inf" exception:wrong-type-arg
--
1.5.6.5
[-- Attachment #6: Simple demonstration of `expt' extensibility --]
[-- Type: text/plain, Size: 3168 bytes --]
(use-modules (oop goops))
(use-modules (srfi srfi-1))
;; A matrix has only one slot, containing a list
;; of rows, where each row is a list of elements.
(define-class <matrix> ()
(rows #:getter matrix-rows #:init-keyword #:rows))
(define (matrix rows)
(make <matrix> #:rows rows))
(define (matrix-ref a i j)
(list-ref (list-ref (matrix-rows a)
i)
j))
(define (display-matrix a port)
(newline port)
(for-each (lambda (row)
(display row port)
(newline port))
(matrix-rows a)))
(define (matrix-cols a)
(apply map list (matrix-rows a)))
(define (matrix-map f a)
(matrix (map (lambda (row) (map f row))
(matrix-rows a))))
(define (scale-matrix s a)
(cond ((= s 1) a)
(else (matrix-map (lambda (elem) (* s elem)) a))))
(define (zero-matrix n m)
(matrix (make-list n (make-list m 0))))
(define (identity-matrix n)
(matrix
(list-tabulate n (lambda (i)
(append (make-list i 0)
(list 1)
(make-list (- n i 1) 0))))))
(define (row-vect . elems) (matrix (list elems)))
(define (col-vect . elems) (matrix (map list elems)))
(define (diagonal-matrix diag)
(let ((n (length diag)))
(matrix (map (lambda (i x)
(append (make-list i 0)
(list x)
(make-list (- n i 1) 0)))
(iota n) diag))))
(define (matrix-mult a b)
(let ((a-rows (matrix-rows a))
(b-cols (matrix-cols b)))
(if (not (= (length (first a-rows))
(length (first b-cols))))
(throw 'matrix-mult a b))
(matrix (map (lambda (a-row)
(map (lambda (b-col)
(apply + (map * a-row b-col)))
b-cols))
a-rows))))
(define-method (display (a <matrix>) port) (display-matrix a port))
(define-method (write (a <matrix>) port) (display-matrix a port))
(define-method (* (s <number>) (a <matrix>)) (scale-matrix s a))
(define-method (* (a <matrix>) (s <number>)) (scale-matrix s a))
(define-method (* (a <matrix>) (b <matrix>)) (matrix-mult a b))
;;
;; Fast fibonacci using matrix exponentiation
;;
;; [0 1] [a] [ b ]
;; [1 1] [b] = [a+b]
;;
;; As shown above, the matrix above transforms a column vector (a b)
;; into (b a+b). This column vector is the state of an iterative
;; algorithm to compute fibonacci, and the matrix is the
;; transformation done at each step.
;;
;; However, instead of multiplying the state vector times the matrix
;; at each step, we exponentiate the matrix to a power, where the
;; expononent N equals the number of steps we want to take. This
;; yields a matrix that advances the state vector by N steps. The
;; trick is to use the same fast exponentiation algorithm used on
;; integers. The number of matrix multiplications performed is
;; proportional to the logarithm of the exponent.
;;
;; With the help of GOOPS, we use expt to exponentiate the matrix. We
;; don't even have to define an expt method for matrices. All that we
;; need is multiplication. expt calls integer-expt, since the
;; exponent is an integer. Thus this depends on integer-expt
;; accepting a base that is not a number.
;;
(define (fib n)
(let* ((m1 (matrix '((0 1)
(1 1))))
(m (expt m1 n))
(v (* m (col-vect 0 1))))
(matrix-ref v 0 0)))
^ permalink raw reply related [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-30 16:27 [PATCH] Fast R6RS div/mod; improved extensibility of numerics Mark H Weaver
@ 2011-01-30 22:24 ` Andy Wingo
2011-01-30 22:55 ` Ludovic Courtès
` (2 more replies)
0 siblings, 3 replies; 12+ messages in thread
From: Andy Wingo @ 2011-01-30 22:24 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
On Sun 30 Jan 2011 17:27, Mark H Weaver <mhw@netris.org> writes:
> Here are another four patches of numerics improvements.
Applied all up to the last. Sweet! And thank you!
> The other three I named `centered/', `centered-quotient', and
> `centered-remainder', because the range of their remainder operation is
> centered on 0 in a twos-complement sort of way. They are almost the
> same as Taylor's `round/' et al, but in my opinion they have superior
> properties, as I explained in a recent commentary on this list.
I agree they do have some good properties, and am happy to see them in
the core, and with non-privileged names (allowing for round_quotient et
al to be added, possibly, in the future).
If you have time though, I wonder would scm_centered_divide be a better
name than scm_centered_quo_and_rem? "Divide" is the pronunciation of
the "/" procedure, I think.
> The last two patches improve the extensibility of numeric operations,
> and as a side effect add documentation strings to many of them.
Much appreciated.
> The last patch might be slightly controversial. Although it does not
> make `integer-expt' a generic function, nonetheless it can now
> exponentiate _any_ scheme object that can be multiplied using `*'.
To me this is fine, in principle; but I wonder about causing `(expt #t
0)' to produce 1. It seems that in any case the range of the output
should be continuous -- i.e. if (expt <my-matrix> 1) produces a matrix,
and (expt <my-matrix> n) for n > 1 makes a matrix, then (expt
<my-matrix> 0) should also produce a matrix. Or no?
Cheers,
Andy
--
http://wingolog.org/
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-30 22:24 ` Andy Wingo
@ 2011-01-30 22:55 ` Ludovic Courtès
2011-01-31 6:19 ` [PATCH] Rework the testing framework for number-theoretic division operators Mark H Weaver
2011-01-31 17:14 ` [PATCH] Fast R6RS div/mod; improved extensibility of numerics Mark H Weaver
2011-01-31 17:35 ` Mark H Weaver
2 siblings, 1 reply; 12+ messages in thread
From: Ludovic Courtès @ 2011-01-30 22:55 UTC (permalink / raw)
To: guile-devel
Hello!
Thanks for working on this!
There are related test failures on i686-linux-gnu:
http://hydra.nixos.org/build/875025
You can see the log and browse the build tree from there.
(I don’t want to be offensive, it’s just that I just feel a bit of
pressure with this release thing. ;-))
Ludo’.
^ permalink raw reply [flat|nested] 12+ messages in thread
* [PATCH] Rework the testing framework for number-theoretic division operators
2011-01-30 22:55 ` Ludovic Courtès
@ 2011-01-31 6:19 ` Mark H Weaver
2011-01-31 8:52 ` Andy Wingo
0 siblings, 1 reply; 12+ messages in thread
From: Mark H Weaver @ 2011-01-31 6:19 UTC (permalink / raw)
To: Ludovic Courtès; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 221 bytes --]
ludo@gnu.org (Ludovic Courtès) writes:
> There are related test failures on i686-linux-gnu:
>
> http://hydra.nixos.org/build/875025
I believe this patch should fix the test failures.
Best,
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Rework the testing framework for number-theoretic division operators --]
[-- Type: text/x-diff, Size: 11848 bytes --]
From 650df0667f16cbbf3c05c8774d4d08c793605f47 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 31 Jan 2011 00:42:35 -0500
Subject: [PATCH] Rework the testing framework for number-theoretic division operators
* test-suite/tests/numbers.test (test-eqv?): Remove special handling of
zeroes. Zeroes are now compared like all other numbers. Exact
numbers are compared with `eqv?' and inexact numbers are compared to
within test-epsilon.
Rework the testing framework for number-theoretic division operators:
`euclidean/', `euclidean-quotient', `euclidean-remainder',
`centered/', `centered-quotient', and `centered-remainder'.
Previously we compared all test results against a simple scheme
implementation of the same operations. However, these operations have
discontinuous jumps where a tiny change in the inputs can lead to a
large change in the outputs, e.g.:
(euclidean/ 130.00000000000 10/7) ==> 91.0 and 0.0
(euclidean/ 129.99999999999 10/7) ==> 90.0 and 1.42857142856141
In the new testing scheme, we compare values against the simple
implementations only if the input arguments contain an infinity or a
NaN. In the common case of two finite arguments, we simply make sure
that the outputs of all three operators (e.g. `euclidean/',
`euclidean-quotient', `euclidean-remainder') equal each other, that
outputs are exact iff both inputs are exact, and that the required
properties of the operator are met: that Q is an integer, that R is
within the specified range, and that N = Q*D + R.
---
test-suite/tests/numbers.test | 194 +++++++++++++++++++----------------------
1 files changed, 91 insertions(+), 103 deletions(-)
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 01bccda..0d711b0 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -18,6 +18,7 @@
(define-module (test-suite test-numbers)
#:use-module (test-suite lib)
#:use-module (ice-9 documentation)
+ #:use-module (srfi srfi-1) ; list library
#:use-module (srfi srfi-11)) ; let-values
;;;
@@ -100,12 +101,9 @@
;;
;; 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.
+;; test-epsilon (1e-10) to be considered equal. 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)
@@ -118,7 +116,7 @@
;; Auxiliary predicate used by test-eqv?
(define (test-real-eqv? x y)
- (cond ((or (exact? x) (zero? x) (nan? x) (inf? x))
+ (cond ((or (exact? x) (nan? x) (inf? x))
(eqv? x y))
(else (and (inexact? y) (> test-epsilon (abs (- x y)))))))
@@ -3551,6 +3549,24 @@
(hi (+ hi test-epsilon)))
(<= lo x hi))))
+ ;; (cartesian-product-map list '(a b) '(1 2))
+ ;; ==> ((a 1) (a 2) (b 1) (b 2))
+ (define (cartesian-product-map f . lsts)
+ (define (cartmap rev-head lsts)
+ (if (null? lsts)
+ (list (apply f (reverse rev-head)))
+ (append-map (lambda (x) (cartmap (cons x rev-head) (cdr lsts)))
+ (car lsts))))
+ (cartmap '() lsts))
+
+ (define (cartesian-product-for-each f . lsts)
+ (define (cartfor rev-head lsts)
+ (if (null? lsts)
+ (apply f (reverse rev-head))
+ (for-each (lambda (x) (cartfor (cons x rev-head) (cdr lsts)))
+ (car lsts))))
+ (cartfor '() lsts))
+
(define (safe-euclidean-quotient x y)
(cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
((zero? y) (throw 'divide-by-zero))
@@ -3560,20 +3576,19 @@
(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))))
+ (let ((q (safe-euclidean-quotient x y)))
+ (- x (* y q))))
+
+ (define (valid-euclidean-answer? x y q r)
+ (if (and (finite? x) (finite? y))
+ (and (eq? (exact? q)
+ (exact? r)
+ (and (exact? x) (exact? y)))
+ (integer? q)
+ (test-eqv? r (- x (* q y)))
+ (test-within-range? 0 (abs y) r))
+ (and (test-eqv? q (safe-euclidean-quotient x y))
+ (test-eqv? r (safe-euclidean-remainder x y)))))
(define (safe-centered-quotient x y)
(cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
@@ -3584,37 +3599,36 @@
(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))))
+ (let ((q (safe-centered-quotient x y)))
+ (- x (* y q))))
+
+ (define (valid-centered-answer? x y q r)
+ (if (and (finite? x) (finite? y))
+ (and (eq? (exact? q)
+ (exact? r)
+ (and (exact? x) (exact? y)))
+ (integer? q)
+ (test-eqv? r (- x (* q y)))
+ (test-within-range? (* -1/2 (abs y))
+ (* +1/2 (abs y))
+ r))
+ (and (test-eqv? q (safe-centered-quotient x y))
+ (test-eqv? r (safe-centered-remainder x y)))))
(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)))))
+ (append (cartesian-product-map * '(1 -1)
+ '(123 125 127 130 3 5 10
+ 123.2 125.0 127.2 130.0
+ 123/7 125/7 127/7 130/7))
+ (cartesian-product-map * '(1 -1)
+ '(123 125 127 130 3 5 10)
+ (list 1
+ (+ 1 most-positive-fixnum)
+ (+ 2 most-positive-fixnum)))
+ (list 0 +0.0 -0.0 +inf.0 -inf.0 +nan.0
+ most-negative-fixnum
+ (1+ most-positive-fixnum)
+ (1- most-negative-fixnum))))
(define test-denominators
(list 10 5 10/7 127/2 10.0 63.5
@@ -3623,58 +3637,32 @@
(+ 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))
-
- (pass-if (documented? euclidean/))
- (pass-if (documented? euclidean-quotient))
- (pass-if (documented? euclidean-remainder))
- (pass-if (documented? centered/))
- (pass-if (documented? centered-quotient))
- (pass-if (documented? centered-remainder))
-
- (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))
+ (pass-if (documented? euclidean/))
+ (pass-if (documented? euclidean-quotient))
+ (pass-if (documented? euclidean-remainder))
+
+ (cartesian-product-for-each
+ (lambda (n d)
+ (run-test (list 'euclidean/ n d) #t
+ (lambda ()
+ (let-values (((q r) (euclidean/ n d)))
+ (and (test-eqv? q (euclidean-quotient n d))
+ (test-eqv? r (euclidean-remainder n d))
+ (valid-euclidean-answer? n d q r))))))
+ test-numerators test-denominators))
+
(with-test-prefix "centered/"
- (do-tests-2 'centered/
- centered/
- safe-centered/)))
+ (pass-if (documented? centered/))
+ (pass-if (documented? centered-quotient))
+ (pass-if (documented? centered-remainder))
+
+ (cartesian-product-for-each
+ (lambda (n d)
+ (run-test (list 'centered/ n d) #t
+ (lambda ()
+ (let-values (((q r) (centered/ n d)))
+ (and (test-eqv? q (centered-quotient n d))
+ (test-eqv? r (centered-remainder n d))
+ (valid-centered-answer? n d q r))))))
+ test-numerators test-denominators)))
--
1.5.6.5
^ permalink raw reply related [flat|nested] 12+ messages in thread
* Re: [PATCH] Rework the testing framework for number-theoretic division operators
2011-01-31 6:19 ` [PATCH] Rework the testing framework for number-theoretic division operators Mark H Weaver
@ 2011-01-31 8:52 ` Andy Wingo
0 siblings, 0 replies; 12+ messages in thread
From: Andy Wingo @ 2011-01-31 8:52 UTC (permalink / raw)
To: Mark H Weaver; +Cc: Ludovic Courtès, guile-devel
On Mon 31 Jan 2011 07:19, Mark H Weaver <mhw@netris.org> writes:
> ludo@gnu.org (Ludovic Courtès) writes:
>> There are related test failures on i686-linux-gnu:
>>
>> http://hydra.nixos.org/build/875025
>
> I believe this patch should fix the test failures.
Applied, thanks.
Andy
--
http://wingolog.org/
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-30 22:24 ` Andy Wingo
2011-01-30 22:55 ` Ludovic Courtès
@ 2011-01-31 17:14 ` Mark H Weaver
2011-01-31 17:35 ` Mark H Weaver
2 siblings, 0 replies; 12+ messages in thread
From: Mark H Weaver @ 2011-01-31 17:14 UTC (permalink / raw)
To: Andy Wingo; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 252 bytes --]
Andy Wingo <wingo@pobox.com> writes:
> If you have time though, I wonder would scm_centered_divide be a better
> name than scm_centered_quo_and_rem? "Divide" is the pronunciation of
> the "/" procedure, I think.
Agreed, here's a patch.
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Rename {euclidean,centered}_quo_rem to {euclidean,centered}_divide --]
[-- Type: text/x-diff, Size: 15157 bytes --]
From b5441b4fc27ee4c9686c69307369d24f7cacd660 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 31 Jan 2011 12:03:02 -0500
Subject: [PATCH] Rename {euclidean,centered}_quo_rem to {euclidean,centered}_divide
* libguile/numbers.c (euclidean_quo_rem): Rename to euclidean_divide.
(centered_quo_rem): Rename to {euclidean,centered}_divide.
* libguile/numbers.h: Rename euclidean_quo_rem to euclidean_divide and
centered_quo_rem to centered_divide.
* doc/ref/api-data.texi: Rename euclidean_quo_rem to euclidean_divide and
centered_quo_rem to centered_divide.
---
doc/ref/api-data.texi | 4 +-
libguile/numbers.c | 126 ++++++++++++++++++++++++------------------------
libguile/numbers.h | 4 +-
3 files changed, 67 insertions(+), 67 deletions(-)
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index b090782..b819fcb 100755
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1244,7 +1244,7 @@ values.
@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_divide (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
@@ -1275,7 +1275,7 @@ Note that these operators are equivalent to the R6RS operators
@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_divide (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
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 3a2244f..41d178b 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -1384,10 +1384,10 @@ scm_i_slow_exact_euclidean_remainder (SCM x, SCM y)
}
-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);
+static SCM scm_i_inexact_euclidean_divide (double x, double y);
+static SCM scm_i_slow_exact_euclidean_divide (SCM x, SCM y);
-SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "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"
@@ -1400,7 +1400,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
"(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
+#define FUNC_NAME s_scm_euclidean_divide
{
if (SCM_LIKELY (SCM_I_INUMP (x)))
{
@@ -1408,7 +1408,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_quo_and_rem);
+ scm_num_overflow (s_scm_euclidean_divide);
else
{
scm_t_inum xx = SCM_I_INUM (x);
@@ -1448,13 +1448,13 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
}
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_quo_and_rem
+ return scm_i_inexact_euclidean_divide
(SCM_I_INUM (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+ return scm_i_slow_exact_euclidean_divide (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
- s_scm_euclidean_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+ s_scm_euclidean_divide);
}
else if (SCM_BIGP (x))
{
@@ -1462,7 +1462,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_quo_and_rem);
+ scm_num_overflow (s_scm_euclidean_divide);
else
{
SCM q = scm_i_mkbig ();
@@ -1496,40 +1496,40 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0,
scm_i_normbig (r)));
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_quo_and_rem
+ return scm_i_inexact_euclidean_divide
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+ return scm_i_slow_exact_euclidean_divide (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2,
- s_scm_euclidean_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+ s_scm_euclidean_divide);
}
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
+ return scm_i_inexact_euclidean_divide
(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);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+ s_scm_euclidean_divide);
}
else if (SCM_FRACTIONP (x))
{
if (SCM_REALP (y))
- return scm_i_inexact_euclidean_quo_and_rem
+ return scm_i_inexact_euclidean_divide
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
else
- return scm_i_slow_exact_euclidean_quo_and_rem (x, y);
+ return scm_i_slow_exact_euclidean_divide (x, y);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG1,
- s_scm_euclidean_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
+ s_scm_euclidean_divide);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_euclidean_quo_and_rem (double x, double y)
+scm_i_inexact_euclidean_divide (double x, double y)
{
double q, r;
@@ -1538,7 +1538,7 @@ scm_i_inexact_euclidean_quo_and_rem (double x, double 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? */
+ scm_num_overflow (s_scm_euclidean_divide); /* or return a NaN? */
else
q = guile_NaN;
r = x - q * y;
@@ -1550,22 +1550,22 @@ scm_i_inexact_euclidean_quo_and_rem (double x, double y)
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_i_slow_exact_euclidean_divide (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);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
+ s_scm_euclidean_divide);
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);
+ SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+ s_scm_euclidean_divide);
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);
+ scm_num_overflow (s_scm_euclidean_divide);
r = scm_difference (x, scm_product (q, y));
return scm_values (scm_list_2 (q, r));
}
@@ -2025,11 +2025,11 @@ scm_i_slow_exact_centered_remainder (SCM x, SCM y)
}
-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);
+static SCM scm_i_inexact_centered_divide (double x, double y);
+static SCM scm_i_bigint_centered_divide (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_divide (SCM x, SCM y);
-SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_centered_divide, "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"
@@ -2042,7 +2042,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
"(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
+#define FUNC_NAME s_scm_centered_divide
{
if (SCM_LIKELY (SCM_I_INUMP (x)))
{
@@ -2050,7 +2050,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_centered_quo_and_rem);
+ scm_num_overflow (s_scm_centered_divide);
else
{
scm_t_inum xx = SCM_I_INUM (x);
@@ -2089,18 +2089,18 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
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
+ can fit in a fixnum) to scm_i_bigint_centered_divide */
+ return scm_i_bigint_centered_divide
(scm_i_long2big (SCM_I_INUM (x)), y);
}
else if (SCM_REALP (y))
- return scm_i_inexact_centered_quo_and_rem
+ return scm_i_inexact_centered_divide
(SCM_I_INUM (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_quo_and_rem (x, y);
+ return scm_i_slow_exact_centered_divide (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
- s_scm_centered_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide);
}
else if (SCM_BIGP (x))
{
@@ -2108,7 +2108,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_centered_quo_and_rem);
+ scm_num_overflow (s_scm_centered_divide);
else
{
SCM q = scm_i_mkbig ();
@@ -2146,42 +2146,42 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0,
}
}
else if (SCM_BIGP (y))
- return scm_i_bigint_centered_quo_and_rem (x, y);
+ return scm_i_bigint_centered_divide (x, y);
else if (SCM_REALP (y))
- return scm_i_inexact_centered_quo_and_rem
+ return scm_i_inexact_centered_divide
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_quo_and_rem (x, y);
+ return scm_i_slow_exact_centered_divide (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2,
- s_scm_centered_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide);
}
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
+ return scm_i_inexact_centered_divide
(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);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide);
}
else if (SCM_FRACTIONP (x))
{
if (SCM_REALP (y))
- return scm_i_inexact_centered_quo_and_rem
+ return scm_i_inexact_centered_divide
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
else
- return scm_i_slow_exact_centered_quo_and_rem (x, y);
+ return scm_i_slow_exact_centered_divide (x, y);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG1,
- s_scm_centered_quo_and_rem);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
+ s_scm_centered_divide);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_centered_quo_and_rem (double x, double y)
+scm_i_inexact_centered_divide (double x, double y)
{
double q, r;
@@ -2190,7 +2190,7 @@ scm_i_inexact_centered_quo_and_rem (double x, double y)
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? */
+ scm_num_overflow (s_scm_centered_divide); /* or return a NaN? */
else
q = guile_NaN;
r = x - q * y;
@@ -2201,7 +2201,7 @@ scm_i_inexact_centered_quo_and_rem (double x, double 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_quo_and_rem (SCM x, SCM y)
+scm_i_bigint_centered_divide (SCM x, SCM y)
{
SCM q, r, min_r;
@@ -2254,16 +2254,16 @@ scm_i_bigint_centered_quo_and_rem (SCM x, SCM y)
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_i_slow_exact_centered_divide (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);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
+ s_scm_centered_divide);
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);
+ SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide);
else if (scm_is_true (scm_positive_p (y)))
q = scm_floor (scm_sum (scm_divide (x, y),
exactly_one_half));
@@ -2271,7 +2271,7 @@ scm_i_slow_exact_centered_quo_and_rem (SCM x, SCM y)
q = scm_ceiling (scm_difference (scm_divide (x, y),
exactly_one_half));
else
- scm_num_overflow (s_scm_centered_quo_and_rem);
+ scm_num_overflow (s_scm_centered_divide);
r = scm_difference (x, scm_product (q, y));
return scm_values (scm_list_2 (q, r));
}
diff --git a/libguile/numbers.h b/libguile/numbers.h
index 2cf3fd7..10a4f17 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -178,10 +178,10 @@ 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_divide (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_divide (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);
--
1.5.6.5
^ permalink raw reply related [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-30 22:24 ` Andy Wingo
2011-01-30 22:55 ` Ludovic Courtès
2011-01-31 17:14 ` [PATCH] Fast R6RS div/mod; improved extensibility of numerics Mark H Weaver
@ 2011-01-31 17:35 ` Mark H Weaver
2011-01-31 19:26 ` Andy Wingo
2011-01-31 20:16 ` Andy Wingo
2 siblings, 2 replies; 12+ messages in thread
From: Mark H Weaver @ 2011-01-31 17:35 UTC (permalink / raw)
To: Andy Wingo; +Cc: guile-devel
Andy Wingo <wingo@pobox.com> writes:
>> The last patch might be slightly controversial. Although it does not
>> make `integer-expt' a generic function, nonetheless it can now
>> exponentiate _any_ scheme object that can be multiplied using `*'.
>
> To me this is fine, in principle; but I wonder about causing `(expt #t
> 0)' to produce 1. It seems that in any case the range of the output
> should be continuous -- i.e. if (expt <my-matrix> 1) produces a matrix,
> and (expt <my-matrix> n) for n > 1 makes a matrix, then (expt
> <my-matrix> 0) should also produce a matrix. Or no?
When using generic arithmetic operators, there's no way to return
objects of the correct type in all corner cases like this.
Since (expt <obj> N) for exact non-negative integer N may be defined as
(* <obj> <obj> <obj> ...) with N arguments, then what you are proposing
is essentially equivalent to this:
(* <my-matrix> <my-matrix>) yields a matrix, and
(* <my-matrix>) yields a matrix, therefore
(*) should also yield a matrix.
Obviously there's no way to do this with a generic `*' operator. In
fact, you can't even do it with a special `matrix*' operator, because it
wouldn't know what size the matrix should be. You'd need a distinct `*'
operator for each size. With tensors it gets even more hairy.
Similarly for `+' and `-'.
To avoid these complications, it is conventional to allow 1 to stand for
any multiplicative identity and for 0 to stand for any additive
identity. Furthermore, it is convenient to allow any scalar to be
interpreted as the relevant multiplicative identity times that scalar.
Best,
Mark
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-31 17:35 ` Mark H Weaver
@ 2011-01-31 19:26 ` Andy Wingo
2011-01-31 20:16 ` Andy Wingo
1 sibling, 0 replies; 12+ messages in thread
From: Andy Wingo @ 2011-01-31 19:26 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
On Mon 31 Jan 2011 18:35, Mark H Weaver <mhw@netris.org> writes:
> Andy Wingo <wingo@pobox.com> writes:
>>> The last patch might be slightly controversial. Although it does not
>>> make `integer-expt' a generic function, nonetheless it can now
>>> exponentiate _any_ scheme object that can be multiplied using `*'.
>>
>> To me this is fine, in principle; but I wonder about causing `(expt #t
>> 0)' to produce 1. It seems that in any case the range of the output
>> should be continuous -- i.e. if (expt <my-matrix> 1) produces a matrix,
>> and (expt <my-matrix> n) for n > 1 makes a matrix, then (expt
>> <my-matrix> 0) should also produce a matrix. Or no?
>
> When using generic arithmetic operators, there's no way to return
> objects of the correct type in all corner cases like this.
Well... OK! I'm feeling reckless, so I'm applying and pushing your
patch. Here goes nothing!
Andy
--
http://wingolog.org/
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-31 17:35 ` Mark H Weaver
2011-01-31 19:26 ` Andy Wingo
@ 2011-01-31 20:16 ` Andy Wingo
2011-01-31 20:30 ` Mark H Weaver
2011-01-31 20:46 ` Mark H Weaver
1 sibling, 2 replies; 12+ messages in thread
From: Andy Wingo @ 2011-01-31 20:16 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
On Mon 31 Jan 2011 18:35, Mark H Weaver <mhw@netris.org> writes:
> Andy Wingo <wingo@pobox.com> writes:
>>> The last patch might be slightly controversial. Although it does not
>>> make `integer-expt' a generic function, nonetheless it can now
>>> exponentiate _any_ scheme object that can be multiplied using `*'.
>>
>> To me this is fine, in principle; but I wonder about causing `(expt #t
>> 0)' to produce 1. It seems that in any case the range of the output
>> should be continuous -- i.e. if (expt <my-matrix> 1) produces a matrix,
>> and (expt <my-matrix> n) for n > 1 makes a matrix, then (expt
>> <my-matrix> 0) should also produce a matrix. Or no?
>
> When using generic arithmetic operators, there's no way to return
> objects of the correct type in all corner cases like this.
I did apply it, but it made the test cases start to fail, because
numbers.test is loaded after ecmascript.test, and:
scheme@(guile-user)> ,L ecmascript
Happy hacking with ECMAScript! To switch back, type `,L scheme'.
ecmascript@(guile-user)> true * true;
$1 = 1
ecmascript@(guile-user)> ,L scheme
Happy hacking with Scheme! To switch back, type `,L ecmascript'.
scheme@(guile-user)> (* #t #t)
$2 = 1
I think that's the thing, anyway. It just sounds really squirrelly.
Could you not extend `expt' or `integer-expt' instead?
Andy
--
http://wingolog.org/
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-31 20:16 ` Andy Wingo
@ 2011-01-31 20:30 ` Mark H Weaver
2011-01-31 20:46 ` Andy Wingo
2011-01-31 20:46 ` Mark H Weaver
1 sibling, 1 reply; 12+ messages in thread
From: Mark H Weaver @ 2011-01-31 20:30 UTC (permalink / raw)
To: Andy Wingo; +Cc: guile-devel
Andy Wingo <wingo@pobox.com> writes:
> I did apply it, but it made the test cases start to fail, because
> numbers.test is loaded after ecmascript.test, and:
>
> scheme@(guile-user)> ,L ecmascript
> Happy hacking with ECMAScript! To switch back, type `,L scheme'.
> ecmascript@(guile-user)> true * true;
> $1 = 1
> ecmascript@(guile-user)> ,L scheme
> Happy hacking with Scheme! To switch back, type `,L ecmascript'.
> scheme@(guile-user)> (* #t #t)
> $2 = 1
How could changes in `expt' possibly affect (true * true) or (* #t #t)?
I don't understand what this has to do with my patch.
Mark
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-31 20:30 ` Mark H Weaver
@ 2011-01-31 20:46 ` Andy Wingo
0 siblings, 0 replies; 12+ messages in thread
From: Andy Wingo @ 2011-01-31 20:46 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
On Mon 31 Jan 2011 21:30, Mark H Weaver <mhw@netris.org> writes:
> Andy Wingo <wingo@pobox.com> writes:
>> I did apply it, but it made the test cases start to fail, because
>> numbers.test is loaded after ecmascript.test, and:
>>
>> scheme@(guile-user)> ,L ecmascript
>> Happy hacking with ECMAScript! To switch back, type `,L scheme'.
>> ecmascript@(guile-user)> true * true;
>> $1 = 1
>> ecmascript@(guile-user)> ,L scheme
>> Happy hacking with Scheme! To switch back, type `,L ecmascript'.
>> scheme@(guile-user)> (* #t #t)
>> $2 = 1
>
> How could changes in `expt' possibly affect (true * true) or (* #t #t)?
> I don't understand what this has to do with my patch.
Sorry, I wasn't clear.
Your patch replaced a test that expected `(expt #t 0)' to fail with one
that expected `(expt #t 2)' to fail. I assumed it was because it was
treating `(expt #t 2)' as `(* #t #t)', but I could be mistaken there.
You can see this error if you run the full ./check-guile, not just
./check-guile numbers.test.
Cheers,
Andy
--
http://wingolog.org/
^ permalink raw reply [flat|nested] 12+ messages in thread
* Re: [PATCH] Fast R6RS div/mod; improved extensibility of numerics
2011-01-31 20:16 ` Andy Wingo
2011-01-31 20:30 ` Mark H Weaver
@ 2011-01-31 20:46 ` Mark H Weaver
1 sibling, 0 replies; 12+ messages in thread
From: Mark H Weaver @ 2011-01-31 20:46 UTC (permalink / raw)
To: Andy Wingo; +Cc: guile-devel
Andy Wingo <wingo@pobox.com> writes:
> I did apply it, but it made the test cases start to fail, because
> numbers.test is loaded after ecmascript.test, and:
>
> scheme@(guile-user)> ,L ecmascript
> Happy hacking with ECMAScript! To switch back, type `,L scheme'.
> ecmascript@(guile-user)> true * true;
> $1 = 1
> ecmascript@(guile-user)> ,L scheme
> Happy hacking with Scheme! To switch back, type `,L ecmascript'.
> scheme@(guile-user)> (* #t #t)
> $2 = 1
Oh, I see. Your implementation of ecmascript extends `*' to work on
booleans, is that right? If so, then I think it naturally follows that
(expt <bool> K) should also work.
If there were a test case in numbers.test to verify that (* #t #t)
throws an exception (a test no less reasonable than (expt #t 2)), then
ecmascript would break that too, even without my recent patch.
I think it would be okay to simply remove the (expt #t 2) and
(integer-expt #t 2) tests.
> I think that's the thing, anyway. It just sounds really squirrelly.
> Could you not extend `expt' or `integer-expt' instead?
In a sense there is no need to extend any of the arithmetic operators.
We could just force people to make their own operators like `g+', `g*',
etc.
The advantage of being able to extend the normal operators is that
existing code that uses `+', `*', `expt', etc, can work with any new
data type you later add. It allows the whole system to be extensible.
I think that's a huge win.
Mark
^ permalink raw reply [flat|nested] 12+ messages in thread
end of thread, other threads:[~2011-01-31 20:46 UTC | newest]
Thread overview: 12+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2011-01-30 16:27 [PATCH] Fast R6RS div/mod; improved extensibility of numerics Mark H Weaver
2011-01-30 22:24 ` Andy Wingo
2011-01-30 22:55 ` Ludovic Courtès
2011-01-31 6:19 ` [PATCH] Rework the testing framework for number-theoretic division operators Mark H Weaver
2011-01-31 8:52 ` Andy Wingo
2011-01-31 17:14 ` [PATCH] Fast R6RS div/mod; improved extensibility of numerics Mark H Weaver
2011-01-31 17:35 ` Mark H Weaver
2011-01-31 19:26 ` Andy Wingo
2011-01-31 20:16 ` Andy Wingo
2011-01-31 20:30 ` Mark H Weaver
2011-01-31 20:46 ` Andy Wingo
2011-01-31 20:46 ` Mark H Weaver
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).