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