>From 847ba69f3cb9c4b6537e899159416d1c916706dc Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 13 Feb 2011 09:16:27 -0500 Subject: [PATCH 4/9] 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'. * libguile/numbers.h: Add function prototypes. * test-suite/tests/numbers.test: Add tests. * doc/ref/api-data.texi (Arithmetic): Add documentation. * NEWS: Add NEWS entry. --- NEWS | 28 + doc/ref/api-data.texi | 148 +++- libguile/numbers.c | 2187 +++++++++++++++++++++++++++++++++++++++++ libguile/numbers.h | 16 + test-suite/tests/numbers.test | 75 ++- 5 files changed, 2451 insertions(+), 3 deletions(-) diff --git a/NEWS b/NEWS index df44517..6bebbf6 100644 --- a/NEWS +++ b/NEWS @@ -23,6 +23,34 @@ instead. `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 6e37e9e..418dc85 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 address@hidden floor/ address@hidden floor-quotient address@hidden floor-remainder address@hidden ceiling/ address@hidden ceiling-quotient address@hidden ceiling-remainder address@hidden truncate/ address@hidden truncate-quotient address@hidden truncate-remainder @rnindex centered/ @rnindex centered-quotient @rnindex centered-remainder address@hidden round/ address@hidden round-quotient address@hidden round-remainder The C arithmetic functions below always takes two arguments, while the Scheme functions can take an arbitrary number. When you need to @@ -1281,6 +1293,93 @@ Note that these operators are equivalent to the R6RS operators @end lisp @end deftypefn address@hidden {Scheme Procedure} {} floor/ @var{x} @var{y} address@hidden {Scheme Procedure} {} floor-quotient @var{x} @var{y} address@hidden {Scheme Procedure} {} floor-remainder @var{x} @var{y} address@hidden {C Function} void scm_floor_divide (SCM @var{x}, SCM @var{y}, SCM address@hidden, SCM address@hidden) address@hidden {C Function} SCM scm_floor_quotient (@var{x}, @var{y}) address@hidden {C Function} SCM scm_floor_remainder (@var{x}, @var{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 address@hidden such that @address@hidden = floor(@var{x}/@var{y})} and address@hidden@var{x} = @address@hidden + @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}. + address@hidden +(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 address@hidden lisp address@hidden deftypefn + address@hidden {Scheme Procedure} {} ceiling/ @var{x} @var{y} address@hidden {Scheme Procedure} {} ceiling-quotient @var{x} @var{y} address@hidden {Scheme Procedure} {} ceiling-remainder @var{x} @var{y} address@hidden {C Function} void scm_ceiling_divide (SCM @var{x}, SCM @var{y}, SCM address@hidden, SCM address@hidden) address@hidden {C Function} SCM scm_ceiling_quotient (@var{x}, @var{y}) address@hidden {C Function} SCM scm_ceiling_remainder (@var{x}, @var{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 address@hidden such that @address@hidden = ceiling(@var{x}/@var{y})} and address@hidden@var{x} = @address@hidden + @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}. + address@hidden +(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 address@hidden lisp address@hidden deftypefn + address@hidden {Scheme Procedure} {} truncate/ @var{x} @var{y} address@hidden {Scheme Procedure} {} truncate-quotient @var{x} @var{y} address@hidden {Scheme Procedure} {} truncate-remainder @var{x} @var{y} address@hidden {C Function} void scm_truncate_divide (SCM @var{x}, SCM @var{y}, SCM address@hidden, SCM address@hidden) address@hidden {C Function} SCM scm_truncate_quotient (@var{x}, @var{y}) address@hidden {C Function} SCM scm_truncate_remainder (@var{x}, @var{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 address@hidden such that @var{q} is @address@hidden/@var{y}} rounded toward zero, +and @address@hidden = @address@hidden + @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}. + address@hidden +(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 address@hidden lisp address@hidden deftypefn + @deftypefn {Scheme Procedure} {} centered/ @var{x} @var{y} @deftypefnx {Scheme Procedure} {} centered-quotient @var{x} @var{y} @deftypefnx {Scheme Procedure} {} centered-remainder @var{x} @var{y} @@ -1313,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 deftypefn address@hidden {Scheme Procedure} {} round/ @var{x} @var{y} address@hidden {Scheme Procedure} {} round-quotient @var{x} @var{y} address@hidden {Scheme Procedure} {} round-remainder @var{x} @var{y} address@hidden {C Function} void scm_round_divide (SCM @var{x}, SCM @var{y}, SCM address@hidden, SCM address@hidden) address@hidden {C Function} SCM scm_round_quotient (@var{x}, @var{y}) address@hidden {C Function} SCM scm_round_remainder (@var{x}, @var{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 address@hidden such that @address@hidden = @address@hidden + @var{r}} and address@hidden is @address@hidden/@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 @address@hidden/@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/}, address@hidden|@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 address@hidden is even. + address@hidden +(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 address@hidden lisp address@hidden deftypefn + @node Scientific @subsubsection Scientific Functions diff --git a/libguile/numbers.c b/libguile/numbers.c index e779c42..9107c81 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -1616,6 +1616,1537 @@ scm_i_exact_rational_euclidean_divide (SCM x, SCM y, SCM *qp, SCM *rp) *rp = scm_divide (r1, scm_product (xd, yd)); } +static SCM scm_i_inexact_floor_quotient (double x, double 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), + "Return the floor of @address@hidden / @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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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)); +} + +static SCM +scm_i_exact_rational_floor_quotient (SCM x, SCM 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_exact_rational_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" + "@address@hidden = @address@hidden + @var{r}}\n" + "where @address@hidden = 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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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)); +} + +static SCM +scm_i_exact_rational_floor_remainder (SCM x, SCM 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 void scm_i_inexact_floor_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_floor_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_floor_divide, "floor/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @address@hidden = @address@hidden + @var{r}}\n" + "and @address@hidden = 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_i_floor_divide +{ + SCM q, r; + + scm_floor_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_floor_divide s_scm_i_floor_divide +#define g_scm_floor_divide g_scm_i_floor_divide + +void +scm_floor_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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; + + 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))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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); + *qp = SCM_I_MAKINUM (-1); + *rp = scm_i_normbig (r); + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + } + else if (xx <= 0) + { + *qp = SCM_INUM0; + *rp = 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); + *qp = SCM_I_MAKINUM (-1); + *rp = scm_i_normbig (r); + } + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + } + return; + } + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + 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), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG1, + s_scm_floor_divide, qp, rp); +} + +static void +scm_i_inexact_floor_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_divide); /* or return a NaN? */ + else + { + double q = floor (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_floor_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_floor_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + +static SCM scm_i_inexact_ceiling_quotient (double x, double 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), + "Return the ceiling of @address@hidden / @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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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)); +} + +static SCM +scm_i_exact_rational_ceiling_quotient (SCM x, SCM 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_exact_rational_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" + "@address@hidden = @address@hidden + @var{r}}\n" + "where @address@hidden = 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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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)); +} + +static SCM +scm_i_exact_rational_ceiling_remainder (SCM x, SCM 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 void scm_i_inexact_ceiling_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_ceiling_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_ceiling_divide, "ceiling/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @address@hidden = @address@hidden + @var{r}}\n" + "and @address@hidden = 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_i_ceiling_divide +{ + SCM q, r; + + scm_ceiling_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_ceiling_divide s_scm_i_ceiling_divide +#define g_scm_ceiling_divide g_scm_i_ceiling_divide + +void +scm_ceiling_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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; + + 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))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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)); + *qp = SCM_INUM1; + *rp = 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); + *qp = SCM_I_MAKINUM (-1); + *rp = SCM_INUM0; + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + } + else if (xx >= 0) + { + *qp = SCM_INUM0; + *rp = 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)); + *qp = SCM_INUM1; + *rp = scm_i_normbig (r); + } + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + } + return; + } + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + 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), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG1, + s_scm_ceiling_divide, qp, rp); +} + +static void +scm_i_inexact_ceiling_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_divide); /* or return a NaN? */ + else + { + double q = ceil (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_ceiling_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_ceiling_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + +static SCM scm_i_inexact_truncate_quotient (double x, double 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), + "Return @address@hidden / @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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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 (scm_c_truncate (x / y)); +} + +static SCM +scm_i_exact_rational_truncate_quotient (SCM x, SCM 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_exact_rational_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" + "@address@hidden = @address@hidden + @var{r}}\n" + "where @address@hidden = 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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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 * scm_c_truncate (x / y)); +} + +static SCM +scm_i_exact_rational_truncate_remainder (SCM x, SCM 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 void scm_i_inexact_truncate_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_truncate_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_truncate_divide, "truncate/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @address@hidden = @address@hidden + @var{r}}\n" + "and @address@hidden = 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_i_truncate_divide +{ + SCM q, r; + + scm_truncate_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_truncate_divide s_scm_i_truncate_divide +#define g_scm_truncate_divide g_scm_i_truncate_divide + +void +scm_truncate_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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); + *qp = SCM_I_MAKINUM (-1); + *rp = SCM_INUM0; + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + 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); + *qp = scm_i_normbig (q); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + 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), qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_truncate_divide, x, y, SCM_ARG1, + s_scm_truncate_divide, qp, rp); +} + +static void +scm_i_inexact_truncate_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_divide); /* or return a NaN? */ + else + { + double q, r, q1; + /* FIXME: Use trunc, after it has been imported from gnulib */ + q1 = x / y; + q = (q1 >= 0) ? floor (q1) : ceil (q1); + r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_truncate_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_truncate_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + 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_exact_rational_centered_quotient (SCM x, SCM y); @@ -2313,6 +3844,662 @@ scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp) *rp = scm_divide (r1, scm_product (xd, yd)); } +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_exact_rational_round_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return @address@hidden / @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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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); +} + +static SCM +scm_i_exact_rational_round_quotient (SCM x, SCM 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_exact_rational_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" + "@address@hidden = @address@hidden + @var{r}}, where\n" + "@var{q} is @address@hidden / @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_exact_rational_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_exact_rational_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 if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (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); + } + 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); +} + +static SCM +scm_i_exact_rational_round_remainder (SCM x, SCM 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 void scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp); +static void scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp); +static void scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_round_divide, "round/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @address@hidden = @address@hidden + @var{r}}\n" + "and @var{q} is @address@hidden / @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_i_round_divide +{ + SCM q, r; + + scm_round_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_round_divide s_scm_i_round_divide +#define g_scm_round_divide g_scm_i_round_divide + +void +scm_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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; + + 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))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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, qp, rp); + } + else if (SCM_REALP (y)) + return scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); + } + 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; + } + *qp = scm_i_normbig (q); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + else if (SCM_BIGP (y)) + return scm_i_bigint_round_divide (x, y, qp, rp); + else if (SCM_REALP (y)) + return scm_i_inexact_round_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); + } + 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), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_round_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG1, + s_scm_round_divide, qp, rp); +} + +static void +scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp) +{ + 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; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static void +scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); +} + +static void +scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_round_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1, (SCM x, SCM y, SCM rest), diff --git a/libguile/numbers.h b/libguile/numbers.h index b8529a3..bc2d4f2 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -181,9 +181,21 @@ SCM_API SCM scm_modulo (SCM x, SCM y); SCM_API void scm_euclidean_divide (SCM x, SCM y, SCM *q, SCM *r); SCM_API SCM scm_euclidean_quotient (SCM x, SCM y); SCM_API SCM scm_euclidean_remainder (SCM x, SCM y); +SCM_API void scm_floor_divide (SCM x, SCM y, SCM *q, SCM *r); +SCM_API SCM scm_floor_quotient (SCM x, SCM y); +SCM_API SCM scm_floor_remainder (SCM x, SCM y); +SCM_API void scm_ceiling_divide (SCM x, SCM y, SCM *q, SCM *r); +SCM_API SCM scm_ceiling_quotient (SCM x, SCM y); +SCM_API SCM scm_ceiling_remainder (SCM x, SCM y); +SCM_API void scm_truncate_divide (SCM x, SCM y, SCM *q, SCM *r); +SCM_API SCM scm_truncate_quotient (SCM x, SCM y); +SCM_API SCM scm_truncate_remainder (SCM x, SCM y); SCM_API void scm_centered_divide (SCM x, SCM y, SCM *q, SCM *r); SCM_API SCM scm_centered_quotient (SCM x, SCM y); SCM_API SCM scm_centered_remainder (SCM x, SCM y); +SCM_API void scm_round_divide (SCM x, SCM y, SCM *q, SCM *r); +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); @@ -200,7 +212,11 @@ SCM_API SCM scm_logcount (SCM n); SCM_API SCM scm_integer_length (SCM n); SCM_INTERNAL SCM scm_i_euclidean_divide (SCM x, SCM y); +SCM_INTERNAL SCM scm_i_floor_divide (SCM x, SCM y); +SCM_INTERNAL SCM scm_i_ceiling_divide (SCM x, SCM y); +SCM_INTERNAL SCM scm_i_truncate_divide (SCM x, SCM y); SCM_INTERNAL SCM scm_i_centered_divide (SCM x, SCM y); +SCM_INTERNAL SCM scm_i_round_divide (SCM x, SCM y); SCM_INTERNAL SCM scm_i_gcd (SCM x, SCM y, SCM rest); SCM_INTERNAL SCM scm_i_lcm (SCM x, SCM y, SCM rest); 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.7.1