From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.org!not-for-mail From: Andy Wingo Newsgroups: gmane.lisp.guile.devel Subject: Re: Simplified code for fast R6RS division operators Date: Sun, 30 Jan 2011 22:36:33 +0100 Message-ID: References: <87bp2y14zu.fsf@yeeloong.netris.org> NNTP-Posting-Host: lo.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: dough.gmane.org 1296423160 22985 80.91.229.12 (30 Jan 2011 21:32:40 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Sun, 30 Jan 2011 21:32:40 +0000 (UTC) Cc: guile-devel To: Mark H Weaver Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Sun Jan 30 22:32:34 2011 Return-path: Envelope-to: guile-devel@m.gmane.org Original-Received: from lists.gnu.org ([199.232.76.165]) by lo.gmane.org with esmtp (Exim 4.69) (envelope-from ) id 1PjesF-00063E-Mf for guile-devel@m.gmane.org; Sun, 30 Jan 2011 22:32:34 +0100 Original-Received: from localhost ([127.0.0.1]:59858 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PjesD-0008Bv-KV for guile-devel@m.gmane.org; Sun, 30 Jan 2011 16:31:53 -0500 Original-Received: from [140.186.70.92] (port=42982 helo=eggs.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1Pjes2-0008BI-Kz for guile-devel@gnu.org; Sun, 30 Jan 2011 16:31:47 -0500 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1Pjerx-0003oR-NK for guile-devel@gnu.org; Sun, 30 Jan 2011 16:31:42 -0500 Original-Received: from a-pb-sasl-sd.pobox.com ([64.74.157.62]:55094 helo=sasl.smtp.pobox.com) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1Pjerx-0003oK-G5 for guile-devel@gnu.org; Sun, 30 Jan 2011 16:31:37 -0500 Original-Received: from sasl.smtp.pobox.com (unknown [127.0.0.1]) by a-pb-sasl-sd.pobox.com (Postfix) with ESMTP id 8198643C2; Sun, 30 Jan 2011 16:32:27 -0500 (EST) DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=pobox.com; h=from:to:cc :subject:references:date:in-reply-to:message-id:mime-version :content-type; s=sasl; bh=xgNXlQnlE1RLvNtkY2gaE232Xa4=; b=q1a06C Vnd4ON5wGV0YN4AxlcY7cdcbXk3qaCHfkiUwGmlEUw1sEM3jiLjvClVfVqZByhQH mkni2OOiDsGMqDuoUXtlSA594nlQ1HKQs2vxl3424z+1eToWF2OOSgO3jhdGbMRw IZSeevgUvTg5WFifD7cMTr2wjiVQr0/IXRtxU= DomainKey-Signature: a=rsa-sha1; c=nofws; d=pobox.com; h=from:to:cc :subject:references:date:in-reply-to:message-id:mime-version :content-type; q=dns; s=sasl; b=SOelbK8O/V+GM1Q9tGEWgUFWhTtDIDVA tByLDT+PDIxaTDhTZxBYG2rQTI4r5gx0CFk7nfS9Blh+jcF5P/VuZVMg2O3MvcKw vgSh5Z2VGGJ6TnGZQ/8iOOoTFjda4hFl0mohNbEkxPx+84XLmHSgI/OqMGGQIu5E YT+OsFDNlRc= Original-Received: from a-pb-sasl-sd.pobox.com (unknown [127.0.0.1]) by a-pb-sasl-sd.pobox.com (Postfix) with ESMTP id 6B31D43C1; Sun, 30 Jan 2011 16:32:26 -0500 (EST) Original-Received: from unquote.localdomain (unknown [90.164.198.39]) (using TLSv1 with cipher DHE-RSA-AES256-SHA (256/256 bits)) (No client certificate requested) by a-pb-sasl-sd.pobox.com (Postfix) with ESMTPSA id 778F743C0; Sun, 30 Jan 2011 16:32:23 -0500 (EST) In-Reply-To: <87bp2y14zu.fsf@yeeloong.netris.org> (Mark H. Weaver's message of "Sun, 30 Jan 2011 16:00:37 -0500") User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/23.2 (gnu/linux) X-Pobox-Relay-ID: 6EB2D1F0-2CB8-11E0-B7A8-BC4EF3E828EC-02397024!a-pb-sasl-sd.pobox.com X-detected-operating-system: by eggs.gnu.org: Solaris 10 (beta) X-Received-From: 64.74.157.62 X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Developers list for Guile, the GNU extensibility library" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Original-Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: news.gmane.org gmane.lisp.guile.devel:11438 Archived-At: --=-=-= On Sun 30 Jan 2011 22:00, Mark H Weaver writes: > I'm sending this to you directly instead of the list, > because I've probably been posting too much :) The list is still probably the right thing :) It's not too much traffic for people that care about development _of_ guile. People who just care about development _with_ guile have guile-users :) > Regarding my code for the fast R6RS division operators: I realized in > the last couple hours that I could cut out some code and simplify other > bits. So here's a new slightly simpler version of that patch in case > you haven't yet started reviewed the code. Thanks for the note, I'll use the new patch (attached, for the list's benefit :). Cheers, Andy --=-=-= Content-Type: text/x-patch Content-Disposition: attachment; filename=0002-Add-two-new-sets-of-fast-quotient-and-remainder-oper.patch Content-Description: r6rs division patch >From 235e4643f63c6d0ec105d8eb606c04c5c3168f94 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 30 Jan 2011 08:48:28 -0500 Subject: [PATCH] Add two new sets of fast quotient and remainder operators * libguile/numbers.c (scm_euclidean_quo_and_rem, scm_euclidean_quotient, scm_euclidean_remainder, scm_centered_quo_and_rem, scm_centered_quotient, scm_centered_remainder): New extensible procedures `euclidean/', `euclidean-quotient', `euclidean-remainder', `centered/', `centered-quotient', `centered-remainder'. * libguile/numbers.h: Add function prototypes. * module/rnrs/base.scm: Remove incorrect stub implementations of `div', `mod', `div-and-mod', `div0', `mod0', and `div0-and-mod0'. Instead do renaming imports of `euclidean-quotient', `euclidean-remainder', `euclidean/', `centered-quotient', `centered-remainder', and `centered/', which are equivalent to the R6RS operators. * module/rnrs/arithmetic/fixnums.scm (fxdiv, fxmod, fxdiv-and-mod, fxdiv0, fxmod0, fxdiv0-and-mod0): Remove redundant checks for division by zero and unnecessary complexity. (fx+/carry): Remove unneeded calls to `inexact->exact'. * module/rnrs/arithmetic/flonums.scm (fldiv, flmod, fldiv-and-mod, fldiv0, flmod0, fldiv0-and-mod0): Remove redundant checks for division by zero and unnecessary complexity. Remove unneeded calls to `inexact->exact' and `exact->inexact' * test-suite/tests/numbers.test: (test-eqv?): New internal predicate for comparing numerical outputs with expected values. Add extensive test code for `euclidean/', `euclidean-quotient', `euclidean-remainder', `centered/', `centered-quotient', `centered-remainder'. * test-suite/tests/r6rs-arithmetic-fixnums.test: Fix some broken test cases, and remove `unresolved' test markers for `fxdiv', `fxmod', `fxdiv-and-mod', `fxdiv0', `fxmod0', and `fxdiv0-and-mod0'. * test-suite/tests/r6rs-arithmetic-flonums.test: Remove `unresolved' test markers for `fldiv', `flmod', `fldiv-and-mod', `fldiv0', `flmod0', and `fldiv0-and-mod0'. * doc/ref/api-data.texi (Arithmetic): Document `euclidean/', `euclidean-quotient', `euclidean-remainder', `centered/', `centered-quotient', and `centered-remainder'. (Operations on Integer Values): Add cross-references to `euclidean/' et al, from `quotient', `remainder', and `modulo'. * doc/ref/r6rs.texi (rnrs base): Improve documentation for `div', `mod', `div-and-mod', `div0', `mod0', and `div0-and-mod0'. Add cross-references to `euclidean/' et al. * NEWS: Add NEWS entry. --- NEWS | 29 + doc/ref/api-data.texi | 79 ++ doc/ref/r6rs.texi | 70 ++- libguile/numbers.c | 1227 ++++++++++++++++++++++++- libguile/numbers.h | 6 + module/rnrs/arithmetic/fixnums.scm | 23 +- module/rnrs/arithmetic/flonums.scm | 31 +- module/rnrs/base.scm | 23 +- test-suite/tests/numbers.test | 173 ++++- test-suite/tests/r6rs-arithmetic-fixnums.test | 23 +- test-suite/tests/r6rs-arithmetic-flonums.test | 9 +- 11 files changed, 1603 insertions(+), 90 deletions(-) diff --git a/NEWS b/NEWS index f45795e..33e49a4 100644 --- a/NEWS +++ b/NEWS @@ -12,6 +12,29 @@ Changes in 1.9.15 (since the 1.9.14 prerelease): ** Changes and bugfixes in numerics code +*** Added two new sets of fast quotient and remainder operators + +Added two new sets of fast quotient and remainder operator pairs with +different semantics than the R5RS operators. They support not only +integers, but all reals, including exact rationals and inexact +floating point numbers. + +These procedures accept two real numbers N and D, where the divisor D +must be non-zero. `euclidean-quotient' returns the integer Q and +`euclidean-remainder' returns the real R such that N = Q*D + R and +0 <= R < |D|. `euclidean/' returns both Q and R, and is more +efficient than computing each separately. Note that when D > 0, +`euclidean-quotient' returns floor(N/D), and when D < 0 it returns +ceiling(N/D). + +`centered-quotient', `centered-remainder', and `centered/' are similar +except that the range of remainders is -abs(D/2) <= R < abs(D/2), and +`centered-quotient' rounds N/D to the nearest integer. + +Note that these operators are equivalent to the R6RS integer division +operators `div', `mod', `div-and-mod', `div0', `mod0', and +`div0-and-mod0'. + *** `eqv?' and `equal?' now compare numbers equivalently scm_equal_p `equal?' now behaves equivalently to scm_eqv_p `eqv?' for @@ -64,6 +87,12 @@ NaNs are neither finite nor infinite. *** R6RS base library changes +**** `div', `mod', `div-and-mod', `div0', `mod0', `div0-and-mod0' + +Efficient versions of these R6RS division operators are now supported. +See the NEWS entry entitled `Added two new sets of fast quotient and +remainder operators' for more information. + **** `infinite?' changes `infinite?' now returns #t for non-real complex infinities, and throws diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi index 4256e18..b090782 100755 --- a/doc/ref/api-data.texi +++ b/doc/ref/api-data.texi @@ -897,6 +897,9 @@ sign as @var{n}. In all cases quotient and remainder satisfy (remainder 13 4) @result{} 1 (remainder -13 4) @result{} -1 @end lisp + +See also @code{euclidean-quotient}, @code{euclidean-remainder} and +related operations in @ref{Arithmetic}. @end deffn @c begin (texi-doc-string "guile" "modulo") @@ -911,6 +914,9 @@ sign as @var{d}. (modulo 13 -4) @result{} -3 (modulo -13 -4) @result{} -1 @end lisp + +See also @code{euclidean-quotient}, @code{euclidean-remainder} and +related operations in @ref{Arithmetic}. @end deffn @c begin (texi-doc-string "guile" "gcd") @@ -1130,6 +1136,12 @@ Returns the magnitude or angle of @var{z} as a @code{double}. @rnindex ceiling @rnindex truncate @rnindex round +@rnindex euclidean/ +@rnindex euclidean-quotient +@rnindex euclidean-remainder +@rnindex centered/ +@rnindex centered-quotient +@rnindex centered-remainder The C arithmetic functions below always takes two arguments, while the Scheme functions can take an arbitrary number. When you need to @@ -1229,6 +1241,73 @@ respectively, but these functions take and return @code{double} values. @end deftypefn +@deffn {Scheme Procedure} euclidean/ x y +@deffnx {Scheme Procedure} euclidean-quotient x y +@deffnx {Scheme Procedure} euclidean-remainder x y +@deffnx {C Function} scm_euclidean_quo_and_rem (x y) +@deffnx {C Function} scm_euclidean_quotient (x y) +@deffnx {C Function} scm_euclidean_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{euclidean-quotient} returns the +integer @var{q} and @code{euclidean-remainder} returns the real number +@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and +@math{0 <= @var{r} < abs(@var{y})}. @code{euclidean/} returns both @var{q} and +@var{r}, and is more efficient than computing each separately. Note +that when @math{@var{y} > 0}, @code{euclidean-quotient} returns +@math{floor(@var{x}/@var{y})}, otherwise it returns +@math{ceiling(@var{x}/@var{y})}. + +Note that these operators are equivalent to the R6RS operators +@code{div}, @code{mod}, and @code{div-and-mod}. + +@lisp +(euclidean-quotient 123 10) @result{} 12 +(euclidean-remainder 123 10) @result{} 3 +(euclidean/ 123 10) @result{} 12 and 3 +(euclidean/ 123 -10) @result{} -12 and 3 +(euclidean/ -123 10) @result{} -13 and 7 +(euclidean/ -123 -10) @result{} 13 and 7 +(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8 +(euclidean/ 16/3 -10/7) @result{} -3 and 22/21 +@end lisp +@end deffn + +@deffn {Scheme Procedure} centered/ x y +@deffnx {Scheme Procedure} centered-quotient x y +@deffnx {Scheme Procedure} centered-remainder x y +@deffnx {C Function} scm_centered_quo_and_rem (x y) +@deffnx {C Function} scm_centered_quotient (x y) +@deffnx {C Function} scm_centered_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{centered-quotient} returns the +integer @var{q} and @code{centered-remainder} returns the real number +@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and +@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}. @code{centered/} +returns both @var{q} and @var{r}, and is more efficient than computing +each separately. + +Note that @code{centered-quotient} returns @math{@var{x}/@var{y}} +rounded to the nearest integer. When @math{@var{x}/@var{y}} lies +exactly half-way between two integers, the tie is broken according to +the sign of @var{y}. If @math{@var{y} > 0}, ties are rounded toward +positive infinity, otherwise they are rounded toward negative infinity. +This is a consequence of the requirement that @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}. + +Note that these operators are equivalent to the R6RS operators +@code{div0}, @code{mod0}, and @code{div0-and-mod0}. + +@lisp +(centered-quotient 123 10) @result{} 12 +(centered-remainder 123 10) @result{} 3 +(centered/ 123 10) @result{} 12 and 3 +(centered/ 123 -10) @result{} -12 and 3 +(centered/ -123 10) @result{} -12 and -3 +(centered/ -123 -10) @result{} 12 and -3 +(centered/ -123.2 -63.5) @result{} 2.0 and 3.8 +(centered/ 16/3 -10/7) @result{} -4 and -8/21 +@end lisp +@end deffn + @node Scientific @subsubsection Scientific Functions diff --git a/doc/ref/r6rs.texi b/doc/ref/r6rs.texi index 5fee65f..d6d9d9f 100644 --- a/doc/ref/r6rs.texi +++ b/doc/ref/r6rs.texi @@ -1,6 +1,6 @@ @c -*-texinfo-*- @c This is part of the GNU Guile Reference Manual. -@c Copyright (C) 2010 +@c Copyright (C) 2010, 2011 @c Free Software Foundation, Inc. @c See the file guile.texi for copying conditions. @@ -464,21 +464,65 @@ grouped below by the existing manual sections to which they correspond. @xref{Arithmetic}, for documentation. @end deffn -@deffn {Scheme Procedure} div x1 x2 -@deffnx {Scheme Procedure} mod x1 x2 -@deffnx {Scheme Procedure} div-and-mod x1 x2 -These procedures implement number-theoretic division. +@rnindex div +@rnindex mod +@rnindex div-and-mod +@deffn {Scheme Procedure} div x y +@deffnx {Scheme Procedure} mod x y +@deffnx {Scheme Procedure} div-and-mod x y +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{div} returns the integer @var{q} +and @code{mod} returns the real number @var{r} such that +@math{@var{x} = @var{q}*@var{y} + @var{r}} and @math{0 <= @var{r} < abs(@var{y})}. +@code{div-and-mod} returns both @var{q} and @var{r}, and is more +efficient than computing each separately. Note that when @math{@var{y} > 0}, +@code{div} returns @math{floor(@var{x}/@var{y})}, otherwise +it returns @math{ceiling(@var{x}/@var{y})}. -@code{div-and-mod} returns two values, the respective results of -@code{(div x1 x2)} and @code{(mod x1 x2)}. +@lisp +(div 123 10) @result{} 12 +(mod 123 10) @result{} 3 +(div-and-mod 123 10) @result{} 12 and 3 +(div-and-mod 123 -10) @result{} -12 and 3 +(div-and-mod -123 10) @result{} -13 and 7 +(div-and-mod -123 -10) @result{} 13 and 7 +(div-and-mod -123.2 -63.5) @result{} 2.0 and 3.8 +(div-and-mod 16/3 -10/7) @result{} -3 and 22/21 +@end lisp @end deffn -@deffn {Scheme Procedure} div0 x1 x2 -@deffnx {Scheme Procedure} mod0 x1 x2 -@deffnx {Scheme Procedure} div0-and-mod0 x1 x2 -These procedures are similar to @code{div}, @code{mod}, and -@code{div-and-mod}, except that @code{mod0} returns values that lie -within a half-open interval centered on zero. +@rnindex div0 +@rnindex mod0 +@rnindex div0-and-mod0 +@deffn {Scheme Procedure} div0 x y +@deffnx {Scheme Procedure} mod0 x y +@deffnx {Scheme Procedure} div0-and-mod0 x y +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{div0} returns the +integer @var{q} and @code{rem0} returns the real number +@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and +@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}. @code{div0-and-mod0} +returns both @var{q} and @var{r}, and is more efficient than computing +each separately. + +Note that @code{div0} returns @math{@var{x}/@var{y}} rounded to the +nearest integer. When @math{@var{x}/@var{y}} lies exactly half-way +between two integers, the tie is broken according to the sign of +@var{y}. If @math{@var{y} > 0}, ties are rounded toward positive +infinity, otherwise they are rounded toward negative infinity. +This is a consequence of the requirement that +@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}. + +@lisp +(div0 123 10) @result{} 12 +(mod0 123 10) @result{} 3 +(div0-and-mod0 123 10) @result{} 12 and 3 +(div0-and-mod0 123 -10) @result{} -12 and 3 +(div0-and-mod0 -123 10) @result{} -12 and -3 +(div0-and-mod0 -123 -10) @result{} 12 and -3 +(div0-and-mod0 -123.2 -63.5) @result{} 2.0 and 3.8 +(div0-and-mod0 16/3 -10/7) @result{} -4 and -8/21 +@end lisp @end deffn @deffn {Scheme Procedure} exact-integer-sqrt k diff --git a/libguile/numbers.c b/libguile/numbers.c index e15a6ac..4515dc9 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -105,6 +105,7 @@ typedef scm_t_signed_bits scm_t_inum; static SCM flo0; +static SCM exactly_one_half; #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0) @@ -1054,6 +1055,1230 @@ scm_modulo (SCM x, SCM y) SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo); } +static SCM scm_i_inexact_euclidean_quotient (double x, double y); +static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{0 <= @var{r} < abs(@var{y})}.\n" + "@lisp\n" + "(euclidean-quotient 123 10) @result{} 12\n" + "(euclidean-quotient 123 -10) @result{} -12\n" + "(euclidean-quotient -123 10) @result{} -13\n" + "(euclidean-quotient -123 -10) @result{} 13\n" + "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n" + "(euclidean-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_euclidean_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_quotient); + else + { + scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum qq = xx / yy; + if (xx < qq * yy) + { + if (yy > 0) + qq--; + else + qq++; + } + return SCM_I_MAKINUM (qq); + } + } + else if (SCM_BIGP (y)) + { + if (SCM_I_INUM (x) >= 0) + return SCM_INUM0; + else + return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y))); + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quotient + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, + s_scm_euclidean_quotient); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_quotient); + else + { + SCM q = scm_i_mkbig (); + if (yy > 0) + mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy); + else + { + mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + } + scm_remember_upto_here_1 (x); + return scm_i_normbig (q); + } + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + mpz_fdiv_q (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + else + mpz_cdiv_q (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (q); + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, + s_scm_euclidean_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_euclidean_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, + s_scm_euclidean_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_euclidean_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1, + s_scm_euclidean_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_euclidean_quotient (double x, double y) +{ + if (SCM_LIKELY (y > 0)) + return scm_from_double (floor (x / y)); + else if (SCM_LIKELY (y < 0)) + return scm_from_double (ceil (x / y)); + else if (y == 0) + scm_num_overflow (s_scm_euclidean_quotient); /* or return a NaN? */ + else + return scm_nan (); +} + +/* Compute exact euclidean_quotient the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_euclidean_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1, + s_scm_euclidean_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, + s_scm_euclidean_quotient); + else if (scm_is_true (scm_positive_p (y))) + return scm_floor (scm_divide (x, y)); + else if (scm_is_true (scm_negative_p (y))) + return scm_ceiling (scm_divide (x, y)); + else + scm_num_overflow (s_scm_euclidean_quotient); +} + +static SCM scm_i_inexact_euclidean_remainder (double x, double y); +static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{0 <= @var{r} < abs(@var{y})} and\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "for some integer @var{q}.\n" + "@lisp\n" + "(euclidean-remainder 123 10) @result{} 3\n" + "(euclidean-remainder 123 -10) @result{} 3\n" + "(euclidean-remainder -123 10) @result{} 7\n" + "(euclidean-remainder -123 -10) @result{} 7\n" + "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n" + "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_euclidean_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_remainder); + else + { + scm_t_inum rr = SCM_I_INUM (x) % yy; + if (rr >= 0) + return SCM_I_MAKINUM (rr); + else if (yy > 0) + return SCM_I_MAKINUM (rr + yy); + else + return SCM_I_MAKINUM (rr - yy); + } + } + else if (SCM_BIGP (y)) + { + scm_t_inum xx = SCM_I_INUM (x); + if (xx >= 0) + return x; + else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + return scm_i_normbig (r); + } + else + { + SCM r = scm_i_mkbig (); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + return scm_i_normbig (r); + } + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_remainder + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, + s_scm_euclidean_remainder); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_remainder); + else + { + scm_t_inum rr; + if (yy < 0) + yy = -yy; + rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + SCM r = scm_i_mkbig (); + mpz_mod (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (r); + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, + s_scm_euclidean_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_euclidean_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, + s_scm_euclidean_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_euclidean_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_euclidean_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1, + s_scm_euclidean_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_euclidean_remainder (double x, double y) +{ + double q; + + /* Although it would be more efficient to use fmod here, we can't + because it would in some cases produce results inconsistent with + scm_i_inexact_euclidean_quotient, such that x != q * y + r (not + even close). In particular, when x is very close to a multiple of + y, then r might be either 0.0 or abs(y)-epsilon, but those two + cases must correspond to different choices of q. If r = 0.0 then q + must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient + chooses one and remainder chooses the other, it would be bad. This + problem was observed with x = 130.0 and y = 10/7. */ + if (SCM_LIKELY (y > 0)) + q = floor (x / y); + else if (SCM_LIKELY (y < 0)) + q = ceil (x / y); + else if (y == 0) + scm_num_overflow (s_scm_euclidean_remainder); /* or return a NaN? */ + else + return scm_nan (); + return scm_from_double (x - q * y); +} + +/* Compute exact euclidean_remainder the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_euclidean_remainder (SCM x, SCM y) +{ + SCM q; + + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1, + s_scm_euclidean_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, + s_scm_euclidean_remainder); + else if (scm_is_true (scm_positive_p (y))) + q = scm_floor (scm_divide (x, y)); + else if (scm_is_true (scm_negative_p (y))) + q = scm_ceiling (scm_divide (x, y)); + else + scm_num_overflow (s_scm_euclidean_remainder); + return scm_difference (x, scm_product (y, q)); +} + + +static SCM scm_i_inexact_euclidean_quo_and_rem (double x, double y); +static SCM scm_i_slow_exact_euclidean_quo_and_rem (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem, "euclidean/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{0 <= @var{r} < abs(@var{y})}.\n" + "@lisp\n" + "(euclidean/ 123 10) @result{} 12 and 3\n" + "(euclidean/ 123 -10) @result{} -12 and 3\n" + "(euclidean/ -123 10) @result{} -13 and 7\n" + "(euclidean/ -123 -10) @result{} 13 and 7\n" + "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_euclidean_quo_and_rem +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_quo_and_rem); + else + { + scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx - qq * yy; + if (rr < 0) + { + if (yy > 0) + { rr += yy; qq--; } + else + { rr -= yy; qq++; } + } + return scm_values (scm_list_2 (SCM_I_MAKINUM (qq), + SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + scm_t_inum xx = SCM_I_INUM (x); + if (xx >= 0) + return scm_values (scm_list_2 (SCM_INUM0, x)); + else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + return scm_values + (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r))); + } + else + { + SCM r = scm_i_mkbig (); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + return scm_values (scm_list_2 (SCM_INUM1, scm_i_normbig (r))); + } + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quo_and_rem + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_quo_and_rem (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2, + s_scm_euclidean_quo_and_rem); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_euclidean_quo_and_rem); + else + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + if (yy > 0) + mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), yy); + else + { + mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), -yy); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + } + scm_remember_upto_here_1 (x); + return scm_values (scm_list_2 (scm_i_normbig (q), + scm_i_normbig (r))); + } + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + else + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_values (scm_list_2 (scm_i_normbig (q), + scm_i_normbig (r))); + } + else if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quo_and_rem + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_euclidean_quo_and_rem (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2, + s_scm_euclidean_quo_and_rem); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_euclidean_quo_and_rem + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2, + s_scm_euclidean_quo_and_rem); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_euclidean_quo_and_rem + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_euclidean_quo_and_rem (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG1, + s_scm_euclidean_quo_and_rem); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_euclidean_quo_and_rem (double x, double y) +{ + double q, r; + + if (SCM_LIKELY (y > 0)) + q = floor (x / y); + else if (SCM_LIKELY (y < 0)) + q = ceil (x / y); + else if (y == 0) + scm_num_overflow (s_scm_euclidean_quo_and_rem); /* or return a NaN? */ + else + q = guile_NaN; + r = x - q * y; + return scm_values (scm_list_2 (scm_from_double (q), + scm_from_double (r))); +} + +/* Compute exact euclidean quotient and remainder the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_euclidean_quo_and_rem (SCM x, SCM y) +{ + SCM q, r; + + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG1, + s_scm_euclidean_quo_and_rem); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem, x, y, SCM_ARG2, + s_scm_euclidean_quo_and_rem); + else if (scm_is_true (scm_positive_p (y))) + q = scm_floor (scm_divide (x, y)); + else if (scm_is_true (scm_negative_p (y))) + q = scm_ceiling (scm_divide (x, y)); + else + scm_num_overflow (s_scm_euclidean_quo_and_rem); + r = scm_difference (x, scm_product (q, y)); + return scm_values (scm_list_2 (q, r)); +} + +static SCM scm_i_inexact_centered_quotient (double x, double y); +static SCM scm_i_bigint_centered_quotient (SCM x, SCM y); +static SCM scm_i_slow_exact_centered_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n" + "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "@lisp\n" + "(centered-quotient 123 10) @result{} 12\n" + "(centered-quotient 123 -10) @result{} -12\n" + "(centered-quotient -123 10) @result{} -12\n" + "(centered-quotient -123 -10) @result{} 12\n" + "(centered-quotient -123.2 -63.5) @result{} 2.0\n" + "(centered-quotient 16/3 -10/7) @result{} -4\n" + "@end lisp") +#define FUNC_NAME s_scm_centered_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_quotient); + else + { + scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx - qq * yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + qq++; + } + else + { + if (rr >= (1 - yy) / 2) + qq--; + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + qq--; + } + else + { + if (rr < yy / 2) + qq++; + } + } + return SCM_I_MAKINUM (qq); + } + } + else if (SCM_BIGP (y)) + { + /* Pass a denormalized bignum version of x (even though it + can fit in a fixnum) to scm_i_bigint_centered_quotient */ + return scm_i_bigint_centered_quotient + (scm_i_long2big (SCM_I_INUM (x)), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_quotient); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + else + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + if (rr < yy / 2) + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + return scm_i_normbig (q); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_quotient (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_centered_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1, + s_scm_centered_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_centered_quotient (double x, double y) +{ + if (SCM_LIKELY (y > 0)) + return scm_from_double (floor (x/y + 0.5)); + else if (SCM_LIKELY (y < 0)) + return scm_from_double (ceil (x/y - 0.5)); + else if (y == 0) + scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */ + else + return scm_nan (); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_centered_quotient (SCM x, SCM y) +{ + SCM q, r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y)/2 */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + else + { + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + scm_remember_upto_here_2 (r, min_r); + return scm_i_normbig (q); +} + +/* Compute exact centered quotient the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_centered_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1, + s_scm_centered_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + else if (scm_is_true (scm_positive_p (y))) + return scm_floor (scm_sum (scm_divide (x, y), + exactly_one_half)); + else if (scm_is_true (scm_negative_p (y))) + return scm_ceiling (scm_difference (scm_divide (x, y), + exactly_one_half)); + else + scm_num_overflow (s_scm_centered_quotient); +} + +static SCM scm_i_inexact_centered_remainder (double x, double y); +static SCM scm_i_bigint_centered_remainder (SCM x, SCM y); +static SCM scm_i_slow_exact_centered_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n" + "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "for some integer @var{q}.\n" + "@lisp\n" + "(centered-remainder 123 10) @result{} 3\n" + "(centered-remainder 123 -10) @result{} 3\n" + "(centered-remainder -123 10) @result{} -3\n" + "(centered-remainder -123 -10) @result{} -3\n" + "(centered-remainder -123.2 -63.5) @result{} 3.8\n" + "(centered-remainder 16/3 -10/7) @result{} -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_centered_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_remainder); + else + { + scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum rr = xx % yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + rr -= yy; + } + else + { + if (rr >= (1 - yy) / 2) + rr += yy; + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + rr += yy; + } + else + { + if (rr < yy / 2) + rr -= yy; + } + } + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + /* Pass a denormalized bignum version of x (even though it + can fit in a fixnum) to scm_i_bigint_centered_remainder */ + return scm_i_bigint_centered_remainder + (scm_i_long2big (SCM_I_INUM (x)), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_remainder); + else + { + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + rr += yy; + } + else + { + rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + if (rr < yy / 2) + rr -= yy; + } + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_remainder (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_centered_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1, + s_scm_centered_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_centered_remainder (double x, double y) +{ + double q; + + /* Although it would be more efficient to use fmod here, we can't + because it would in some cases produce results inconsistent with + scm_i_inexact_centered_quotient, such that x != r + q * y (not even + close). In particular, when x-y/2 is very close to a multiple of + y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those + two cases must correspond to different choices of q. If quotient + chooses one and remainder chooses the other, it would be bad. */ + if (SCM_LIKELY (y > 0)) + q = floor (x/y + 0.5); + else if (SCM_LIKELY (y < 0)) + q = ceil (x/y - 0.5); + else if (y == 0) + scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */ + else + return scm_nan (); + return scm_from_double (x - q * y); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_centered_remainder (SCM x, SCM y) +{ + SCM r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y)/2 */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_add (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + else + { + mpz_fdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_sub (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (r); +} + +/* Compute exact centered_remainder the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_centered_remainder (SCM x, SCM y) +{ + SCM q; + + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1, + s_scm_centered_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + else if (scm_is_true (scm_positive_p (y))) + q = scm_floor (scm_sum (scm_divide (x, y), exactly_one_half)); + else if (scm_is_true (scm_negative_p (y))) + q = scm_ceiling (scm_difference (scm_divide (x, y), exactly_one_half)); + else + scm_num_overflow (s_scm_centered_remainder); + return scm_difference (x, scm_product (y, q)); +} + + +static SCM scm_i_inexact_centered_quo_and_rem (double x, double y); +static SCM scm_i_bigint_centered_quo_and_rem (SCM x, SCM y); +static SCM scm_i_slow_exact_centered_quo_and_rem (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem, "centered/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "@lisp\n" + "(centered/ 123 10) @result{} 12 and 3\n" + "(centered/ 123 -10) @result{} -12 and 3\n" + "(centered/ -123 10) @result{} -12 and -3\n" + "(centered/ -123 -10) @result{} 12 and -3\n" + "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_centered_quo_and_rem +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_quo_and_rem); + else + { + scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx - qq * yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + { qq++; rr -= yy; } + } + else + { + if (rr >= (1 - yy) / 2) + { qq--; rr += yy; } + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + { qq--; rr += yy; } + } + else + { + if (rr < yy / 2) + { qq++; rr -= yy; } + } + } + return scm_values (scm_list_2 (SCM_I_MAKINUM (qq), + SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + /* Pass a denormalized bignum version of x (even though it + can fit in a fixnum) to scm_i_bigint_centered_quo_and_rem */ + return scm_i_bigint_centered_quo_and_rem + (scm_i_long2big (SCM_I_INUM (x)), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quo_and_rem + (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_quo_and_rem (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2, + s_scm_centered_quo_and_rem); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_quo_and_rem); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + { + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + rr += yy; + } + } + else + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + if (rr < yy / 2) + { + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + rr -= yy; + } + } + return scm_values (scm_list_2 (scm_i_normbig (q), + SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_quo_and_rem (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quo_and_rem + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_centered_quo_and_rem (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2, + s_scm_centered_quo_and_rem); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_quo_and_rem + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2, + s_scm_centered_quo_and_rem); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_quo_and_rem + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_centered_quo_and_rem (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG1, + s_scm_centered_quo_and_rem); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_centered_quo_and_rem (double x, double y) +{ + double q, r; + + if (SCM_LIKELY (y > 0)) + q = floor (x/y + 0.5); + else if (SCM_LIKELY (y < 0)) + q = ceil (x/y - 0.5); + else if (y == 0) + scm_num_overflow (s_scm_centered_quo_and_rem); /* or return a NaN? */ + else + q = guile_NaN; + r = x - q * y; + return scm_values (scm_list_2 (scm_from_double (q), + scm_from_double (r))); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_centered_quo_and_rem (SCM x, SCM y) +{ + SCM q, r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y/2) */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + { + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + mpz_add (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + } + else + { + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + { + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + mpz_sub (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + } + scm_remember_upto_here_2 (x, y); + return scm_values (scm_list_2 (scm_i_normbig (q), + scm_i_normbig (r))); +} + +/* Compute exact centered quotient and remainder the slow way. + We use this only if both arguments are exact, + and at least one of them is a fraction */ +static SCM +scm_i_slow_exact_centered_quo_and_rem (SCM x, SCM y) +{ + SCM q, r; + + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG1, + s_scm_centered_quo_and_rem); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem, x, y, SCM_ARG2, + s_scm_centered_quo_and_rem); + else if (scm_is_true (scm_positive_p (y))) + q = scm_floor (scm_sum (scm_divide (x, y), + exactly_one_half)); + else if (scm_is_true (scm_negative_p (y))) + q = scm_ceiling (scm_difference (scm_divide (x, y), + exactly_one_half)); + else + scm_num_overflow (s_scm_centered_quo_and_rem); + r = scm_difference (x, scm_product (q, y)); + return scm_values (scm_list_2 (q, r)); +} + + SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1, (SCM x, SCM y, SCM rest), "Return the greatest common divisor of all parameter values.\n" @@ -5356,8 +6581,6 @@ SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0, } #undef FUNC_NAME -static SCM exactly_one_half; - SCM_DEFINE (scm_round_number, "round", 1, 0, 0, (SCM x), "Round the number @var{x} towards the nearest integer. " diff --git a/libguile/numbers.h b/libguile/numbers.h index 740dc80..76d2972 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -177,6 +177,12 @@ SCM_API SCM scm_abs (SCM x); SCM_API SCM scm_quotient (SCM x, SCM y); SCM_API SCM scm_remainder (SCM x, SCM y); SCM_API SCM scm_modulo (SCM x, SCM y); +SCM_API SCM scm_euclidean_quo_and_rem (SCM x, SCM y); +SCM_API SCM scm_euclidean_quotient (SCM x, SCM y); +SCM_API SCM scm_euclidean_remainder (SCM x, SCM y); +SCM_API SCM scm_centered_quo_and_rem (SCM x, SCM y); +SCM_API SCM scm_centered_quotient (SCM x, SCM y); +SCM_API SCM scm_centered_remainder (SCM x, SCM y); SCM_API SCM scm_gcd (SCM x, SCM y); SCM_API SCM scm_lcm (SCM n1, SCM n2); SCM_API SCM scm_logand (SCM n1, SCM n2); diff --git a/module/rnrs/arithmetic/fixnums.scm b/module/rnrs/arithmetic/fixnums.scm index c1f3571..befbe9d 100644 --- a/module/rnrs/arithmetic/fixnums.scm +++ b/module/rnrs/arithmetic/fixnums.scm @@ -1,6 +1,6 @@ ;;; fixnums.scm --- The R6RS fixnums arithmetic library -;; Copyright (C) 2010 Free Software Foundation, Inc. +;; Copyright (C) 2010, 2011 Free Software Foundation, Inc. ;; ;; This library is free software; you can redistribute it and/or ;; modify it under the terms of the GNU Lesser General Public @@ -175,40 +175,33 @@ (define (fxdiv fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) - (let ((r (div fx1 fx2))) r)) + (div fx1 fx2)) (define (fxmod fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) - (let ((r (mod fx1 fx2))) r)) + (mod fx1 fx2)) (define (fxdiv-and-mod fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) (div-and-mod fx1 fx2)) (define (fxdiv0 fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) - (let ((r (div0 fx1 fx2))) r)) + (div0 fx1 fx2)) (define (fxmod0 fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) - (let ((r (mod0 fx1 fx2))) r)) + (mod0 fx1 fx2)) (define (fxdiv0-and-mod0 fx1 fx2) (assert-fixnum fx1 fx2) - (if (zero? fx2) (raise (make-assertion-violation))) - (call-with-values (lambda () (div0-and-mod0 fx1 fx2)) - (lambda (q r) (values q r)))) + (div0-and-mod0 fx1 fx2)) (define (fx+/carry fx1 fx2 fx3) (assert-fixnum fx1 fx2 fx3) (let* ((s (+ fx1 fx2 fx3)) - (s0 (mod0 s (inexact->exact (expt 2 (fixnum-width))))) - (s1 (div0 s (inexact->exact (expt 2 (fixnum-width)))))) + (s0 (mod0 s (expt 2 (fixnum-width)))) + (s1 (div0 s (expt 2 (fixnum-width))))) (values s0 s1))) (define (fx-/carry fx1 fx2 fx3) diff --git a/module/rnrs/arithmetic/flonums.scm b/module/rnrs/arithmetic/flonums.scm index 4fadbd0..b65c294 100644 --- a/module/rnrs/arithmetic/flonums.scm +++ b/module/rnrs/arithmetic/flonums.scm @@ -1,6 +1,6 @@ ;;; flonums.scm --- The R6RS flonums arithmetic library -;; Copyright (C) 2010 Free Software Foundation, Inc. +;; Copyright (C) 2010, 2011 Free Software Foundation, Inc. ;; ;; This library is free software; you can redistribute it and/or ;; modify it under the terms of the GNU Lesser General Public @@ -127,40 +127,27 @@ (define (fldiv-and-mod fl1 fl2) (assert-iflonum fl1 fl2) - (if (zero? fl2) (raise (make-assertion-violation))) - (let ((fx1 (inexact->exact fl1)) - (fx2 (inexact->exact fl2))) - (call-with-values (lambda () (div-and-mod fx1 fx2)) - (lambda (div mod) (values (exact->inexact div) - (exact->inexact mod)))))) + (div-and-mod fl1 fl2)) (define (fldiv fl1 fl2) (assert-iflonum fl1 fl2) - (if (zero? fl2) (raise (make-assertion-violation))) - (let ((fx1 (inexact->exact fl1)) - (fx2 (inexact->exact fl2))) - (exact->inexact (quotient fx1 fx2)))) + (div fl1 fl2)) (define (flmod fl1 fl2) (assert-iflonum fl1 fl2) - (if (zero? fl2) (raise (make-assertion-violation))) - (let ((fx1 (inexact->exact fl1)) - (fx2 (inexact->exact fl2))) - (exact->inexact (modulo fx1 fx2)))) + (mod fl1 fl2)) (define (fldiv0-and-mod0 fl1 fl2) (assert-iflonum fl1 fl2) - (if (zero? fl2) (raise (make-assertion-violation))) - (let* ((fx1 (inexact->exact fl1)) - (fx2 (inexact->exact fl2))) - (call-with-values (lambda () (div0-and-mod0 fx1 fx2)) - (lambda (q r) (values (real->flonum q) (real->flonum r)))))) + (div0-and-mod0 fl1 fl2)) (define (fldiv0 fl1 fl2) - (call-with-values (lambda () (fldiv0-and-mod0 fl1 fl2)) (lambda (q r) q))) + (assert-iflonum fl1 fl2) + (div0 fl1 fl2)) (define (flmod0 fl1 fl2) - (call-with-values (lambda () (fldiv0-and-mod0 fl1 fl2)) (lambda (q r) r))) + (assert-iflonum fl1 fl2) + (mod0 fl1 fl2)) (define (flnumerator fl) (assert-flonum fl) diff --git a/module/rnrs/base.scm b/module/rnrs/base.scm index 04a7e23..c81ded1 100644 --- a/module/rnrs/base.scm +++ b/module/rnrs/base.scm @@ -74,8 +74,12 @@ syntax-rules identifier-syntax) (import (rename (except (guile) error raise) - (quotient div) - (modulo mod) + (euclidean-quotient div) + (euclidean-remainder mod) + (euclidean/ div-and-mod) + (centered-quotient div0) + (centered-remainder mod0) + (centered/ div0-and-mod0) (inf? infinite?) (exact->inexact inexact) (inexact->exact exact)) @@ -119,21 +123,6 @@ (define (vector-map proc . vecs) (list->vector (apply map (cons proc (map vector->list vecs))))) - (define (div-and-mod x y) (let ((q (div x y)) (r (mod x y))) (values q r))) - - (define (div0 x y) - (call-with-values (lambda () (div0-and-mod0 x y)) (lambda (q r) q))) - - (define (mod0 x y) - (call-with-values (lambda () (div0-and-mod0 x y)) (lambda (q r) r))) - - (define (div0-and-mod0 x y) - (call-with-values (lambda () (div-and-mod x y)) - (lambda (q r) - (cond ((< r (abs (/ y 2))) (values q r)) - ((negative? y) (values (- q 1) (+ r y))) - (else (values (+ q 1) (+ r y))))))) - (define raise (@ (rnrs exceptions) raise)) (define condition diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 36e3128..9cf9202 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -17,7 +17,8 @@ (define-module (test-suite test-numbers) #:use-module (test-suite lib) - #:use-module (ice-9 documentation)) + #:use-module (ice-9 documentation) + #:use-module (srfi srfi-11)) ; let-values ;;; ;;; miscellaneous @@ -92,6 +93,35 @@ (negative? obj) (inf? obj))) +;; +;; Tolerance used by test-eqv? for inexact numbers. +;; +(define test-epsilon 1e-10) + +;; +;; Like eqv?, except that inexact finite numbers need only be within +;; test-epsilon (1e-10) to be considered equal. An exception is made +;; for zeroes, however. If X is zero, then it is tested using eqv? +;; without any allowance for imprecision. In particular, 0.0 is +;; considered distinct from -0.0. For non-real complex numbers, +;; each component is tested according to these rules. The intent +;; is that the known-correct value will be the first parameter. +;; +(define (test-eqv? x y) + (cond ((real? x) + (and (real? y) (test-real-eqv? x y))) + ((complex? x) + (and (not (real? y)) + (test-real-eqv? (real-part x) (real-part y)) + (test-real-eqv? (imag-part x) (imag-part y)))) + (else (eqv? x y)))) + +;; Auxiliary predicate used by test-eqv? +(define (test-real-eqv? x y) + (cond ((or (exact? x) (zero? x) (nan? x) (inf? x)) + (eqv? x y)) + (else (and (inexact? y) (> test-epsilon (abs (- x y))))))) + (define const-e 2.7182818284590452354) (define const-e^2 7.3890560989306502274) (define const-1/e 0.3678794411714423215) @@ -3480,3 +3510,144 @@ (pass-if "-100i swings back to 45deg down" (eqv-loosely? +7.071-7.071i (sqrt -100.0i)))) +;;; +;;; euclidean/ +;;; euclidean-quotient +;;; euclidean-remainder +;;; centered/ +;;; centered-quotient +;;; centered-remainder +;;; + +(with-test-prefix "Number-theoretic division" + + ;; Tests that (lo <= x < hi), + ;; but allowing for imprecision + ;; if x is inexact. + (define (test-within-range? lo hi x) + (if (exact? x) + (and (<= lo x) (< x hi)) + (let ((lo (- lo test-epsilon)) + (hi (+ hi test-epsilon))) + (<= lo x hi)))) + + (define (safe-euclidean-quotient x y) + (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg)) + ((zero? y) (throw 'divide-by-zero)) + ((nan? y) (nan)) + ((positive? y) (floor (/ x y))) + ((negative? y) (ceiling (/ x y))) + (else (throw 'unknown-problem)))) + + (define (safe-euclidean-remainder x y) + (- x (* y (safe-euclidean-quotient x y)))) + + (define (safe-euclidean/ x y) + (let ((q (safe-euclidean-quotient x y)) + (r (safe-euclidean-remainder x y))) + (if (not (and (eq? (exact? q) (exact? r)) + (eq? (exact? q) (and (exact? x) (exact? y))) + (test-real-eqv? r (- x (* q y))) + (or (and (integer? q) + (test-within-range? 0 (abs y) r)) + (not (finite? x)) + (not (finite? y))))) + (throw 'safe-euclidean/-is-broken (list x y q r)) + (values q r)))) + + (define (safe-centered-quotient x y) + (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg)) + ((zero? y) (throw 'divide-by-zero)) + ((nan? y) (nan)) + ((positive? y) (floor (+ 1/2 (/ x y)))) + ((negative? y) (ceiling (+ -1/2 (/ x y)))) + (else (throw 'unknown-problem)))) + + (define (safe-centered-remainder x y) + (- x (* y (safe-centered-quotient x y)))) + + (define (safe-centered/ x y) + (let ((q (safe-centered-quotient x y)) + (r (safe-centered-remainder x y))) + (if (not (and (eq? (exact? q) (exact? r)) + (eq? (exact? q) (and (exact? x) (exact? y))) + (test-real-eqv? r (- x (* q y))) + (or (and (integer? q) + (test-within-range? (* -1/2 (abs y)) + (* +1/2 (abs y)) + r)) + (not (finite? x)) + (not (finite? y))))) + (throw 'safe-centered/-is-broken (list x y q r)) + (values q r)))) + + (define test-numerators + (append + (list 123 125 127 130 3 5 10 123.2 125.0 + -123 -125 -127 -130 -3 -5 -10 -123.2 -125.0 + 127.2 130.0 123/7 125/7 127/7 130/7 + -127.2 -130.0 -123/7 -125/7 -127/7 -130/7 + 0 +0.0 -0.0 +inf.0 -inf.0 +nan.0 + most-negative-fixnum (1+ most-positive-fixnum) + (1- most-negative-fixnum)) + (apply append + (map (lambda (x) (list (* x (+ 1 most-positive-fixnum)) + (* x (+ 2 most-positive-fixnum)))) + '( 123 125 127 130 3 5 10 + -123 -125 -127 -130 -3 -5 -10))))) + + (define test-denominators + (list 10 5 10/7 127/2 10.0 63.5 + -10 -5 -10/7 -127/2 -10.0 -63.5 + +inf.0 -inf.0 +nan.0 most-negative-fixnum + (+ 1 most-positive-fixnum) (+ -1 most-negative-fixnum) + (+ 2 most-positive-fixnum) (+ -2 most-negative-fixnum))) + + (define (do-tests-1 op-name real-op safe-op) + (for-each (lambda (d) + (for-each (lambda (n) + (run-test (list op-name n d) #t + (lambda () + (test-eqv? (real-op n d) + (safe-op n d))))) + test-numerators)) + test-denominators)) + + (define (do-tests-2 op-name real-op safe-op) + (for-each (lambda (d) + (for-each (lambda (n) + (run-test (list op-name n d) #t + (lambda () + (let-values + (((q r) (safe-op n d)) + ((q1 r1) (real-op n d))) + (and (test-eqv? q q1) + (test-eqv? r r1)))))) + test-numerators)) + test-denominators)) + + (with-test-prefix "euclidean-quotient" + (do-tests-1 'euclidean-quotient + euclidean-quotient + safe-euclidean-quotient)) + (with-test-prefix "euclidean-remainder" + (do-tests-1 'euclidean-remainder + euclidean-remainder + safe-euclidean-remainder)) + (with-test-prefix "euclidean/" + (do-tests-2 'euclidean/ + euclidean/ + safe-euclidean/)) + + (with-test-prefix "centered-quotient" + (do-tests-1 'centered-quotient + centered-quotient + safe-centered-quotient)) + (with-test-prefix "centered-remainder" + (do-tests-1 'centered-remainder + centered-remainder + safe-centered-remainder)) + (with-test-prefix "centered/" + (do-tests-2 'centered/ + centered/ + safe-centered/))) diff --git a/test-suite/tests/r6rs-arithmetic-fixnums.test b/test-suite/tests/r6rs-arithmetic-fixnums.test index fed72eb..d39d544 100644 --- a/test-suite/tests/r6rs-arithmetic-fixnums.test +++ b/test-suite/tests/r6rs-arithmetic-fixnums.test @@ -1,6 +1,6 @@ ;;; arithmetic-fixnums.test --- Test suite for R6RS (rnrs arithmetic bitwise) -;; Copyright (C) 2010 Free Software Foundation, Inc. +;; Copyright (C) 2010, 2011 Free Software Foundation, Inc. ;; ;; This library is free software; you can redistribute it and/or ;; modify it under the terms of the GNU Lesser General Public @@ -121,32 +121,25 @@ (pass-if "simple" (call-with-values (lambda () (fxdiv-and-mod 123 10)) (lambda (d m) - (or (and (fx=? d 12) (fx=? m 3)) - (throw 'unresolved)))))) + (and (fx=? d 12) (fx=? m 3)))))) -(with-test-prefix "fxdiv" - (pass-if "simple" (or (fx=? (fxdiv -123 10) -13) (throw 'unresolved)))) - -(with-test-prefix "fxmod" - (pass-if "simple" (or (fx=? (fxmod -123 10) 7) (throw 'unresolved)))) +(with-test-prefix "fxdiv" (pass-if "simple" (fx=? (fxdiv -123 10) -13))) +(with-test-prefix "fxmod" (pass-if "simple" (fx=? (fxmod -123 10) 7))) (with-test-prefix "fxdiv0-and-mod0" (pass-if "simple" (call-with-values (lambda () (fxdiv0-and-mod0 -123 10)) (lambda (d m) - (or (and (fx=? d 12) (fx=? m -3)) - (throw 'unresolved)))))) - -(with-test-prefix "fxdiv0" - (pass-if "simple" (or (fx=? (fxdiv0 -123 10) 12) (throw 'unresolved)))) + (and (fx=? d -12) (fx=? m -3)))))) -(with-test-prefix "fxmod0" - (pass-if "simple" (or (fx=? (fxmod0 -123 10) -3) (throw 'unresolved)))) +(with-test-prefix "fxdiv0" (pass-if "simple" (fx=? (fxdiv0 -123 10) -12))) +(with-test-prefix "fxmod0" (pass-if "simple" (fx=? (fxmod0 -123 10) -3))) ;; Without working div and mod implementations and without any example results ;; from the spec, I have no idea what the results of these functions should ;; be. -juliang +;; UPDATE: div and mod implementations are now working properly -mhw (with-test-prefix "fx+/carry" (pass-if "simple" (throw 'unresolved))) diff --git a/test-suite/tests/r6rs-arithmetic-flonums.test b/test-suite/tests/r6rs-arithmetic-flonums.test index 873447b..af9dbbf 100644 --- a/test-suite/tests/r6rs-arithmetic-flonums.test +++ b/test-suite/tests/r6rs-arithmetic-flonums.test @@ -1,6 +1,6 @@ ;;; arithmetic-flonums.test --- Test suite for R6RS (rnrs arithmetic flonums) -;; Copyright (C) 2010 Free Software Foundation, Inc. +;; Copyright (C) 2010, 2011 Free Software Foundation, Inc. ;; ;; This library is free software; you can redistribute it and/or ;; modify it under the terms of the GNU Lesser General Public @@ -195,14 +195,13 @@ (pass-if "simple" (call-with-values (lambda () (fldiv0-and-mod0 -123.0 10.0)) (lambda (div mod) - (or (and (fl=? div -12.0) (fl=? mod -3.0)) - (throw 'unresolved)))))) + (and (fl=? div -12.0) (fl=? mod -3.0)))))) (with-test-prefix "fldiv0" - (pass-if "simple" (or (fl=? (fldiv0 -123.0 10.0) -12.0) (throw 'unresolved)))) + (pass-if "simple" (fl=? (fldiv0 -123.0 10.0) -12.0))) (with-test-prefix "flmod0" - (pass-if "simple" (or (fl=? (flmod0 -123.0 10.0) -3.0) (throw 'unresolved)))) + (pass-if "simple" (fl=? (flmod0 -123.0 10.0) -3.0))) (with-test-prefix "flnumerator" (pass-if "simple" (fl=? (flnumerator 0.5) 1.0)) -- 1.5.6.5 --=-=-= -- http://wingolog.org/ --=-=-=--