From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.org!not-for-mail From: Mark H Weaver Newsgroups: gmane.lisp.guile.devel Subject: [PATCH] New division operators, and optimization for fractions Date: Thu, 10 Feb 2011 18:42:18 -0500 Message-ID: <87sjvvv4md.fsf@yeeloong.netris.org> NNTP-Posting-Host: lo.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: dough.gmane.org 1297382238 7218 80.91.229.12 (10 Feb 2011 23:57:18 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Thu, 10 Feb 2011 23:57:18 +0000 (UTC) To: guile-devel@gnu.org Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Fri Feb 11 00:57:11 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 1PngM0-0007NY-76 for guile-devel@m.gmane.org; Fri, 11 Feb 2011 00:57:10 +0100 Original-Received: from localhost ([127.0.0.1]:44013 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PngLw-0004Fh-2P for guile-devel@m.gmane.org; Thu, 10 Feb 2011 18:55:12 -0500 Original-Received: from [140.186.70.92] (port=35985 helo=eggs.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PngFv-0000yd-0l for guile-devel@gnu.org; Thu, 10 Feb 2011 18:53:15 -0500 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1PngA8-0000qn-4N for guile-devel@gnu.org; Thu, 10 Feb 2011 18:43:08 -0500 Original-Received: from world.peace.net ([216.204.32.208]:57537) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1PngA6-0000qe-77 for guile-devel@gnu.org; Thu, 10 Feb 2011 18:43:00 -0500 Original-Received: from ip68-9-118-38.ri.ri.cox.net ([68.9.118.38] helo=freedomincluded) by world.peace.net with esmtpa (Exim 4.69) (envelope-from ) id 1Png9U-0004Fs-02; Thu, 10 Feb 2011 18:42:21 -0500 Original-Received: from mhw by freedomincluded with local (Exim 4.69) (envelope-from ) id 1Png9S-00067E-R5; Thu, 10 Feb 2011 18:42:19 -0500 User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/23.1 (gnu/linux) X-detected-operating-system: by eggs.gnu.org: GNU/Linux 2.6 (newer, 3) X-Received-From: 216.204.32.208 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:11550 Archived-At: --=-=-= Hello all, Here are three more patches. The first adds fast implementations of the remaining number-theoretic division operators as described in Taylor Campbell's proposal: floor/, ceiling/, truncate/, round/, as well as the single-valued variants. The third patch optimizes the case where both arguments are exact and at least one is a fraction. The second patch is needed by the third. It adds an internal C function that extracts two values from a values object. I needed it because my optimization for fractions works by calling the same division operator recursively on a subproblem involving only integers, and I need a way to extract the two values. Initially, I began to make internal versions of the divide functions that returned their two values via parameters of type (SCM *), but then I realized that in the cases where I dispatch into GOOPS, I'd still need to extract from a values object. If I've overlooked some existing mechanism for extracting multiple values from C, or if you have a better idea, please let me know and I'll rework the patch. Note that these patches should be applied after the ones in my recent post entitled "Miscellaneous fixes and improvements". Best, Mark --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0001-Add-four-new-sets-of-fast-quotient-and-remainder-ope.patch Content-Description: Add four new sets of fast quotient and remainder operators >From c5cb51d82ed23c557721a6a68a38a1deb060c0ca Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Thu, 10 Feb 2011 17:04:50 -0500 Subject: [PATCH] Add four new sets of fast quotient and remainder operators * libguile/numbers.c (scm_floor_divide, scm_floor_quotient, scm_floor_remainder, scm_ceiling_divide, scm_ceiling_quotient, scm_ceiling_remainder, scm_truncate_divide, scm_truncate_quotient, scm_truncate_remainder, scm_round_divide, scm_round_quotient, scm_round_remainder): New extensible procedures `floor/', `floor-quotient', `floor-remainder', `ceiling/', `ceiling-quotient', `ceiling-remainder', `truncate/', `truncate-quotient', `truncate-remainder', `round/', `round-quotient', and `round-remainder'. (scm_round_number, scm_floor, scm_ceiling): When applied to a fraction, use scm_round_quotient, scm_floor_quotient and scm_ceiling_quotient, respectively, for improved efficiency. * libguile/numbers.h: Add function prototypes. * test-suite/tests/numbers.test: Add tests. * doc/ref/api-data.texi (Arithmetic): Add documentation. * NEWS: Add NEWS entries for new division operators. --- NEWS | 28 + doc/ref/api-data.texi | 155 +++- libguile/numbers.c | 2222 ++++++++++++++++++++++++++++++++++++++++- libguile/numbers.h | 12 + test-suite/tests/numbers.test | 75 ++- 5 files changed, 2451 insertions(+), 41 deletions(-) diff --git a/NEWS b/NEWS index 7912259..0b47121 100644 --- a/NEWS +++ b/NEWS @@ -15,6 +15,34 @@ Changes since the 1.9.15 prerelease: `define-once' is like Lisp's `defvar': it creates a toplevel binding, but only if one does not exist already. +** Added four new sets of fast quotient and remainder operators + +Added four new sets of fast quotient and remainder operators 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. Each set of operators computes an integer quotient +Q and a real remainder R such that N = Q*D + R and |R| < |D|. They +differ only in how N/D is rounded to produce Q. + +`floor-quotient' and `floor-remainder' compute Q and R, respectively, +where Q has been rounded toward negative infinity. `floor/' returns +both Q and R, and is more efficient than computing each separately. +Note that when applied to integers, `floor-remainder' is equivalent to +the R5RS integer-only `modulo' operator. `ceiling-quotient', +`ceiling-remainder', and `ceiling/' are similar except that Q is +rounded toward positive infinity. + +For `truncate-quotient', `truncate-remainder', and `truncate/', Q is +rounded toward zero. Note that when applied to integers, +`truncate-quotient' and `truncate-remainder' are equivalent to the +R5RS integer-only operators `quotient' and `remainder'. + +For `round-quotient', `round-remainder', and `round/', Q is rounded to +the nearest integer, with ties going to the nearest even integer. + ** Improved exactness handling for complex number parsing When parsing non-real complex numbers, exactness specifiers are now diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi index fd2e7ee..7604c2c 100644 --- a/doc/ref/api-data.texi +++ b/doc/ref/api-data.texi @@ -907,7 +907,7 @@ sign as @var{n}. In all cases quotient and remainder satisfy (remainder -13 4) @result{} -1 @end lisp -See also @code{euclidean-quotient}, @code{euclidean-remainder} and +See also @code{truncate-quotient}, @code{truncate-remainder} and related operations in @ref{Arithmetic}. @end deffn @@ -924,7 +924,7 @@ sign as @var{d}. (modulo -13 -4) @result{} -1 @end lisp -See also @code{euclidean-quotient}, @code{euclidean-remainder} and +See also @code{floor-quotient}, @code{floor-remainder} and related operations in @ref{Arithmetic}. @end deffn @@ -1148,9 +1148,21 @@ Returns the magnitude or angle of @var{z} as a @code{double}. @rnindex euclidean/ @rnindex euclidean-quotient @rnindex euclidean-remainder +@rnindex floor/ +@rnindex floor-quotient +@rnindex floor-remainder +@rnindex ceiling/ +@rnindex ceiling-quotient +@rnindex ceiling-remainder +@rnindex truncate/ +@rnindex truncate-quotient +@rnindex truncate-remainder @rnindex centered/ @rnindex centered-quotient @rnindex centered-remainder +@rnindex round/ +@rnindex round-quotient +@rnindex round-remainder The C arithmetic functions below always takes two arguments, while the Scheme functions can take an arbitrary number. When you need to @@ -1260,7 +1272,7 @@ 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 +@math{0 <= @var{r} < |@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 @@ -1281,6 +1293,93 @@ Note that these operators are equivalent to the R6RS operators @end lisp @end deffn +@deffn {Scheme Procedure} floor/ x y +@deffnx {Scheme Procedure} floor-quotient x y +@deffnx {Scheme Procedure} floor-remainder x y +@deffnx {C Function} scm_floor_divide (x y) +@deffnx {C Function} scm_floor_quotient (x y) +@deffnx {C Function} scm_floor_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{floor-quotient} returns the +integer @var{q} and @code{floor-remainder} returns the real number +@var{r} such that @math{@var{q} = floor(@var{x}/@var{y})} and +@math{@var{x} = @var{q}*@var{y} + @var{r}}. @code{floor/} returns +both @var{q} and @var{r}, and is more efficient than computing each +separately. Note that @var{r}, if non-zero, will have the same sign +as @var{y}. + +When @var{x} and @var{y} are integers, @code{floor-quotient} is +equivalent to the R5RS integer-only operator @code{modulo}. + +@lisp +(floor-quotient 123 10) @result{} 12 +(floor-remainder 123 10) @result{} 3 +(floor/ 123 10) @result{} 12 and 3 +(floor/ 123 -10) @result{} -13 and -7 +(floor/ -123 10) @result{} -13 and 7 +(floor/ -123 -10) @result{} 12 and -3 +(floor/ -123.2 -63.5) @result{} 1.0 and -59.7 +(floor/ 16/3 -10/7) @result{} -4 and -8/21 +@end lisp +@end deffn + +@deffn {Scheme Procedure} ceiling/ x y +@deffnx {Scheme Procedure} ceiling-quotient x y +@deffnx {Scheme Procedure} ceiling-remainder x y +@deffnx {C Function} scm_ceiling_divide (x y) +@deffnx {C Function} scm_ceiling_quotient (x y) +@deffnx {C Function} scm_ceiling_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{ceiling-quotient} returns the +integer @var{q} and @code{ceiling-remainder} returns the real number +@var{r} such that @math{@var{q} = ceiling(@var{x}/@var{y})} and +@math{@var{x} = @var{q}*@var{y} + @var{r}}. @code{ceiling/} returns +both @var{q} and @var{r}, and is more efficient than computing each +separately. Note that @var{r}, if non-zero, will have the opposite sign +of @var{y}. + +@lisp +(ceiling-quotient 123 10) @result{} 13 +(ceiling-remainder 123 10) @result{} -7 +(ceiling/ 123 10) @result{} 13 and -7 +(ceiling/ 123 -10) @result{} -12 and 3 +(ceiling/ -123 10) @result{} -12 and -3 +(ceiling/ -123 -10) @result{} 13 and 7 +(ceiling/ -123.2 -63.5) @result{} 2.0 and 3.8 +(ceiling/ 16/3 -10/7) @result{} -3 and 22/21 +@end lisp +@end deffn + +@deffn {Scheme Procedure} truncate/ x y +@deffnx {Scheme Procedure} truncate-quotient x y +@deffnx {Scheme Procedure} truncate-remainder x y +@deffnx {C Function} scm_truncate_divide (x y) +@deffnx {C Function} scm_truncate_quotient (x y) +@deffnx {C Function} scm_truncate_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{truncate-quotient} returns the +integer @var{q} and @code{truncate-remainder} returns the real number +@var{r} such that @var{q} is @math{@var{x}/@var{y}} rounded toward zero, +and @math{@var{x} = @var{q}*@var{y} + @var{r}}. @code{truncate/} returns +both @var{q} and @var{r}, and is more efficient than computing each +separately. Note that @var{r}, if non-zero, will have the same sign +as @var{x}. + +When @var{x} and @var{y} are integers, these operators are equivalent to +the R5RS integer-only operators @code{quotient} and @code{remainder}. + +@lisp +(truncate-quotient 123 10) @result{} 12 +(truncate-remainder 123 10) @result{} 3 +(truncate/ 123 10) @result{} 12 and 3 +(truncate/ 123 -10) @result{} -12 and 3 +(truncate/ -123 10) @result{} -12 and -3 +(truncate/ -123 -10) @result{} 12 and -3 +(truncate/ -123.2 -63.5) @result{} 1.0 and -59.7 +(truncate/ 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 @@ -1291,7 +1390,7 @@ 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/} +@math{-|@var{y}/2| <= @var{r} < |@var{y}/2|}. @code{centered/} returns both @var{q} and @var{r}, and is more efficient than computing each separately. @@ -1300,7 +1399,8 @@ 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)}. +This is a consequence of the requirement that +@math{-|@var{y}/2| <= @var{r} < |@var{y}/2|}. Note that these operators are equivalent to the R6RS operators @code{div0}, @code{mod0}, and @code{div0-and-mod0}. @@ -1312,11 +1412,56 @@ Note that these operators are equivalent to the R6RS operators (centered/ 123 -10) @result{} -12 and 3 (centered/ -123 10) @result{} -12 and -3 (centered/ -123 -10) @result{} 12 and -3 +(centered/ 125 10) @result{} 13 and -5 +(centered/ 127 10) @result{} 13 and -3 +(centered/ 135 10) @result{} 14 and -5 (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 +@deffn {Scheme Procedure} round/ x y +@deffnx {Scheme Procedure} round-quotient x y +@deffnx {Scheme Procedure} round-remainder x y +@deffnx {C Function} scm_round_divide (x y) +@deffnx {C Function} scm_round_quotient (x y) +@deffnx {C Function} scm_round_remainder (x y) +These procedures accept two real numbers @var{x} and @var{y}, where the +divisor @var{y} must be non-zero. @code{round-quotient} returns the +integer @var{q} and @code{round-remainder} returns the real number +@var{r} such that @math{@var{x} = @var{q}*@var{y} + @var{r}} and +@var{q} is @math{@var{x}/@var{y}} rounded to the nearest integer, +with ties going to the nearest even integer. @code{round/} +returns both @var{q} and @var{r}, and is more efficient than computing +each separately. + +Note that @code{round/} and @code{centered/} are almost equivalent, but +their behavior differs when @math{@var{x}/@var{y}} lies exactly half-way +between two integers. In this case, @code{round/} chooses the nearest +even integer, whereas @code{centered/} chooses in such a way to satisfy +the constraint @math{-|@var{y}/2| <= @var{r} < |@var{y}/2|}, which +is stronger than the corresponding constraint for @code{round/}, +@math{-|@var{y}/2| <= @var{r} <= |@var{y}/2|}. In particular, +when @var{x} and @var{y} are integers, the number of possible remainders +returned by @code{centered/} is @math{|@var{y}|}, whereas the number of +possible remainders returned by @code{round/} is @math{|@var{y}|+1} when +@var{y} is even. + +@lisp +(round-quotient 123 10) @result{} 12 +(round-remainder 123 10) @result{} 3 +(round/ 123 10) @result{} 12 and 3 +(round/ 123 -10) @result{} -12 and 3 +(round/ -123 10) @result{} -12 and -3 +(round/ -123 -10) @result{} 12 and -3 +(round/ 125 10) @result{} 12 and 5 +(round/ 127 10) @result{} 13 and -3 +(round/ 135 10) @result{} 14 and -5 +(round/ -123.2 -63.5) @result{} 2.0 and 3.8 +(round/ 16/3 -10/7) @result{} -4 and -8/21 +@end lisp +@end deffn + @node Scientific @subsubsection Scientific Functions diff --git a/libguile/numbers.c b/libguile/numbers.c index 05840ef..90e449f 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -1596,6 +1596,1522 @@ scm_i_slow_exact_euclidean_divide (SCM x, SCM y) return scm_values (scm_list_2 (q, r)); } +static SCM scm_i_inexact_floor_quotient (double x, double y); +static SCM scm_i_slow_exact_floor_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the floor of @math{@var{x} / @var{y}}.\n" + "@lisp\n" + "(floor-quotient 123 10) @result{} 12\n" + "(floor-quotient 123 -10) @result{} -13\n" + "(floor-quotient -123 10) @result{} -13\n" + "(floor-quotient -123 -10) @result{} 12\n" + "(floor-quotient -123.2 -63.5) @result{} 1.0\n" + "(floor-quotient 16/3 -10/7) @result{} -4\n" + "@end lisp") +#define FUNC_NAME s_scm_floor_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (x); + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + scm_t_inum xx1 = xx; + scm_t_inum qq; + if (SCM_LIKELY (yy > 0)) + { + if (SCM_UNLIKELY (xx < 0)) + xx1 = xx - yy + 1; + } + else if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_floor_quotient); + else if (xx > 0) + xx1 = xx - yy - 1; + qq = xx1 / yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) + return SCM_I_MAKINUM ((xx < 0) ? -1 : 0); + else + return SCM_I_MAKINUM ((xx > 0) ? -1 : 0); + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_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_floor_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + 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_cdiv_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 (); + mpz_fdiv_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_floor_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_floor_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_floor_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_floor_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1, + s_scm_floor_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_floor_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_quotient); /* or return a NaN? */ + else + return scm_from_double (floor (x / y)); +} + +/* Compute exact floor_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_floor_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1, + s_scm_floor_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_floor_quotient); + else + return scm_floor (scm_divide (x, y)); +} + +static SCM scm_i_inexact_floor_remainder (double x, double y); +static SCM scm_i_slow_exact_floor_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{@var{q} = floor(@var{x} / @var{y})}.\n" + "@lisp\n" + "(floor-remainder 123 10) @result{} 3\n" + "(floor-remainder 123 -10) @result{} -7\n" + "(floor-remainder -123 10) @result{} 7\n" + "(floor-remainder -123 -10) @result{} -3\n" + "(floor-remainder -123.2 -63.5) @result{} -59.7\n" + "(floor-remainder 16/3 -10/7) @result{} -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_floor_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_floor_remainder); + else + { + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr < 0); + else + needs_adjustment = (rr > 0); + + if (needs_adjustment) + rr += yy; + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) + { + if (xx < 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 + return x; + } + else if (xx <= 0) + return x; + 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); + return scm_i_normbig (r); + } + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_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_floor_remainder); + else + { + scm_t_inum rr; + if (yy > 0) + rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy); + else + rr = -mpz_cdiv_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_fdiv_r (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_floor_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_floor_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_floor_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_floor_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1, + s_scm_floor_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_floor_remainder (double x, double y) +{ + /* 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_floor_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 y, but those two cases must + correspond to different choices of q. If r = 0.0 then q must be + x/y, and if r = y then q must be x/y-1. If quotient chooses one + and remainder chooses the other, it would be bad. */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_remainder); /* or return a NaN? */ + else + return scm_from_double (x - y * floor (x / y)); +} + +/* Compute exact floor_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_floor_remainder (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1, + s_scm_floor_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_floor_remainder); + else + return scm_difference + (x, scm_product (y, scm_floor (scm_divide (x, y)))); +} + + +static SCM scm_i_inexact_floor_divide (double x, double y); +static SCM scm_i_slow_exact_floor_divide (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_floor_divide, "floor/", 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{@var{q} = floor(@var{x} / @var{y})}.\n" + "@lisp\n" + "(floor/ 123 10) @result{} 12 and 3\n" + "(floor/ 123 -10) @result{} -13 and -7\n" + "(floor/ -123 10) @result{} -13 and 7\n" + "(floor/ -123 -10) @result{} 12 and -3\n" + "(floor/ -123.2 -63.5) @result{} 1.0 and -59.7\n" + "(floor/ 16/3 -10/7) @result{} -4 and -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_floor_divide +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_floor_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + int needs_adjustment; + SCM q; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr < 0); + else + needs_adjustment = (rr > 0); + + if (needs_adjustment) + { + rr += yy; + qq--; + } + + if (SCM_LIKELY (SCM_FIXABLE (qq))) + q = SCM_I_MAKINUM (qq); + else + q = scm_i_inum2big (qq); + return scm_values (scm_list_2 (q, SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) + { + if (xx < 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 + return scm_values (scm_list_2 (SCM_INUM0, x)); + } + else if (xx <= 0) + return scm_values (scm_list_2 (SCM_INUM0, x)); + 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); + return scm_values (scm_list_2 (SCM_I_MAKINUM (-1), + scm_i_normbig (r))); + } + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_divide (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide); + } + 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_floor_divide); + 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_cdiv_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 (); + 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); + return scm_values (scm_list_2 (scm_i_normbig (q), + scm_i_normbig (r))); + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_floor_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_floor_divide + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_floor_divide (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG1, + s_scm_floor_divide); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_floor_divide (double x, double y) +{ + double q, r; + + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_divide); /* or return a NaN? */ + else + q = floor (x / y); + r = x - q * y; + return scm_values (scm_list_2 (scm_from_double (q), + scm_from_double (r))); +} + +/* Compute exact floor 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_floor_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_floor_divide, x, y, SCM_ARG1, + s_scm_floor_divide); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_floor_divide); + else + q = scm_floor (scm_divide (x, y)); + r = scm_difference (x, scm_product (q, y)); + return scm_values (scm_list_2 (q, r)); +} + +static SCM scm_i_inexact_ceiling_quotient (double x, double y); +static SCM scm_i_slow_exact_ceiling_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the ceiling of @math{@var{x} / @var{y}}.\n" + "@lisp\n" + "(ceiling-quotient 123 10) @result{} 13\n" + "(ceiling-quotient 123 -10) @result{} -12\n" + "(ceiling-quotient -123 10) @result{} -12\n" + "(ceiling-quotient -123 -10) @result{} 13\n" + "(ceiling-quotient -123.2 -63.5) @result{} 2.0\n" + "(ceiling-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_ceiling_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_ceiling_quotient); + else + { + scm_t_inum xx1 = xx; + scm_t_inum qq; + if (SCM_LIKELY (yy > 0)) + { + if (SCM_LIKELY (xx >= 0)) + xx1 = xx + yy - 1; + } + else if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_quotient); + else if (xx < 0) + xx1 = xx + yy + 1; + qq = xx1 / yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 0)) + return SCM_INUM1; + else if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_I_MAKINUM (-1); + } + else + return SCM_INUM0; + } + else if (xx >= 0) + return SCM_INUM0; + else + return SCM_INUM1; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_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_ceiling_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + else + { + SCM q = scm_i_mkbig (); + if (yy > 0) + mpz_cdiv_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 (); + 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_ceiling_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_ceiling_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_ceiling_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, + s_scm_ceiling_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_ceiling_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_quotient); /* or return a NaN? */ + else + return scm_from_double (ceil (x / y)); +} + +/* Compute exact ceiling_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_ceiling_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, + s_scm_ceiling_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_ceiling_quotient); + else + return scm_ceiling (scm_divide (x, y)); +} + +static SCM scm_i_inexact_ceiling_remainder (double x, double y); +static SCM scm_i_slow_exact_ceiling_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{@var{q} = ceiling(@var{x} / @var{y})}.\n" + "@lisp\n" + "(ceiling-remainder 123 10) @result{} -7\n" + "(ceiling-remainder 123 -10) @result{} 3\n" + "(ceiling-remainder -123 10) @result{} -3\n" + "(ceiling-remainder -123 -10) @result{} 7\n" + "(ceiling-remainder -123.2 -63.5) @result{} 3.8\n" + "(ceiling-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_ceiling_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_ceiling_remainder); + else + { + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr > 0); + else + needs_adjustment = (rr < 0); + + if (needs_adjustment) + rr -= yy; + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 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); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + return scm_i_normbig (r); + } + else if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_INUM0; + } + else + return x; + } + else if (xx >= 0) + return x; + 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_ceiling_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_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_ceiling_remainder); + else + { + scm_t_inum rr; + if (yy > 0) + rr = -mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy); + else + 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_cdiv_r (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_ceiling_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_ceiling_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_ceiling_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, + s_scm_ceiling_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_ceiling_remainder (double x, double y) +{ + /* 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_ceiling_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 -y, but those two cases must + correspond to different choices of q. If r = 0.0 then q must be + x/y, and if r = -y then q must be x/y+1. If quotient chooses one + and remainder chooses the other, it would be bad. */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_remainder); /* or return a NaN? */ + else + return scm_from_double (x - y * ceil (x / y)); +} + +/* Compute exact ceiling_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_ceiling_remainder (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, + s_scm_ceiling_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_ceiling_remainder); + else + return scm_difference + (x, scm_product (y, scm_ceiling (scm_divide (x, y)))); +} + +static SCM scm_i_inexact_ceiling_divide (double x, double y); +static SCM scm_i_slow_exact_ceiling_divide (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_ceiling_divide, "ceiling/", 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{@var{q} = ceiling(@var{x} / @var{y})}.\n" + "@lisp\n" + "(ceiling/ 123 10) @result{} 13 and -7\n" + "(ceiling/ 123 -10) @result{} -12 and 3\n" + "(ceiling/ -123 10) @result{} -12 and -3\n" + "(ceiling/ -123 -10) @result{} 13 and 7\n" + "(ceiling/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(ceiling/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_ceiling_divide +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_ceiling_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + int needs_adjustment; + SCM q; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr > 0); + else + needs_adjustment = (rr < 0); + + if (needs_adjustment) + { + rr -= yy; + qq++; + } + if (SCM_LIKELY (SCM_FIXABLE (qq))) + q = SCM_I_MAKINUM (qq); + else + q = scm_i_inum2big (qq); + return scm_values (scm_list_2 (q, SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 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); + 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_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* 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 + return scm_values (scm_list_2 (SCM_INUM0, x)); + } + else if (xx >= 0) + return scm_values (scm_list_2 (SCM_INUM0, x)); + 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_ceiling_divide (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide); + } + 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_ceiling_divide); + else + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + if (yy > 0) + mpz_cdiv_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 (); + 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_ceiling_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_ceiling_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_ceiling_divide + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_ceiling_divide (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG1, + s_scm_ceiling_divide); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_ceiling_divide (double x, double y) +{ + double q, r; + + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_divide); /* or return a NaN? */ + else + q = ceil (x / y); + r = x - q * y; + return scm_values (scm_list_2 (scm_from_double (q), + scm_from_double (r))); +} + +/* Compute exact ceiling 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_ceiling_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_ceiling_divide, x, y, SCM_ARG1, + s_scm_ceiling_divide); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_ceiling_divide); + else + q = scm_ceiling (scm_divide (x, y)); + r = scm_difference (x, scm_product (q, y)); + return scm_values (scm_list_2 (q, r)); +} + +static SCM scm_i_inexact_truncate_quotient (double x, double y); +static SCM scm_i_slow_exact_truncate_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return @math{@var{x} / @var{y}} rounded toward zero.\n" + "@lisp\n" + "(truncate-quotient 123 10) @result{} 12\n" + "(truncate-quotient 123 -10) @result{} -12\n" + "(truncate-quotient -123 10) @result{} -12\n" + "(truncate-quotient -123 -10) @result{} 12\n" + "(truncate-quotient -123.2 -63.5) @result{} 1.0\n" + "(truncate-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_truncate_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_truncate_quotient); + else + { + scm_t_inum qq = xx / yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_I_MAKINUM (-1); + } + else + return SCM_INUM0; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_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_truncate_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + else + { + SCM q = scm_i_mkbig (); + if (yy > 0) + mpz_tdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy); + else + { + mpz_tdiv_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 (); + mpz_tdiv_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_truncate_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_truncate_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, + s_scm_truncate_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_truncate_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_quotient); /* or return a NaN? */ + else + return scm_from_double (trunc (x / y)); +} + +/* Compute exact truncate_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_truncate_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, + s_scm_truncate_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_truncate_quotient); + else + return scm_truncate_number (scm_divide (x, y)); +} + +static SCM scm_i_inexact_truncate_remainder (double x, double y); +static SCM scm_i_slow_exact_truncate_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{@var{q} = truncate(@var{x} / @var{y})}.\n" + "@lisp\n" + "(truncate-remainder 123 10) @result{} 3\n" + "(truncate-remainder 123 -10) @result{} 3\n" + "(truncate-remainder -123 10) @result{} -3\n" + "(truncate-remainder -123 -10) @result{} -3\n" + "(truncate-remainder -123.2 -63.5) @result{} -59.7\n" + "(truncate-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_truncate_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_truncate_remainder); + else + return SCM_I_MAKINUM (xx % yy); + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_INUM0; + } + else + return x; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_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_truncate_remainder); + else + { + scm_t_inum rr = (mpz_tdiv_ui (SCM_I_BIG_MPZ (x), + (yy > 0) ? yy : -yy) + * mpz_sgn (SCM_I_BIG_MPZ (x))); + scm_remember_upto_here_1 (x); + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + SCM r = scm_i_mkbig (); + mpz_tdiv_r (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_truncate_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_truncate_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, + s_scm_truncate_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_truncate_remainder (double x, double y) +{ + /* 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_truncate_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 sgn(x)*|y|, 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_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_remainder); /* or return a NaN? */ + else + return scm_from_double (x - y * trunc (x / y)); +} + +/* Compute exact truncate_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_truncate_remainder (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, + s_scm_truncate_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_truncate_remainder); + else + return scm_difference + (x, scm_product (y, scm_truncate_number (scm_divide (x, y)))); +} + + +static SCM scm_i_inexact_truncate_divide (double x, double y); +static SCM scm_i_slow_exact_truncate_divide (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_truncate_divide, "truncate/", 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{@var{q} = truncate(@var{x} / @var{y})}.\n" + "@lisp\n" + "(truncate/ 123 10) @result{} 12 and 3\n" + "(truncate/ 123 -10) @result{} -12 and 3\n" + "(truncate/ -123 10) @result{} -12 and -3\n" + "(truncate/ -123 -10) @result{} 12 and -3\n" + "(truncate/ -123.2 -63.5) @result{} 1.0 and -59.7\n" + "(truncate/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_truncate_divide +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_truncate_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + SCM q; + + if (SCM_LIKELY (SCM_FIXABLE (qq))) + q = SCM_I_MAKINUM (qq); + else + q = scm_i_inum2big (qq); + return scm_values (scm_list_2 (q, SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* 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 + return scm_values (scm_list_2 (SCM_INUM0, x)); + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide); + } + 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_truncate_divide); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + if (yy > 0) + rr = mpz_tdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + else + { + rr = mpz_tdiv_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)); + } + rr *= mpz_sgn (SCM_I_BIG_MPZ (x)); + scm_remember_upto_here_1 (x); + return scm_values (scm_list_2 (scm_i_normbig (q), + SCM_I_MAKINUM (rr))); + } + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + mpz_tdiv_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_truncate_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_truncate_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_divide + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_truncate_divide (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG1, + s_scm_truncate_divide); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_truncate_divide (double x, double y) +{ + double q, r; + + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_divide); /* or return a NaN? */ + else + q = trunc (x / y); + r = x - q * y; + return scm_values (scm_list_2 (scm_from_double (q), + scm_from_double (r))); +} + +/* Compute exact truncate 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_truncate_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_truncate_divide, x, y, SCM_ARG1, + s_scm_truncate_divide); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_truncate_divide); + else + q = scm_truncate_number (scm_divide (x, y)); + 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); @@ -2306,6 +3822,670 @@ scm_i_slow_exact_centered_divide (SCM x, SCM y) return scm_values (scm_list_2 (q, r)); } +static SCM scm_i_inexact_round_quotient (double x, double y); +static SCM scm_i_bigint_round_quotient (SCM x, SCM y); +static SCM scm_i_slow_exact_round_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return @math{@var{x} / @var{y}} to the nearest integer,\n" + "with ties going to the nearest even integer.\n" + "@lisp\n" + "(round-quotient 123 10) @result{} 12\n" + "(round-quotient 123 -10) @result{} -12\n" + "(round-quotient -123 10) @result{} -12\n" + "(round-quotient -123 -10) @result{} 12\n" + "(round-quotient 125 10) @result{} 12\n" + "(round-quotient 127 10) @result{} 13\n" + "(round-quotient 135 10) @result{} 14\n" + "(round-quotient -123.2 -63.5) @result{} 2.0\n" + "(round-quotient 16/3 -10/7) @result{} -4\n" + "@end lisp") +#define FUNC_NAME s_scm_round_quotient +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_round_quotient); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + + if (SCM_LIKELY (yy < 0)) + { + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + qq++; + else if (r2 <= -ay) + qq--; + } + else + { + if (r2 > ay) + qq++; + else if (r2 < -ay) + qq--; + } + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (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_round_quotient */ + return scm_i_bigint_round_quotient (scm_i_long2big (xx), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_round_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_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_round_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + int needs_adjustment; + + if (yy > 0) + { + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); + } + else + { + rr = - mpz_cdiv_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)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); + } + scm_remember_upto_here_1 (x); + if (needs_adjustment) + 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_round_quotient (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_round_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_quotient (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_round_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_round_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_round_quotient (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1, + s_scm_round_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_round_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_quotient); /* or return a NaN? */ + else + return scm_from_double (scm_c_round (x / y)); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_round_quotient (SCM x, SCM y) +{ + SCM q, r, r2; + int cmp, needs_adjustment; + + /* 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 (); + r2 = scm_i_mkbig (); + + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ + scm_remember_upto_here_2 (x, r); + + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + scm_remember_upto_here_2 (r2, y); + + if (needs_adjustment) + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + + return scm_i_normbig (q); +} + +/* Compute exact round 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_round_quotient (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1, + s_scm_round_quotient); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_round_quotient); + else + return scm_round_number (scm_divide (x, y)); +} + +static SCM scm_i_inexact_round_remainder (double x, double y); +static SCM scm_i_bigint_round_remainder (SCM x, SCM y); +static SCM scm_i_slow_exact_round_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}, where\n" + "@var{q} is @math{@var{x} / @var{y}} rounded to the\n" + "nearest integer, with ties going to the nearest\n" + "even integer.\n" + "@lisp\n" + "(round-remainder 123 10) @result{} 3\n" + "(round-remainder 123 -10) @result{} 3\n" + "(round-remainder -123 10) @result{} -3\n" + "(round-remainder -123 -10) @result{} -3\n" + "(round-remainder 125 10) @result{} 5\n" + "(round-remainder 127 10) @result{} -3\n" + "(round-remainder 135 10) @result{} -5\n" + "(round-remainder -123.2 -63.5) @result{} 3.8\n" + "(round-remainder 16/3 -10/7) @result{} -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_round_remainder +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_round_remainder); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + + if (SCM_LIKELY (yy < 0)) + { + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + rr -= yy; + else if (r2 <= -ay) + rr += yy; + } + else + { + if (r2 > ay) + rr -= yy; + else if (r2 < -ay) + 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_round_remainder */ + return scm_i_bigint_round_remainder + (scm_i_long2big (xx), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_round_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_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_round_remainder); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + int needs_adjustment; + + if (yy > 0) + { + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); + } + else + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); + } + scm_remember_upto_here_2 (x, q); + if (needs_adjustment) + rr -= yy; + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_round_remainder (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_round_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_remainder (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_round_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_round_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_round_remainder (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG1, + s_scm_round_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_round_remainder (double x, double y) +{ + /* 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_round_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), 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_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_remainder); /* or return a NaN? */ + else + { + double q = scm_c_round (x / y); + 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_round_remainder (SCM x, SCM y) +{ + SCM q, r, r2; + int cmp, needs_adjustment; + + /* 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 (); + r2 = scm_i_mkbig (); + + 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_1 (x); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ + + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + scm_remember_upto_here_2 (q, r2); + + if (needs_adjustment) + mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y)); + + scm_remember_upto_here_1 (y); + return scm_i_normbig (r); +} + +/* Compute exact round_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_round_remainder (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG1, + s_scm_round_remainder); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_round_remainder); + else + { + SCM q = scm_round_number (scm_divide (x, y)); + return scm_difference (x, scm_product (q, y)); + } +} + + +static SCM scm_i_inexact_round_divide (double x, double y); +static SCM scm_i_bigint_round_divide (SCM x, SCM y); +static SCM scm_i_slow_exact_round_divide (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 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 @var{q} is @math{@var{x} / @var{y}} rounded to the\n" + "nearest integer, with ties going to the nearest even integer.\n" + "@lisp\n" + "(round/ 123 10) @result{} 12 and 3\n" + "(round/ 123 -10) @result{} -12 and 3\n" + "(round/ -123 10) @result{} -12 and -3\n" + "(round/ -123 -10) @result{} 12 and -3\n" + "(round/ 125 10) @result{} 12 and 5\n" + "(round/ 127 10) @result{} 13 and -3\n" + "(round/ 135 10) @result{} 14 and -5\n" + "(round/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(round/ 16/3 -10/7) @result{} -4 and -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_round_divide +{ + if (SCM_LIKELY (SCM_I_INUMP (x))) + { + scm_t_inum xx = SCM_I_INUM (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_round_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + SCM q; + + if (SCM_LIKELY (yy < 0)) + { + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + { qq++; rr -= yy; } + else if (r2 <= -ay) + { qq--; rr += yy; } + } + else + { + if (r2 > ay) + { qq++; rr -= yy; } + else if (r2 < -ay) + { qq--; rr += yy; } + } + if (SCM_LIKELY (SCM_FIXABLE (qq))) + q = SCM_I_MAKINUM (qq); + else + q = scm_i_inum2big (qq); + return scm_values (scm_list_2 (q, 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_round_divide */ + return scm_i_bigint_round_divide + (scm_i_long2big (SCM_I_INUM (x)), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide); + } + 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_round_divide); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + int needs_adjustment; + + if (yy > 0) + { + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); + } + else + { + rr = - mpz_cdiv_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)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); + } + scm_remember_upto_here_1 (x); + if (needs_adjustment) + { + 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_round_divide (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_round_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_slow_exact_round_divide (x, y); + else + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_round_divide + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_round_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else + return scm_i_slow_exact_round_divide (x, y); + } + else + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG1, + s_scm_round_divide); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_round_divide (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_divide); /* or return a NaN? */ + else + { + double q = scm_c_round (x / y); + double 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_round_divide (SCM x, SCM y) +{ + SCM q, r, r2; + int cmp, needs_adjustment; + + /* 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 (); + r2 = scm_i_mkbig (); + + 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_1 (x); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ + + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + + if (needs_adjustment) + { + 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 (r2, y); + return scm_values (scm_list_2 (scm_i_normbig (q), + scm_i_normbig (r))); +} + +/* Compute exact round 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_round_divide (SCM x, SCM y) +{ + if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG1, + s_scm_round_divide); + else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide); + else if (scm_is_true (scm_zero_p (y))) + scm_num_overflow (s_scm_round_divide); + else + { + SCM q = scm_round_number (scm_divide (x, y)); + SCM 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), @@ -6727,11 +8907,11 @@ SCM_DEFINE (scm_round_number, "round", 1, 0, 0, return x; else if (SCM_REALP (x)) return scm_from_double (scm_c_round (SCM_REAL_VALUE (x))); + else if (SCM_FRACTIONP (x)) + return scm_round_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else { - /* OPTIMIZE-ME: Fraction case could be done more efficiently by a - single quotient+remainder division then examining to see which way - the rounding should go. */ SCM plus_half = scm_sum (x, exactly_one_half); SCM result = scm_floor (plus_half); /* Adjust so that the rounding is towards even. */ @@ -6754,22 +8934,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0, else if (SCM_REALP (x)) return scm_from_double (floor (SCM_REAL_VALUE (x))); else if (SCM_FRACTIONP (x)) - { - SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x), - SCM_FRACTION_DENOMINATOR (x)); - if (scm_is_false (scm_negative_p (x))) - { - /* For positive x, rounding towards zero is correct. */ - return q; - } - else - { - /* For negative x, we need to return q-1 unless x is an - integer. But fractions are never integer, per our - assumptions. */ - return scm_difference (q, SCM_INUM1); - } - } + return scm_floor_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else SCM_WTA_DISPATCH_1 (g_scm_floor, x, 1, s_scm_floor); } @@ -6785,22 +8951,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0, else if (SCM_REALP (x)) return scm_from_double (ceil (SCM_REAL_VALUE (x))); else if (SCM_FRACTIONP (x)) - { - SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x), - SCM_FRACTION_DENOMINATOR (x)); - if (scm_is_false (scm_positive_p (x))) - { - /* For negative x, rounding towards zero is correct. */ - return q; - } - else - { - /* For positive x, we need to return q+1 unless x is an - integer. But fractions are never integer, per our - assumptions. */ - return scm_sum (q, SCM_INUM1); - } - } + return scm_ceiling_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else SCM_WTA_DISPATCH_1 (g_scm_ceiling, x, 1, s_scm_ceiling); } diff --git a/libguile/numbers.h b/libguile/numbers.h index 10a4f17..66cbe63 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -181,9 +181,21 @@ SCM_API SCM scm_modulo (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_floor_divide (SCM x, SCM y); +SCM_API SCM scm_floor_quotient (SCM x, SCM y); +SCM_API SCM scm_floor_remainder (SCM x, SCM y); +SCM_API SCM scm_ceiling_divide (SCM x, SCM y); +SCM_API SCM scm_ceiling_quotient (SCM x, SCM y); +SCM_API SCM scm_ceiling_remainder (SCM x, SCM y); +SCM_API SCM scm_truncate_divide (SCM x, SCM y); +SCM_API SCM scm_truncate_quotient (SCM x, SCM y); +SCM_API SCM scm_truncate_remainder (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_round_divide (SCM x, SCM y); +SCM_API SCM scm_round_quotient (SCM x, SCM y); +SCM_API SCM scm_round_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/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index f738189..ef59a02 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4148,6 +4148,42 @@ (test-within-range? 0 <= r < (abs y))) (test-eqv? q (/ x y))))) + (define (valid-floor-answer? x y q r) + (and (eq? (exact? q) + (exact? r) + (and (exact? x) (exact? y))) + (test-eqv? r (- x (* q y))) + (if (and (finite? x) (finite? y)) + (and (integer? q) + (if (> y 0) + (test-within-range? 0 <= r < y) + (test-within-range? y < r <= 0))) + (test-eqv? q (/ x y))))) + + (define (valid-ceiling-answer? x y q r) + (and (eq? (exact? q) + (exact? r) + (and (exact? x) (exact? y))) + (test-eqv? r (- x (* q y))) + (if (and (finite? x) (finite? y)) + (and (integer? q) + (if (> y 0) + (test-within-range? (- y) < r <= 0) + (test-within-range? 0 <= r < (- y)))) + (test-eqv? q (/ x y))))) + + (define (valid-truncate-answer? x y q r) + (and (eq? (exact? q) + (exact? r) + (and (exact? x) (exact? y))) + (test-eqv? r (- x (* q y))) + (if (and (finite? x) (finite? y)) + (and (integer? q) + (if (> x 0) + (test-within-range? 0 <= r < (abs y)) + (test-within-range? (- (abs y)) < r <= 0))) + (test-eqv? q (/ x y))))) + (define (valid-centered-answer? x y q r) (and (eq? (exact? q) (exact? r) @@ -4159,6 +4195,19 @@ (* -1/2 (abs y)) <= r < (* +1/2 (abs y)))) (test-eqv? q (/ x y))))) + (define (valid-round-answer? x y q r) + (and (eq? (exact? q) + (exact? r) + (and (exact? x) (exact? y))) + (test-eqv? r (- x (* q y))) + (if (and (finite? x) (finite? y)) + (and (integer? q) + (let ((ay/2 (/ (abs y) 2))) + (if (even? q) + (test-within-range? (- ay/2) <= r <= ay/2) + (test-within-range? (- ay/2) < r < ay/2)))) + (test-eqv? q (/ x y))))) + (define (for lsts f) (apply for-each f lsts)) (define big (expt 10 (1+ (inexact->exact (ceiling (log10 fixnum-max)))))) @@ -4284,8 +4333,32 @@ euclidean-remainder valid-euclidean-answer?)) + (with-test-prefix "floor/" + (run-division-tests floor/ + floor-quotient + floor-remainder + valid-floor-answer?)) + + (with-test-prefix "ceiling/" + (run-division-tests ceiling/ + ceiling-quotient + ceiling-remainder + valid-ceiling-answer?)) + + (with-test-prefix "truncate/" + (run-division-tests truncate/ + truncate-quotient + truncate-remainder + valid-truncate-answer?)) + (with-test-prefix "centered/" (run-division-tests centered/ centered-quotient centered-remainder - valid-centered-answer?))) + valid-centered-answer?)) + + (with-test-prefix "round/" + (run-division-tests round/ + round-quotient + round-remainder + valid-round-answer?))) -- 1.5.6.5 --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0002-Added-internal-C-function-to-extract-from-values-obj.patch Content-Description: Added internal C function to extract from values object >From a729ef1b8b1905cb6c2ec73cc9b54e2e9453a2b3 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Thu, 10 Feb 2011 18:03:14 -0500 Subject: [PATCH] Added internal C function to extract from values object * libguile/values.c (scm_i_extract_values_2): New internal function that extracts two values from a values object. * libguile/values.h: Added prototype. --- libguile/values.c | 18 ++++++++++++++++++ libguile/values.h | 2 ++ 2 files changed, 20 insertions(+), 0 deletions(-) diff --git a/libguile/values.c b/libguile/values.c index 967fcd6..b9fe4cc 100644 --- a/libguile/values.c +++ b/libguile/values.c @@ -35,6 +35,24 @@ SCM scm_values_vtable; +/* OBJ must be a values object containing exactly two values. + scm_i_extract_values_2 puts those two values into *p1 and *p2. */ +void +scm_i_extract_values_2 (SCM obj, SCM *p1, SCM *p2) +{ + SCM values; + + SCM_ASSERT_TYPE (SCM_VALUESP (obj), obj, SCM_ARG1, + "scm_i_extract_values_2", "values"); + values = scm_struct_ref (obj, SCM_INUM0); + if (!scm_is_null_or_nil (SCM_CDDR (values))) + scm_wrong_type_arg_msg + ("scm_i_extract_values_2", SCM_ARG1, obj, + "a values object containing exactly two values"); + *p1 = SCM_CAR (values); + *p2 = SCM_CADR (values); +} + static SCM print_values (SCM obj, SCM pwps) { diff --git a/libguile/values.h b/libguile/values.h index 0750aec..65ad8a1 100644 --- a/libguile/values.h +++ b/libguile/values.h @@ -30,6 +30,8 @@ SCM_API SCM scm_values_vtable; #define SCM_VALUESP(x) (SCM_STRUCTP (x)\ && scm_is_eq (scm_struct_vtable (x), scm_values_vtable)) +SCM_INTERNAL void scm_i_extract_values_2 (SCM obj, SCM *p1, SCM *p2); + SCM_API SCM scm_values (SCM args); SCM_INTERNAL void scm_init_values (void); -- 1.5.6.5 --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0003-Optimize-division-operators-handling-of-fractions.patch Content-Description: Optimize division operators handling of fractions >From 5a0677e5daafa6dd10dd7612468d1ea1e38692b7 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Thu, 10 Feb 2011 18:18:34 -0500 Subject: [PATCH] Optimize division operators handling of fractions * libguile/numbers.c: (scm_euclidean_quotient, scm_euclidean_remainder, scm_euclidean_divide, scm_floor_quotient, scm_floor_remainder, scm_floor_divide, scm_ceiling_quotient, scm_ceiling_remainder, scm_ceiling_divide, scm_truncate_quotient, scm_truncate_remainder, scm_truncate_divide, scm_centered_quotient, scm_centered_remainder, scm_centered_divide, scm_round_quotient, scm_round_remainder, scm_round_divide): Optimize case where both arguments are exact and at least one is a fraction, by reducing to a subproblem involving only integers, and then adjusting the resulting remainder as needed. --- libguile/numbers.c | 622 +++++++++++++++++++++------------------------------- 1 files changed, 249 insertions(+), 373 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 90e449f..3ef8fe5 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -1070,7 +1070,7 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0, #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); +static SCM scm_i_exact_rational_euclidean_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, (SCM x, SCM y), @@ -1125,7 +1125,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_euclidean_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_quotient (x, y); + return scm_i_exact_rational_euclidean_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, s_scm_euclidean_quotient); @@ -1171,7 +1171,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, 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); + return scm_i_exact_rational_euclidean_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, s_scm_euclidean_quotient); @@ -1191,8 +1191,11 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_euclidean_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_euclidean_quotient (x, y); else - return scm_i_slow_exact_euclidean_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, + s_scm_euclidean_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1, @@ -1213,28 +1216,16 @@ scm_i_inexact_euclidean_quotient (double x, double y) 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) +scm_i_exact_rational_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); + return scm_euclidean_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } static SCM scm_i_inexact_euclidean_remainder (double x, double y); -static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_euclidean_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, (SCM x, SCM y), @@ -1294,7 +1285,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_euclidean_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_remainder (x, y); + return scm_i_exact_rational_euclidean_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, s_scm_euclidean_remainder); @@ -1329,7 +1320,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, 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); + return scm_i_exact_rational_euclidean_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, s_scm_euclidean_remainder); @@ -1349,8 +1340,11 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_euclidean_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_euclidean_remainder (x, y); else - return scm_i_slow_exact_euclidean_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, + s_scm_euclidean_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1, @@ -1383,32 +1377,19 @@ scm_i_inexact_euclidean_remainder (double x, double y) 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_i_exact_rational_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)); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_euclidean_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } static SCM scm_i_inexact_euclidean_divide (double x, double y); -static SCM scm_i_slow_exact_euclidean_divide (SCM x, SCM y); +static SCM scm_i_exact_rational_euclidean_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, (SCM x, SCM y), @@ -1477,7 +1458,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_euclidean_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_divide (x, y); + return scm_i_exact_rational_euclidean_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, s_scm_euclidean_divide); @@ -1525,7 +1506,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, 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_divide (x, y); + return scm_i_exact_rational_euclidean_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, s_scm_euclidean_divide); @@ -1545,8 +1526,11 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_euclidean_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_euclidean_divide (x, y); else - return scm_i_slow_exact_euclidean_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, + s_scm_euclidean_divide); } else SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1, @@ -1572,32 +1556,23 @@ scm_i_inexact_euclidean_divide (double x, double y) 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_divide (SCM x, SCM y) +scm_i_exact_rational_euclidean_divide (SCM x, SCM y) { - SCM q, r; + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - 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_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_divide); - r = scm_difference (x, scm_product (q, y)); + scm_i_extract_values_2 + (scm_euclidean_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); return scm_values (scm_list_2 (q, r)); } static SCM scm_i_inexact_floor_quotient (double x, double y); -static SCM scm_i_slow_exact_floor_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_floor_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, (SCM x, SCM y), @@ -1647,7 +1622,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_floor_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_quotient (x, y); + return scm_i_exact_rational_floor_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, s_scm_floor_quotient); @@ -1688,7 +1663,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, return scm_i_inexact_floor_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_quotient (x, y); + return scm_i_exact_rational_floor_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, s_scm_floor_quotient); @@ -1708,8 +1683,11 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_floor_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_quotient (x, y); else - return scm_i_slow_exact_floor_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1, @@ -1726,26 +1704,16 @@ scm_i_inexact_floor_quotient (double x, double y) return scm_from_double (floor (x / y)); } -/* Compute exact floor_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_floor_quotient (SCM x, SCM y) +scm_i_exact_rational_floor_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1, - s_scm_floor_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, - s_scm_floor_quotient); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_floor_quotient); - else - return scm_floor (scm_divide (x, y)); + return scm_floor_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } static SCM scm_i_inexact_floor_remainder (double x, double y); -static SCM scm_i_slow_exact_floor_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_floor_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, (SCM x, SCM y), @@ -1814,7 +1782,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_floor_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_remainder (x, y); + return scm_i_exact_rational_floor_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, s_scm_floor_remainder); @@ -1850,7 +1818,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, return scm_i_inexact_floor_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_remainder (x, y); + return scm_i_exact_rational_floor_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, s_scm_floor_remainder); @@ -1870,8 +1838,11 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_floor_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_remainder (x, y); else - return scm_i_slow_exact_floor_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1, @@ -1896,28 +1867,19 @@ scm_i_inexact_floor_remainder (double x, double y) return scm_from_double (x - y * floor (x / y)); } -/* Compute exact floor_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_floor_remainder (SCM x, SCM y) +scm_i_exact_rational_floor_remainder (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1, - s_scm_floor_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, - s_scm_floor_remainder); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_floor_remainder); - else - return scm_difference - (x, scm_product (y, scm_floor (scm_divide (x, y)))); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_floor_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } static SCM scm_i_inexact_floor_divide (double x, double y); -static SCM scm_i_slow_exact_floor_divide (SCM x, SCM y); +static SCM scm_i_exact_rational_floor_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_floor_divide, "floor/", 2, 0, 0, (SCM x, SCM y), @@ -1998,7 +1960,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_divide, "floor/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_floor_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_divide (x, y); + return scm_i_exact_rational_floor_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, s_scm_floor_divide); @@ -2042,7 +2004,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor_divide, "floor/", 2, 0, 0, return scm_i_inexact_floor_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_floor_divide (x, y); + return scm_i_exact_rational_floor_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, s_scm_floor_divide); @@ -2062,8 +2024,11 @@ SCM_PRIMITIVE_GENERIC (scm_floor_divide, "floor/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_floor_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y); else - return scm_i_slow_exact_floor_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide); } else SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG1, @@ -2085,30 +2050,23 @@ scm_i_inexact_floor_divide (double x, double y) scm_from_double (r))); } -/* Compute exact floor 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_floor_divide (SCM x, SCM y) +scm_i_exact_rational_floor_divide (SCM x, SCM y) { - SCM q, r; + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG1, - s_scm_floor_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_floor_divide, x, y, SCM_ARG2, - s_scm_floor_divide); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_floor_divide); - else - q = scm_floor (scm_divide (x, y)); - r = scm_difference (x, scm_product (q, y)); + scm_i_extract_values_2 + (scm_floor_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); return scm_values (scm_list_2 (q, r)); } static SCM scm_i_inexact_ceiling_quotient (double x, double y); -static SCM scm_i_slow_exact_ceiling_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_ceiling_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, (SCM x, SCM y), @@ -2178,7 +2136,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_ceiling_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_quotient (x, y); + return scm_i_exact_rational_ceiling_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, s_scm_ceiling_quotient); @@ -2219,7 +2177,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, return scm_i_inexact_ceiling_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_quotient (x, y); + return scm_i_exact_rational_ceiling_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, s_scm_ceiling_quotient); @@ -2239,8 +2197,11 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_ceiling_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_quotient (x, y); else - return scm_i_slow_exact_ceiling_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, @@ -2257,26 +2218,16 @@ scm_i_inexact_ceiling_quotient (double x, double y) return scm_from_double (ceil (x / y)); } -/* Compute exact ceiling_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_ceiling_quotient (SCM x, SCM y) +scm_i_exact_rational_ceiling_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, - s_scm_ceiling_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, - s_scm_ceiling_quotient); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_ceiling_quotient); - else - return scm_ceiling (scm_divide (x, y)); + return scm_ceiling_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } static SCM scm_i_inexact_ceiling_remainder (double x, double y); -static SCM scm_i_slow_exact_ceiling_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_ceiling_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, (SCM x, SCM y), @@ -2355,7 +2306,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_ceiling_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_remainder (x, y); + return scm_i_exact_rational_ceiling_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, s_scm_ceiling_remainder); @@ -2391,7 +2342,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, return scm_i_inexact_ceiling_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_remainder (x, y); + return scm_i_exact_rational_ceiling_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, s_scm_ceiling_remainder); @@ -2411,8 +2362,11 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_ceiling_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_remainder (x, y); else - return scm_i_slow_exact_ceiling_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, @@ -2437,27 +2391,18 @@ scm_i_inexact_ceiling_remainder (double x, double y) return scm_from_double (x - y * ceil (x / y)); } -/* Compute exact ceiling_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_ceiling_remainder (SCM x, SCM y) +scm_i_exact_rational_ceiling_remainder (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, - s_scm_ceiling_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, - s_scm_ceiling_remainder); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_ceiling_remainder); - else - return scm_difference - (x, scm_product (y, scm_ceiling (scm_divide (x, y)))); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_ceiling_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } static SCM scm_i_inexact_ceiling_divide (double x, double y); -static SCM scm_i_slow_exact_ceiling_divide (SCM x, SCM y); +static SCM scm_i_exact_rational_ceiling_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_ceiling_divide, "ceiling/", 2, 0, 0, (SCM x, SCM y), @@ -2548,7 +2493,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_divide, "ceiling/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_ceiling_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_divide (x, y); + return scm_i_exact_rational_ceiling_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, s_scm_ceiling_divide); @@ -2592,7 +2537,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_divide, "ceiling/", 2, 0, 0, return scm_i_inexact_ceiling_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_ceiling_divide (x, y); + return scm_i_exact_rational_ceiling_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, s_scm_ceiling_divide); @@ -2612,8 +2557,11 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_divide, "ceiling/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_ceiling_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y); else - return scm_i_slow_exact_ceiling_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide); } else SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG1, @@ -2635,30 +2583,23 @@ scm_i_inexact_ceiling_divide (double x, double y) scm_from_double (r))); } -/* Compute exact ceiling 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_ceiling_divide (SCM x, SCM y) +scm_i_exact_rational_ceiling_divide (SCM x, SCM y) { - SCM q, r; + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG1, - s_scm_ceiling_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, - s_scm_ceiling_divide); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_ceiling_divide); - else - q = scm_ceiling (scm_divide (x, y)); - r = scm_difference (x, scm_product (q, y)); + scm_i_extract_values_2 + (scm_ceiling_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); return scm_values (scm_list_2 (q, r)); } static SCM scm_i_inexact_truncate_quotient (double x, double y); -static SCM scm_i_slow_exact_truncate_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_truncate_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, (SCM x, SCM y), @@ -2706,7 +2647,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_truncate_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_quotient (x, y); + return scm_i_exact_rational_truncate_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, s_scm_truncate_quotient); @@ -2747,7 +2688,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, return scm_i_inexact_truncate_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_quotient (x, y); + return scm_i_exact_rational_truncate_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, s_scm_truncate_quotient); @@ -2767,8 +2708,11 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_truncate_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_quotient (x, y); else - return scm_i_slow_exact_truncate_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, @@ -2785,26 +2729,16 @@ scm_i_inexact_truncate_quotient (double x, double y) return scm_from_double (trunc (x / y)); } -/* Compute exact truncate_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_truncate_quotient (SCM x, SCM y) +scm_i_exact_rational_truncate_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, - s_scm_truncate_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, - s_scm_truncate_quotient); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_truncate_quotient); - else - return scm_truncate_number (scm_divide (x, y)); + return scm_truncate_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } static SCM scm_i_inexact_truncate_remainder (double x, double y); -static SCM scm_i_slow_exact_truncate_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_truncate_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, (SCM x, SCM y), @@ -2848,7 +2782,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_truncate_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_remainder (x, y); + return scm_i_exact_rational_truncate_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, s_scm_truncate_remainder); @@ -2882,7 +2816,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, return scm_i_inexact_truncate_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_remainder (x, y); + return scm_i_exact_rational_truncate_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, s_scm_truncate_remainder); @@ -2902,8 +2836,11 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_truncate_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_remainder (x, y); else - return scm_i_slow_exact_truncate_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, @@ -2927,28 +2864,19 @@ scm_i_inexact_truncate_remainder (double x, double y) return scm_from_double (x - y * trunc (x / y)); } -/* Compute exact truncate_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_truncate_remainder (SCM x, SCM y) +scm_i_exact_rational_truncate_remainder (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, - s_scm_truncate_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, - s_scm_truncate_remainder); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_truncate_remainder); - else - return scm_difference - (x, scm_product (y, scm_truncate_number (scm_divide (x, y)))); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_truncate_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } static SCM scm_i_inexact_truncate_divide (double x, double y); -static SCM scm_i_slow_exact_truncate_divide (SCM x, SCM y); +static SCM scm_i_exact_rational_truncate_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_truncate_divide, "truncate/", 2, 0, 0, (SCM x, SCM y), @@ -3002,7 +2930,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_divide, "truncate/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_truncate_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_divide (x, y); + return scm_i_exact_rational_truncate_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, s_scm_truncate_divide); @@ -3047,7 +2975,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_divide, "truncate/", 2, 0, 0, return scm_i_inexact_truncate_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_truncate_divide (x, y); + return scm_i_exact_rational_truncate_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, s_scm_truncate_divide); @@ -3067,8 +2995,11 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_divide, "truncate/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_truncate_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y); else - return scm_i_slow_exact_truncate_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide); } else SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG1, @@ -3090,31 +3021,24 @@ scm_i_inexact_truncate_divide (double x, double y) scm_from_double (r))); } -/* Compute exact truncate 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_truncate_divide (SCM x, SCM y) +scm_i_exact_rational_truncate_divide (SCM x, SCM y) { - SCM q, r; + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG1, - s_scm_truncate_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_truncate_divide, x, y, SCM_ARG2, - s_scm_truncate_divide); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_truncate_divide); - else - q = scm_truncate_number (scm_divide (x, y)); - r = scm_difference (x, scm_product (q, y)); + scm_i_extract_values_2 + (scm_truncate_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); 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); +static SCM scm_i_exact_rational_centered_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, (SCM x, SCM y), @@ -3184,7 +3108,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_centered_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_quotient (x, y); + return scm_i_exact_rational_centered_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, s_scm_centered_quotient); @@ -3233,7 +3157,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, 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); + return scm_i_exact_rational_centered_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, s_scm_centered_quotient); @@ -3253,8 +3177,11 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_centered_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_quotient (x, y); else - return scm_i_slow_exact_centered_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1, @@ -3318,31 +3245,17 @@ scm_i_bigint_centered_quotient (SCM x, SCM y) 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) +scm_i_exact_rational_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); + return scm_centered_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } 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); +static SCM scm_i_exact_rational_centered_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, (SCM x, SCM y), @@ -3409,7 +3322,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_centered_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_remainder (x, y); + return scm_i_exact_rational_centered_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, s_scm_centered_remainder); @@ -3450,7 +3363,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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); + return scm_i_exact_rational_centered_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, s_scm_centered_remainder); @@ -3470,8 +3383,11 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_centered_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_remainder (x, y); else - return scm_i_slow_exact_centered_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1, @@ -3544,33 +3460,20 @@ scm_i_bigint_centered_remainder (SCM x, SCM 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_i_exact_rational_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)); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_centered_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } 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); +static SCM scm_i_exact_rational_centered_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, (SCM x, SCM y), @@ -3643,7 +3546,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_divide (x, y); + return scm_i_exact_rational_centered_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, s_scm_centered_divide); @@ -3697,7 +3600,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, 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_divide (x, y); + return scm_i_exact_rational_centered_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, s_scm_centered_divide); @@ -3708,7 +3611,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_inexact_centered_divide (SCM_REAL_VALUE (x), scm_to_double (y)); - else + else SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, s_scm_centered_divide); } @@ -3717,8 +3620,11 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_centered_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_divide (x, y); else - return scm_i_slow_exact_centered_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, + s_scm_centered_divide); } else SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1, @@ -3796,35 +3702,24 @@ scm_i_bigint_centered_divide (SCM x, SCM y) 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_divide (SCM x, SCM y) +scm_i_exact_rational_centered_divide (SCM x, SCM y) { - SCM q, r; + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - 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_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)); - 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_divide); - r = scm_difference (x, scm_product (q, y)); + scm_i_extract_values_2 + (scm_centered_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); return scm_values (scm_list_2 (q, r)); } static SCM scm_i_inexact_round_quotient (double x, double y); static SCM scm_i_bigint_round_quotient (SCM x, SCM y); -static SCM scm_i_slow_exact_round_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_round_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, (SCM x, SCM y), @@ -3893,7 +3788,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_round_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_quotient (x, y); + return scm_i_exact_rational_round_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, s_scm_round_quotient); @@ -3944,7 +3839,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, return scm_i_inexact_round_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_quotient (x, y); + return scm_i_exact_rational_round_quotient (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, s_scm_round_quotient); @@ -3964,8 +3859,11 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_round_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_quotient (x, y); else - return scm_i_slow_exact_round_quotient (x, y); + SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1, @@ -4014,27 +3912,17 @@ scm_i_bigint_round_quotient (SCM x, SCM y) return scm_i_normbig (q); } -/* Compute exact round 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_round_quotient (SCM x, SCM y) +scm_i_exact_rational_round_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1, - s_scm_round_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, - s_scm_round_quotient); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_round_quotient); - else - return scm_round_number (scm_divide (x, y)); + return scm_round_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } static SCM scm_i_inexact_round_remainder (double x, double y); static SCM scm_i_bigint_round_remainder (SCM x, SCM y); -static SCM scm_i_slow_exact_round_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_round_remainder (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, (SCM x, SCM y), @@ -4104,7 +3992,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_round_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_remainder (x, y); + return scm_i_exact_rational_round_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, s_scm_round_remainder); @@ -4152,7 +4040,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, return scm_i_inexact_round_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_remainder (x, y); + return scm_i_exact_rational_round_remainder (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, s_scm_round_remainder); @@ -4172,8 +4060,11 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_round_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_remainder (x, y); else - return scm_i_slow_exact_round_remainder (x, y); + SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG1, @@ -4234,31 +4125,20 @@ scm_i_bigint_round_remainder (SCM x, SCM y) return scm_i_normbig (r); } -/* Compute exact round_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_round_remainder (SCM x, SCM y) +scm_i_exact_rational_round_remainder (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG1, - s_scm_round_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_round_remainder, x, y, SCM_ARG2, - s_scm_round_remainder); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_round_remainder); - else - { - SCM q = scm_round_number (scm_divide (x, y)); - return scm_difference (x, scm_product (q, y)); - } + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_round_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } static SCM scm_i_inexact_round_divide (double x, double y); static SCM scm_i_bigint_round_divide (SCM x, SCM y); -static SCM scm_i_slow_exact_round_divide (SCM x, SCM y); +static SCM scm_i_exact_rational_round_divide (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 2, 0, 0, (SCM x, SCM y), @@ -4332,7 +4212,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 2, 0, 0, else if (SCM_REALP (y)) return scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_divide (x, y); + return scm_i_exact_rational_round_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, s_scm_round_divide); @@ -4385,7 +4265,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 2, 0, 0, return scm_i_inexact_round_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_round_divide (x, y); + return scm_i_exact_rational_round_divide (x, y); else SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, s_scm_round_divide); @@ -4396,7 +4276,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 2, 0, 0, SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_inexact_round_divide (SCM_REAL_VALUE (x), scm_to_double (y)); - else + else SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, s_scm_round_divide); } @@ -4405,8 +4285,11 @@ SCM_PRIMITIVE_GENERIC (scm_round_divide, "round/", 2, 0, 0, if (SCM_REALP (y)) return scm_i_inexact_round_divide (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_divide (x, y); else - return scm_i_slow_exact_round_divide (x, y); + SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide); } else SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG1, @@ -4464,26 +4347,19 @@ scm_i_bigint_round_divide (SCM x, SCM y) scm_i_normbig (r))); } -/* Compute exact round 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_round_divide (SCM x, SCM y) +scm_i_exact_rational_round_divide (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG1, - s_scm_round_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_round_divide, x, y, SCM_ARG2, - s_scm_round_divide); - else if (scm_is_true (scm_zero_p (y))) - scm_num_overflow (s_scm_round_divide); - else - { - SCM q = scm_round_number (scm_divide (x, y)); - SCM r = scm_difference (x, scm_product (q, y)); - return scm_values (scm_list_2 (q, r)); - } + SCM q, r, r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_i_extract_values_2 + (scm_round_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)), + &q, &r1); + r = scm_divide (r1, scm_product (xd, yd)); + return scm_values (scm_list_2 (q, r)); } -- 1.5.6.5 --=-=-=--