summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMark H Weaver <mhw@netris.org>2011-02-13 09:16:27 -0500
committerAndy Wingo <wingo@pobox.com>2011-02-14 20:30:53 +0100
commit8f9da3406b507058475db2e8dde6fbb50d12faf7 (patch)
tree2f885066ba42c049773d1a93349126fece818990
parent03ddd15baed49dcfed0d7fc76648778b76df56d3 (diff)
downloadguile-8f9da3406b507058475db2e8dde6fbb50d12faf7.tar.gz
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.
-rw-r--r--NEWS28
-rw-r--r--doc/ref/api-data.texi148
-rw-r--r--libguile/numbers.c2187
-rw-r--r--libguile/numbers.h16
-rw-r--r--test-suite/tests/numbers.test75
5 files changed, 2451 insertions, 3 deletions
diff --git a/NEWS b/NEWS
index df44517e5..6bebbf62d 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 2faeb127a..a8cce240b 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
@@ -1281,6 +1293,93 @@ Note that these operators are equivalent to the R6RS operators
@end lisp
@end deftypefn
+@deftypefn {Scheme Procedure} {} floor/ @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} floor-quotient @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} floor-remainder @var{x} @var{y}
+@deftypefnx {C Function} void scm_floor_divide (SCM @var{x}, SCM @var{y}, SCM *@var{q}, SCM *@var{r})
+@deftypefnx {C Function} SCM scm_floor_quotient (@var{x}, @var{y})
+@deftypefnx {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
+@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 deftypefn
+
+@deftypefn {Scheme Procedure} {} ceiling/ @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} ceiling-quotient @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} ceiling-remainder @var{x} @var{y}
+@deftypefnx {C Function} void scm_ceiling_divide (SCM @var{x}, SCM @var{y}, SCM *@var{q}, SCM *@var{r})
+@deftypefnx {C Function} SCM scm_ceiling_quotient (@var{x}, @var{y})
+@deftypefnx {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
+@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 deftypefn
+
+@deftypefn {Scheme Procedure} {} truncate/ @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} truncate-quotient @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} truncate-remainder @var{x} @var{y}
+@deftypefnx {C Function} void scm_truncate_divide (SCM @var{x}, SCM @var{y}, SCM *@var{q}, SCM *@var{r})
+@deftypefnx {C Function} SCM scm_truncate_quotient (@var{x}, @var{y})
+@deftypefnx {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
+@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 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
+@deftypefn {Scheme Procedure} {} round/ @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} round-quotient @var{x} @var{y}
+@deftypefnx {Scheme Procedure} {} round-remainder @var{x} @var{y}
+@deftypefnx {C Function} void scm_round_divide (SCM @var{x}, SCM @var{y}, SCM *@var{q}, SCM *@var{r})
+@deftypefnx {C Function} SCM scm_round_quotient (@var{x}, @var{y})
+@deftypefnx {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
+@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 deftypefn
+
@node Scientific
@subsubsection Scientific Functions
diff --git a/libguile/numbers.c b/libguile/numbers.c
index e779c4220..9107c8191 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 @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_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"
+ "@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_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 @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_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 @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_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"
+ "@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_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 @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_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 @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_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"
+ "@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_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 @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_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 @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_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"
+ "@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_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 @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_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 b8529a363..bc2d4f243 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 f73818921..ef59a02ca 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?)))