# include <config.h>
#endif
+#include <verify.h>
+
#include <math.h>
#include <string.h>
#include <unicase.h>
#ifndef M_LOG10E
#define M_LOG10E 0.43429448190325182765
#endif
+#ifndef M_LN2
+#define M_LN2 0.69314718055994530942
+#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/* the macro above will not work as is with fractions */
+/* Default to 1, because as we used to hard-code `free' as the
+ deallocator, we know that overriding these functions with
+ instrumented `malloc' / `free' is OK. */
+int scm_install_gmp_memory_functions = 1;
static SCM flo0;
static SCM exactly_one_half;
+static SCM flo_log10e;
#define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
static double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); }
#endif
-/* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
- an explicit check. In some future gmp (don't know what version number),
- mpz_cmp_d is supposed to do this itself. */
+/* mpz_cmp_d in GMP before 4.2 didn't recognise infinities, so
+ xmpz_cmp_d uses an explicit check. Starting with GMP 4.2 (released
+ in March 2006), mpz_cmp_d now handles infinities properly. */
#if 1
#define xmpz_cmp_d(z, d) \
(isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
#if defined (GUILE_I)
-#if HAVE_COMPLEX_DOUBLE
+#if defined HAVE_COMPLEX_DOUBLE
/* For an SCM object Z which is a complex number (ie. satisfies
SCM_COMPLEXP), return its value as a C level "complex double". */
static mpz_t z_negative_one;
\f
+
/* Clear the `mpz_t' embedded in bignum PTR. */
static void
finalize_bignum (GC_PTR ptr, GC_PTR data)
mpz_clear (SCM_I_BIG_MPZ (bignum));
}
+/* The next three functions (custom_libgmp_*) are passed to
+ mp_set_memory_functions (in GMP) so that memory used by the digits
+ themselves is known to the garbage collector. This is needed so
+ that GC will be run at appropriate times. Otherwise, a program which
+ creates many large bignums would malloc a huge amount of memory
+ before the GC runs. */
+static void *
+custom_gmp_malloc (size_t alloc_size)
+{
+ return scm_malloc (alloc_size);
+}
+
+static void *
+custom_gmp_realloc (void *old_ptr, size_t old_size, size_t new_size)
+{
+ return scm_realloc (old_ptr, new_size);
+}
+
+static void
+custom_gmp_free (void *ptr, size_t size)
+{
+ free (ptr);
+}
+
+
/* Return a new uninitialized bignum. */
static inline SCM
make_bignum (void)
we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
negatives as twos complement.
- In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
- following the hardware rounding mode, but applied to the absolute value
- of the mpz_t operand. This is not what we want so we put the high
- DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
- mpz_get_d is supposed to always truncate towards zero.
+ In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up
+ following the hardware rounding mode, but applied to the absolute
+ value of the mpz_t operand. This is not what we want so we put the
+ high DBL_MANT_DIG bits into a temporary. Starting with GMP 4.2
+ (released in March 2006) mpz_get_d now always truncates towards zero.
- ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
- is a slowdown. It'd be faster to pick out the relevant high bits with
- mpz_getlimbn if we could be bothered coding that, and if the new
- truncating gmp doesn't come out. */
+ ENHANCE-ME: The temporary init+clear to force the rounding in GMP
+ before 4.2 is a slowdown. It'd be faster to pick out the relevant
+ high bits with mpz_getlimbn. */
double
scm_i_big2dbl (SCM b)
#if 1
{
- /* Current GMP, eg. 4.1.3, force truncation towards zero */
+ /* For GMP earlier than 4.2, force truncation towards zero */
+
+ /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits,
+ _not_ the number of bits, so this code will break badly on a
+ system with non-binary doubles. */
+
mpz_t tmp;
if (bits > DBL_MANT_DIG)
{
}
}
#else
- /* Future GMP */
+ /* GMP 4.2 or later */
result = mpz_get_d (SCM_I_BIG_MPZ (b));
#endif
}
#undef FUNC_NAME
+int
+scm_is_exact (SCM val)
+{
+ return scm_is_true (scm_exact_p (val));
+}
SCM_PRIMITIVE_GENERIC (scm_inexact_p, "inexact?", 1, 0, 0,
(SCM x),
}
#undef FUNC_NAME
+int
+scm_is_inexact (SCM val)
+{
+ return scm_is_true (scm_inexact_p (val));
+}
SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0,
(SCM n),
"Return the quotient of the numbers @var{x} and @var{y}.")
#define FUNC_NAME s_scm_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_quotient);
- else
- {
- scm_t_inum z = xx / yy;
- if (SCM_LIKELY (SCM_FIXABLE (z)))
- return SCM_I_MAKINUM (z);
- else
- return scm_i_inum2big (z);
- }
- }
- else if (SCM_BIGP (y))
- {
- if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
- && (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
- SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
- }
- else if (SCM_BIGP (x))
+ if (SCM_LIKELY (scm_is_integer (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_quotient);
- else if (SCM_UNLIKELY (yy == 1))
- return x;
- else
- {
- SCM result = scm_i_mkbig ();
- if (yy < 0)
- {
- mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- - yy);
- mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
- }
- else
- mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- return scm_i_normbig (result);
- }
- }
- else if (SCM_BIGP (y))
- {
- SCM result = scm_i_mkbig ();
- mpz_tdiv_q (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_2 (x, y);
- return scm_i_normbig (result);
- }
+ if (SCM_LIKELY (scm_is_integer (y)))
+ return scm_truncate_quotient (x, y);
else
SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
}
"@end lisp")
#define FUNC_NAME s_scm_remainder
{
- if (SCM_LIKELY (SCM_I_INUMP (x)))
- {
- if (SCM_LIKELY (SCM_I_INUMP (y)))
- {
- scm_t_inum yy = SCM_I_INUM (y);
- if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_remainder);
- else
- {
- /* C99 specifies that "%" is the remainder corresponding to a
- quotient rounded towards zero, and that's also traditional
- for machine division, so z here should be well defined. */
- scm_t_inum z = SCM_I_INUM (x) % yy;
- return SCM_I_MAKINUM (z);
- }
- }
- else if (SCM_BIGP (y))
- {
- if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
- && (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
- SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
- }
- else if (SCM_BIGP (x))
+ if (SCM_LIKELY (scm_is_integer (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_remainder);
- else
- {
- SCM result = scm_i_mkbig ();
- if (yy < 0)
- yy = - yy;
- mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ(x), yy);
- scm_remember_upto_here_1 (x);
- return scm_i_normbig (result);
- }
- }
- else if (SCM_BIGP (y))
- {
- SCM result = scm_i_mkbig ();
- mpz_tdiv_r (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_2 (x, y);
- return scm_i_normbig (result);
- }
+ if (SCM_LIKELY (scm_is_integer (y)))
+ return scm_truncate_remainder (x, y);
else
SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
}
"@end lisp")
#define FUNC_NAME s_scm_modulo
{
- 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_modulo);
- else
- {
- /* C99 specifies that "%" is the remainder corresponding to a
- quotient rounded towards zero, and that's also traditional
- for machine division, so z here should be well defined. */
- scm_t_inum z = xx % yy;
- scm_t_inum result;
-
- if (yy < 0)
- {
- if (z > 0)
- result = z + yy;
- else
- result = z;
- }
- else
- {
- if (z < 0)
- result = z + yy;
- else
- result = z;
- }
- return SCM_I_MAKINUM (result);
- }
- }
- else if (SCM_BIGP (y))
- {
- int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
- {
- mpz_t z_x;
- SCM result;
-
- if (sgn_y < 0)
- {
- SCM pos_y = scm_i_clonebig (y, 0);
- /* do this after the last scm_op */
- mpz_init_set_si (z_x, xx);
- result = pos_y; /* re-use this bignum */
- mpz_mod (SCM_I_BIG_MPZ (result),
- z_x,
- SCM_I_BIG_MPZ (pos_y));
- scm_remember_upto_here_1 (pos_y);
- }
- else
- {
- result = scm_i_mkbig ();
- /* do this after the last scm_op */
- mpz_init_set_si (z_x, xx);
- mpz_mod (SCM_I_BIG_MPZ (result),
- z_x,
- SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_1 (y);
- }
-
- if ((sgn_y < 0) && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
- mpz_add (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (y),
- SCM_I_BIG_MPZ (result));
- scm_remember_upto_here_1 (y);
- /* and do this before the next one */
- mpz_clear (z_x);
- return scm_i_normbig (result);
- }
- }
- else
- SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
- }
- else if (SCM_BIGP (x))
+ if (SCM_LIKELY (scm_is_integer (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_modulo);
- else
- {
- SCM result = scm_i_mkbig ();
- mpz_mod_ui (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- (yy < 0) ? - yy : yy);
- scm_remember_upto_here_1 (x);
- if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
- mpz_sub_ui (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (result),
- - yy);
- return scm_i_normbig (result);
- }
- }
- else if (SCM_BIGP (y))
- {
- SCM result = scm_i_mkbig ();
- int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
- SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
- mpz_mod (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (pos_y));
-
- scm_remember_upto_here_1 (x);
- if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
- mpz_add (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (y),
- SCM_I_BIG_MPZ (result));
- scm_remember_upto_here_2 (y, pos_y);
- return scm_i_normbig (result);
- }
+ if (SCM_LIKELY (scm_is_integer (y)))
+ return scm_floor_remainder (x, y);
else
SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
}
}
#undef FUNC_NAME
-static SCM scm_i_inexact_euclidean_quotient (double x, double y);
-static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y);
+/* two_valued_wta_dispatch_2 is a version of SCM_WTA_DISPATCH_2 for
+ two-valued functions. It is called from primitive generics that take
+ two arguments and return two values, when the core procedure is
+ unable to handle the given argument types. If there are GOOPS
+ methods for this primitive generic, it dispatches to GOOPS and, if
+ successful, expects two values to be returned, which are placed in
+ *rp1 and *rp2. If there are no GOOPS methods, it throws a
+ wrong-type-arg exception.
+
+ FIXME: This obviously belongs somewhere else, but until we decide on
+ the right API, it is here as a static function, because it is needed
+ by the *_divide functions below.
+*/
+static void
+two_valued_wta_dispatch_2 (SCM gf, SCM a1, SCM a2, int pos,
+ const char *subr, SCM *rp1, SCM *rp2)
+{
+ if (SCM_UNPACK (gf))
+ scm_i_extract_values_2 (scm_call_generic_2 (gf, a1, a2), rp1, rp2);
+ else
+ scm_wrong_type_arg (subr, pos, (pos == SCM_ARG1) ? a1 : a2);
+}
+
+SCM_DEFINE (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the integer @var{q} such that\n"
+ "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "where @math{0 <= @var{r} < abs(@var{y})}.\n"
+ "@lisp\n"
+ "(euclidean-quotient 123 10) @result{} 12\n"
+ "(euclidean-quotient 123 -10) @result{} -12\n"
+ "(euclidean-quotient -123 10) @result{} -13\n"
+ "(euclidean-quotient -123 -10) @result{} 13\n"
+ "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
+ "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_euclidean_quotient
+{
+ if (scm_is_false (scm_negative_p (y)))
+ return scm_floor_quotient (x, y);
+ else
+ return scm_ceiling_quotient (x, y);
+}
+#undef FUNC_NAME
+
+SCM_DEFINE (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the real number @var{r} such that\n"
+ "@math{0 <= @var{r} < abs(@var{y})} and\n"
+ "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "for some integer @var{q}.\n"
+ "@lisp\n"
+ "(euclidean-remainder 123 10) @result{} 3\n"
+ "(euclidean-remainder 123 -10) @result{} 3\n"
+ "(euclidean-remainder -123 10) @result{} 7\n"
+ "(euclidean-remainder -123 -10) @result{} 7\n"
+ "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
+ "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_euclidean_remainder
+{
+ if (scm_is_false (scm_negative_p (y)))
+ return scm_floor_remainder (x, y);
+ else
+ return scm_ceiling_remainder (x, y);
+}
+#undef FUNC_NAME
+
+SCM_DEFINE (scm_i_euclidean_divide, "euclidean/", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the integer @var{q} and the real number @var{r}\n"
+ "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "and @math{0 <= @var{r} < abs(@var{y})}.\n"
+ "@lisp\n"
+ "(euclidean/ 123 10) @result{} 12 and 3\n"
+ "(euclidean/ 123 -10) @result{} -12 and 3\n"
+ "(euclidean/ -123 10) @result{} -13 and 7\n"
+ "(euclidean/ -123 -10) @result{} 13 and 7\n"
+ "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+ "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_i_euclidean_divide
+{
+ if (scm_is_false (scm_negative_p (y)))
+ return scm_i_floor_divide (x, y);
+ else
+ return scm_i_ceiling_divide (x, y);
+}
+#undef FUNC_NAME
+
+void
+scm_euclidean_divide (SCM x, SCM y, SCM *qp, SCM *rp)
+{
+ if (scm_is_false (scm_negative_p (y)))
+ return scm_floor_divide (x, y, qp, rp);
+ else
+ return scm_ceiling_divide (x, y, qp, rp);
+}
+
+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_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0,
(SCM x, SCM y),
- "Return the integer @var{q} such that\n"
- "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
- "where @math{0 <= @var{r} < abs(@var{y})}.\n"
+ "Return the floor of @math{@var{x} / @var{y}}.\n"
"@lisp\n"
- "(euclidean-quotient 123 10) @result{} 12\n"
- "(euclidean-quotient 123 -10) @result{} -12\n"
- "(euclidean-quotient -123 10) @result{} -13\n"
- "(euclidean-quotient -123 -10) @result{} 13\n"
- "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
- "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
+ "(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_euclidean_quotient
+#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);
- if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_quotient);
- else
+ scm_t_inum xx1 = xx;
+ scm_t_inum qq;
+ if (SCM_LIKELY (yy > 0))
{
- scm_t_inum xx = SCM_I_INUM (x);
- scm_t_inum qq = xx / yy;
- if (xx < qq * yy)
- {
- if (yy > 0)
- qq--;
- else
- qq++;
- }
- return SCM_I_MAKINUM (qq);
+ 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))
{
- if (SCM_I_INUM (x) >= 0)
- return SCM_INUM0;
+ 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 (- mpz_sgn (SCM_I_BIG_MPZ (y)));
+ return SCM_I_MAKINUM ((xx > 0) ? -1 : 0);
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_quotient
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_floor_quotient (xx, SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_quotient (x, y);
+ return scm_i_exact_rational_floor_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
- s_scm_euclidean_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2,
+ s_scm_floor_quotient);
}
else if (SCM_BIGP (x))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_quotient);
+ scm_num_overflow (s_scm_floor_quotient);
+ else if (SCM_UNLIKELY (yy == 1))
+ return x;
else
{
SCM q = scm_i_mkbig ();
mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy);
else
{
- mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy);
+ mpz_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);
else if (SCM_BIGP (y))
{
SCM q = scm_i_mkbig ();
- if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
- mpz_fdiv_q (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
- else
- mpz_cdiv_q (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
+ 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_euclidean_quotient
+ return scm_i_inexact_floor_quotient
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_quotient (x, y);
+ return scm_i_exact_rational_floor_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
- s_scm_euclidean_quotient);
+ 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_euclidean_quotient
+ return scm_i_inexact_floor_quotient
(SCM_REAL_VALUE (x), scm_to_double (y));
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
- s_scm_euclidean_quotient);
+ 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_euclidean_quotient
+ return scm_i_inexact_floor_quotient
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_floor_quotient (x, y);
else
- return scm_i_slow_exact_euclidean_quotient (x, y);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2,
+ s_scm_floor_quotient);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
- s_scm_euclidean_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1,
+ s_scm_floor_quotient);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_euclidean_quotient (double x, double y)
+scm_i_inexact_floor_quotient (double x, double y)
{
- if (SCM_LIKELY (y > 0))
- return scm_from_double (floor (x / y));
- else if (SCM_LIKELY (y < 0))
- return scm_from_double (ceil (x / y));
- else if (y == 0)
- scm_num_overflow (s_scm_euclidean_quotient); /* or return a NaN? */
+ if (SCM_UNLIKELY (y == 0))
+ scm_num_overflow (s_scm_floor_quotient); /* or return a NaN? */
else
- return scm_nan ();
+ return scm_from_double (floor (x / y));
}
-/* Compute exact euclidean_quotient the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
static SCM
-scm_i_slow_exact_euclidean_quotient (SCM x, SCM y)
+scm_i_exact_rational_floor_quotient (SCM x, SCM y)
{
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
- s_scm_euclidean_quotient);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
- s_scm_euclidean_quotient);
- else if (scm_is_true (scm_positive_p (y)))
- return scm_floor (scm_divide (x, y));
- else if (scm_is_true (scm_negative_p (y)))
- return scm_ceiling (scm_divide (x, y));
- else
- scm_num_overflow (s_scm_euclidean_quotient);
+ return scm_floor_quotient
+ (scm_product (scm_numerator (x), scm_denominator (y)),
+ scm_product (scm_numerator (y), scm_denominator (x)));
}
-static SCM scm_i_inexact_euclidean_remainder (double x, double y);
-static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y);
+static SCM scm_i_inexact_floor_remainder (double x, double y);
+static SCM scm_i_exact_rational_floor_remainder (SCM x, SCM y);
-SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
+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{0 <= @var{r} < abs(@var{y})} and\n"
"@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
- "for some integer @var{q}.\n"
+ "where @math{@var{q} = floor(@var{x} / @var{y})}.\n"
"@lisp\n"
- "(euclidean-remainder 123 10) @result{} 3\n"
- "(euclidean-remainder 123 -10) @result{} 3\n"
- "(euclidean-remainder -123 10) @result{} 7\n"
- "(euclidean-remainder -123 -10) @result{} 7\n"
- "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
- "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
+ "(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_euclidean_remainder
+#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_euclidean_remainder);
+ scm_num_overflow (s_scm_floor_remainder);
else
{
- scm_t_inum rr = SCM_I_INUM (x) % yy;
- if (rr >= 0)
- return SCM_I_MAKINUM (rr);
- else if (yy > 0)
- return SCM_I_MAKINUM (rr + yy);
+ scm_t_inum rr = xx % yy;
+ int needs_adjustment;
+
+ if (SCM_LIKELY (yy > 0))
+ needs_adjustment = (rr < 0);
else
- return SCM_I_MAKINUM (rr - yy);
+ needs_adjustment = (rr > 0);
+
+ if (needs_adjustment)
+ rr += yy;
+ return SCM_I_MAKINUM (rr);
}
}
else if (SCM_BIGP (y))
{
- scm_t_inum xx = SCM_I_INUM (x);
- if (xx >= 0)
- return x;
- else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ int sign = mpz_sgn (SCM_I_BIG_MPZ (y));
+ scm_remember_upto_here_1 (y);
+ if (sign > 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);
+ 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);
+ mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), xx);
scm_remember_upto_here_1 (y);
- mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r));
return scm_i_normbig (r);
}
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_remainder
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_floor_remainder (xx, SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_remainder (x, y);
+ return scm_i_exact_rational_floor_remainder (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
- s_scm_euclidean_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2,
+ s_scm_floor_remainder);
}
else if (SCM_BIGP (x))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_remainder);
+ scm_num_overflow (s_scm_floor_remainder);
else
{
scm_t_inum rr;
- if (yy < 0)
- yy = -yy;
- rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
+ 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_mod (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
+ 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_euclidean_remainder
+ return scm_i_inexact_floor_remainder
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_remainder (x, y);
+ return scm_i_exact_rational_floor_remainder (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
- s_scm_euclidean_remainder);
+ 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_euclidean_remainder
+ return scm_i_inexact_floor_remainder
(SCM_REAL_VALUE (x), scm_to_double (y));
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
- s_scm_euclidean_remainder);
+ 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_euclidean_remainder
+ return scm_i_inexact_floor_remainder
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_floor_remainder (x, y);
else
- return scm_i_slow_exact_euclidean_remainder (x, y);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2,
+ s_scm_floor_remainder);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
- s_scm_euclidean_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1,
+ s_scm_floor_remainder);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_euclidean_remainder (double x, double y)
+scm_i_inexact_floor_remainder (double x, double y)
{
- double q;
-
/* Although it would be more efficient to use fmod here, we can't
because it would in some cases produce results inconsistent with
- scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
- even close). In particular, when x is very close to a multiple of
- y, then r might be either 0.0 or abs(y)-epsilon, but those two
- cases must correspond to different choices of q. If r = 0.0 then q
- must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
- chooses one and remainder chooses the other, it would be bad. This
- problem was observed with x = 130.0 and y = 10/7. */
- if (SCM_LIKELY (y > 0))
- q = floor (x / y);
- else if (SCM_LIKELY (y < 0))
- q = ceil (x / y);
- else if (y == 0)
- scm_num_overflow (s_scm_euclidean_remainder); /* or return a NaN? */
+ 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_nan ();
- return scm_from_double (x - q * y);
+ return scm_from_double (x - y * floor (x / y));
}
-/* Compute exact euclidean_remainder the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
static SCM
-scm_i_slow_exact_euclidean_remainder (SCM x, SCM y)
+scm_i_exact_rational_floor_remainder (SCM x, SCM y)
{
- SCM q;
-
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
- s_scm_euclidean_remainder);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
- s_scm_euclidean_remainder);
- else if (scm_is_true (scm_positive_p (y)))
- q = scm_floor (scm_divide (x, y));
- else if (scm_is_true (scm_negative_p (y)))
- q = scm_ceiling (scm_divide (x, y));
- else
- scm_num_overflow (s_scm_euclidean_remainder);
- return scm_difference (x, scm_product (y, q));
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
+ SCM r1 = scm_floor_remainder (scm_product (scm_numerator (x), yd),
+ scm_product (scm_numerator (y), xd));
+ return scm_divide (r1, scm_product (xd, yd));
}
-static SCM scm_i_inexact_euclidean_divide (double x, double y);
-static SCM scm_i_slow_exact_euclidean_divide (SCM x, SCM y);
+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_euclidean_divide, "euclidean/", 2, 0, 0,
+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{0 <= @var{r} < abs(@var{y})}.\n"
+ "and @math{@var{q} = floor(@var{x} / @var{y})}.\n"
"@lisp\n"
- "(euclidean/ 123 10) @result{} 12 and 3\n"
- "(euclidean/ 123 -10) @result{} -12 and 3\n"
- "(euclidean/ -123 10) @result{} -13 and 7\n"
- "(euclidean/ -123 -10) @result{} 13 and 7\n"
- "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
- "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
+ "(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_euclidean_divide
+#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_euclidean_divide);
+ scm_num_overflow (s_scm_floor_divide);
else
{
- scm_t_inum xx = SCM_I_INUM (x);
scm_t_inum qq = xx / yy;
- scm_t_inum rr = xx - qq * yy;
- if (rr < 0)
+ 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)
{
- if (yy > 0)
- { rr += yy; qq--; }
- else
- { rr -= yy; qq++; }
+ rr += yy;
+ qq--;
}
- return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
- SCM_I_MAKINUM (rr)));
+
+ 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))
{
- scm_t_inum xx = SCM_I_INUM (x);
- if (xx >= 0)
- return scm_values (scm_list_2 (SCM_INUM0, x));
- else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ int sign = mpz_sgn (SCM_I_BIG_MPZ (y));
+ scm_remember_upto_here_1 (y);
+ if (sign > 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)));
+ 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);
+ 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)));
+ *qp = SCM_I_MAKINUM (-1);
+ *rp = scm_i_normbig (r);
}
+ return;
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_divide
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_floor_divide (xx, SCM_REAL_VALUE (y), qp, rp);
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_euclidean_divide (x, y);
+ return scm_i_exact_rational_floor_divide (x, y, qp, rp);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
- s_scm_euclidean_divide);
+ 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))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_euclidean_divide);
+ scm_num_overflow (s_scm_floor_divide);
else
{
SCM q = scm_i_mkbig ();
SCM_I_BIG_MPZ (x), yy);
else
{
- mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+ 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)));
+ *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 ();
- if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
- mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- else
- mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ 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)));
+ *qp = scm_i_normbig (q);
+ *rp = scm_i_normbig (r);
+ return;
}
else if (SCM_REALP (y))
- return scm_i_inexact_euclidean_divide
- (scm_i_big2dbl (x), SCM_REAL_VALUE (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_slow_exact_euclidean_divide (x, y);
+ return scm_i_exact_rational_floor_divide (x, y, qp, rp);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
- s_scm_euclidean_divide);
+ 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_euclidean_divide
- (SCM_REAL_VALUE (x), scm_to_double (y));
+ return scm_i_inexact_floor_divide
+ (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp);
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
- s_scm_euclidean_divide);
+ 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_euclidean_divide
- (scm_i_fraction2double (x), SCM_REAL_VALUE (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 scm_i_slow_exact_euclidean_divide (x, y);
+ return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2,
+ s_scm_floor_divide, qp, rp);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
- s_scm_euclidean_divide);
+ return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG1,
+ s_scm_floor_divide, qp, rp);
}
-#undef FUNC_NAME
-static SCM
-scm_i_inexact_euclidean_divide (double x, double y)
+static void
+scm_i_inexact_floor_divide (double x, double y, SCM *qp, SCM *rp)
{
- double q, r;
-
- if (SCM_LIKELY (y > 0))
- q = floor (x / y);
- else if (SCM_LIKELY (y < 0))
- q = ceil (x / y);
- else if (y == 0)
- scm_num_overflow (s_scm_euclidean_divide); /* or return a NaN? */
+ if (SCM_UNLIKELY (y == 0))
+ scm_num_overflow (s_scm_floor_divide); /* or return a NaN? */
else
- q = guile_NaN;
- r = x - q * y;
- return scm_values (scm_list_2 (scm_from_double (q),
- scm_from_double (r)));
+ {
+ double q = floor (x / y);
+ double r = x - q * y;
+ *qp = scm_from_double (q);
+ *rp = scm_from_double (r);
+ }
}
-/* Compute exact euclidean quotient and remainder the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
-static SCM
-scm_i_slow_exact_euclidean_divide (SCM x, SCM y)
+static void
+scm_i_exact_rational_floor_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{
- SCM q, r;
+ SCM r1;
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
- s_scm_euclidean_divide);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
- s_scm_euclidean_divide);
- else if (scm_is_true (scm_positive_p (y)))
- q = scm_floor (scm_divide (x, y));
- else if (scm_is_true (scm_negative_p (y)))
- q = scm_ceiling (scm_divide (x, y));
- else
- scm_num_overflow (s_scm_euclidean_divide);
- r = scm_difference (x, scm_product (q, y));
- return scm_values (scm_list_2 (q, r));
+ 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_centered_quotient (double x, double y);
-static SCM scm_i_bigint_centered_quotient (SCM x, SCM y);
-static SCM scm_i_slow_exact_centered_quotient (SCM x, SCM y);
+static SCM scm_i_inexact_ceiling_quotient (double x, double y);
+static SCM scm_i_exact_rational_ceiling_quotient (SCM x, SCM y);
-SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
+SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0,
(SCM x, SCM y),
- "Return the integer @var{q} such that\n"
- "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
- "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+ "Return the ceiling of @math{@var{x} / @var{y}}.\n"
"@lisp\n"
- "(centered-quotient 123 10) @result{} 12\n"
- "(centered-quotient 123 -10) @result{} -12\n"
- "(centered-quotient -123 10) @result{} -12\n"
- "(centered-quotient -123 -10) @result{} 12\n"
- "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
- "(centered-quotient 16/3 -10/7) @result{} -4\n"
+ "(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_centered_quotient
+#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_centered_quotient);
+ scm_num_overflow (s_scm_ceiling_quotient);
else
{
- scm_t_inum xx = SCM_I_INUM (x);
- scm_t_inum qq = xx / yy;
- scm_t_inum rr = xx - qq * yy;
- if (SCM_LIKELY (xx > 0))
+ scm_t_inum xx1 = xx;
+ scm_t_inum qq;
+ if (SCM_LIKELY (yy > 0))
{
- if (SCM_LIKELY (yy > 0))
- {
- if (rr >= (yy + 1) / 2)
- qq++;
- }
- else
- {
- if (rr >= (1 - yy) / 2)
- qq--;
- }
+ if (SCM_LIKELY (xx >= 0))
+ xx1 = xx + yy - 1;
}
+ else if (xx < 0)
+ xx1 = xx + yy + 1;
+ qq = xx1 / yy;
+ if (SCM_LIKELY (SCM_FIXABLE (qq)))
+ return SCM_I_MAKINUM (qq);
else
- {
- if (SCM_LIKELY (yy > 0))
- {
- if (rr < -yy / 2)
- qq--;
- }
- else
- {
- if (rr < yy / 2)
- qq++;
- }
- }
- return SCM_I_MAKINUM (qq);
+ 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_centered_quotient */
- return scm_i_bigint_centered_quotient
- (scm_i_long2big (SCM_I_INUM (x)), 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_centered_quotient
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_ceiling_quotient (xx, SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_quotient (x, y);
+ return scm_i_exact_rational_ceiling_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
- s_scm_centered_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2,
+ s_scm_ceiling_quotient);
}
else if (SCM_BIGP (x))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_centered_quotient);
+ scm_num_overflow (s_scm_ceiling_quotient);
+ else if (SCM_UNLIKELY (yy == 1))
+ return x;
else
{
SCM q = scm_i_mkbig ();
- scm_t_inum rr;
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
if (yy > 0)
- {
- rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- if (rr < -yy / 2)
- mpz_sub_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- }
+ mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy);
else
{
- rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (x), -yy);
- scm_remember_upto_here_1 (x);
+ mpz_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));
- if (rr < yy / 2)
- mpz_add_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
}
+ scm_remember_upto_here_1 (x);
return scm_i_normbig (q);
}
}
else if (SCM_BIGP (y))
- return scm_i_bigint_centered_quotient (x, 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_centered_quotient
+ return scm_i_inexact_ceiling_quotient
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_quotient (x, y);
+ return scm_i_exact_rational_ceiling_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
- s_scm_centered_quotient);
+ 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_centered_quotient
+ return scm_i_inexact_ceiling_quotient
(SCM_REAL_VALUE (x), scm_to_double (y));
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
- s_scm_centered_quotient);
+ 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_centered_quotient
+ return scm_i_inexact_ceiling_quotient
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_ceiling_quotient (x, y);
else
- return scm_i_slow_exact_centered_quotient (x, y);
+ SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2,
+ s_scm_ceiling_quotient);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
- s_scm_centered_quotient);
+ SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1,
+ s_scm_ceiling_quotient);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_centered_quotient (double x, double y)
-{
- if (SCM_LIKELY (y > 0))
- return scm_from_double (floor (x/y + 0.5));
- else if (SCM_LIKELY (y < 0))
- return scm_from_double (ceil (x/y - 0.5));
- else if (y == 0)
- scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */
- else
- return scm_nan ();
-}
-
-/* Assumes that both x and y are bigints, though
- x might be able to fit into a fixnum. */
-static SCM
-scm_i_bigint_centered_quotient (SCM x, SCM y)
+scm_i_inexact_ceiling_quotient (double x, double y)
{
- SCM q, r, min_r;
-
- /* Note that x might be small enough to fit into a
- fixnum, so we must not let it escape into the wild */
- q = scm_i_mkbig ();
- r = scm_i_mkbig ();
-
- /* min_r will eventually become -abs(y)/2 */
- min_r = scm_i_mkbig ();
- mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
- SCM_I_BIG_MPZ (y), 1);
-
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
- if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
- {
- mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_2 (x, y);
- mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- mpz_sub_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- }
+ if (SCM_UNLIKELY (y == 0))
+ scm_num_overflow (s_scm_ceiling_quotient); /* or return a NaN? */
else
- {
- mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_2 (x, y);
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- mpz_add_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- }
- scm_remember_upto_here_2 (r, min_r);
- return scm_i_normbig (q);
+ return scm_from_double (ceil (x / y));
}
-/* Compute exact centered quotient the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
static SCM
-scm_i_slow_exact_centered_quotient (SCM x, SCM y)
+scm_i_exact_rational_ceiling_quotient (SCM x, SCM y)
{
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
- s_scm_centered_quotient);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
- s_scm_centered_quotient);
- else if (scm_is_true (scm_positive_p (y)))
- return scm_floor (scm_sum (scm_divide (x, y),
- exactly_one_half));
- else if (scm_is_true (scm_negative_p (y)))
- return scm_ceiling (scm_difference (scm_divide (x, y),
- exactly_one_half));
- else
- scm_num_overflow (s_scm_centered_quotient);
+ return scm_ceiling_quotient
+ (scm_product (scm_numerator (x), scm_denominator (y)),
+ scm_product (scm_numerator (y), scm_denominator (x)));
}
-static SCM scm_i_inexact_centered_remainder (double x, double y);
-static SCM scm_i_bigint_centered_remainder (SCM x, SCM y);
-static SCM scm_i_slow_exact_centered_remainder (SCM x, SCM y);
+static SCM scm_i_inexact_ceiling_remainder (double x, double y);
+static SCM scm_i_exact_rational_ceiling_remainder (SCM x, SCM y);
-SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
+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{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
- "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
- "for some integer @var{q}.\n"
+ "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "where @math{@var{q} = ceiling(@var{x} / @var{y})}.\n"
"@lisp\n"
- "(centered-remainder 123 10) @result{} 3\n"
- "(centered-remainder 123 -10) @result{} 3\n"
- "(centered-remainder -123 10) @result{} -3\n"
- "(centered-remainder -123 -10) @result{} -3\n"
- "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
- "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
+ "(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_centered_remainder
+#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_centered_remainder);
+ scm_num_overflow (s_scm_ceiling_remainder);
else
{
- scm_t_inum xx = SCM_I_INUM (x);
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))
{
- if (SCM_LIKELY (yy > 0))
- {
- if (rr >= (yy + 1) / 2)
- rr -= yy;
- }
- else
- {
- if (rr >= (1 - yy) / 2)
- rr += yy;
- }
+ 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 (trunc (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 * trunc (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 = trunc (x / y);
+ double 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);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the integer @var{q} such that\n"
+ "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
+ "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+ "@lisp\n"
+ "(centered-quotient 123 10) @result{} 12\n"
+ "(centered-quotient 123 -10) @result{} -12\n"
+ "(centered-quotient -123 10) @result{} -12\n"
+ "(centered-quotient -123 -10) @result{} 12\n"
+ "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
+ "(centered-quotient 16/3 -10/7) @result{} -4\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_centered_quotient
+{
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
+ {
+ 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_centered_quotient);
+ else
+ {
+ scm_t_inum qq = xx / yy;
+ scm_t_inum rr = xx % yy;
+ if (SCM_LIKELY (xx > 0))
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr >= (yy + 1) / 2)
+ qq++;
+ }
+ else
+ {
+ if (rr >= (1 - yy) / 2)
+ qq--;
+ }
+ }
+ else
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr < -yy / 2)
+ qq--;
+ }
+ else
+ {
+ if (rr < yy / 2)
+ qq++;
+ }
+ }
+ 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_centered_quotient */
+ return scm_i_bigint_centered_quotient (scm_i_long2big (xx), y);
+ }
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_quotient (xx, SCM_REAL_VALUE (y));
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_quotient (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+ s_scm_centered_quotient);
+ }
+ else if (SCM_BIGP (x))
+ {
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
+ {
+ scm_t_inum yy = SCM_I_INUM (y);
+ if (SCM_UNLIKELY (yy == 0))
+ scm_num_overflow (s_scm_centered_quotient);
+ else if (SCM_UNLIKELY (yy == 1))
+ return x;
+ else
+ {
+ SCM q = scm_i_mkbig ();
+ scm_t_inum rr;
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (yy > 0)
+ {
+ rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (x), yy);
+ scm_remember_upto_here_1 (x);
+ if (rr < -yy / 2)
+ mpz_sub_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ }
+ else
+ {
+ rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (x), -yy);
+ scm_remember_upto_here_1 (x);
+ mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+ if (rr < yy / 2)
+ mpz_add_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ }
+ return scm_i_normbig (q);
+ }
+ }
+ else if (SCM_BIGP (y))
+ return scm_i_bigint_centered_quotient (x, y);
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_quotient
+ (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_quotient (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+ s_scm_centered_quotient);
+ }
+ else if (SCM_REALP (x))
+ {
+ if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+ SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_inexact_centered_quotient
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+ s_scm_centered_quotient);
+ }
+ else if (SCM_FRACTIONP (x))
+ {
+ if (SCM_REALP (y))
+ return scm_i_inexact_centered_quotient
+ (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_quotient (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+ s_scm_centered_quotient);
+ }
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
+ s_scm_centered_quotient);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_quotient (double x, double y)
+{
+ if (SCM_LIKELY (y > 0))
+ return scm_from_double (floor (x/y + 0.5));
+ else if (SCM_LIKELY (y < 0))
+ return scm_from_double (ceil (x/y - 0.5));
+ else if (y == 0)
+ scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */
+ else
+ return scm_nan ();
+}
+
+/* Assumes that both x and y are bigints, though
+ x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_quotient (SCM x, SCM y)
+{
+ SCM q, r, min_r;
+
+ /* Note that x might be small enough to fit into a
+ fixnum, so we must not let it escape into the wild */
+ q = scm_i_mkbig ();
+ r = scm_i_mkbig ();
+
+ /* min_r will eventually become -abs(y)/2 */
+ min_r = scm_i_mkbig ();
+ mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+ SCM_I_BIG_MPZ (y), 1);
+
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ {
+ mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ scm_remember_upto_here_2 (x, y);
+ mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ mpz_sub_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ }
+ else
+ {
+ mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ scm_remember_upto_here_2 (x, y);
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ mpz_add_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ }
+ scm_remember_upto_here_2 (r, min_r);
+ return scm_i_normbig (q);
+}
+
+static SCM
+scm_i_exact_rational_centered_quotient (SCM x, SCM y)
+{
+ return scm_centered_quotient
+ (scm_product (scm_numerator (x), scm_denominator (y)),
+ scm_product (scm_numerator (y), scm_denominator (x)));
+}
+
+static SCM scm_i_inexact_centered_remainder (double x, double y);
+static SCM scm_i_bigint_centered_remainder (SCM x, SCM y);
+static SCM scm_i_exact_rational_centered_remainder (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the real number @var{r} such that\n"
+ "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
+ "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "for some integer @var{q}.\n"
+ "@lisp\n"
+ "(centered-remainder 123 10) @result{} 3\n"
+ "(centered-remainder 123 -10) @result{} 3\n"
+ "(centered-remainder -123 10) @result{} -3\n"
+ "(centered-remainder -123 -10) @result{} -3\n"
+ "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
+ "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_centered_remainder
+{
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
+ {
+ 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_centered_remainder);
+ else
+ {
+ scm_t_inum rr = xx % yy;
+ if (SCM_LIKELY (xx > 0))
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr >= (yy + 1) / 2)
+ rr -= yy;
+ }
+ else
+ {
+ if (rr >= (1 - yy) / 2)
+ rr += yy;
+ }
+ }
+ else
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr < -yy / 2)
+ rr += yy;
+ }
+ else
+ {
+ if (rr < yy / 2)
+ rr -= yy;
+ }
+ }
+ return SCM_I_MAKINUM (rr);
+ }
+ }
+ else if (SCM_BIGP (y))
+ {
+ /* Pass a denormalized bignum version of x (even though it
+ can fit in a fixnum) to scm_i_bigint_centered_remainder */
+ return scm_i_bigint_centered_remainder (scm_i_long2big (xx), y);
+ }
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_remainder (xx, SCM_REAL_VALUE (y));
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_remainder (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+ s_scm_centered_remainder);
+ }
+ else if (SCM_BIGP (x))
+ {
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
+ {
+ scm_t_inum yy = SCM_I_INUM (y);
+ if (SCM_UNLIKELY (yy == 0))
+ scm_num_overflow (s_scm_centered_remainder);
+ else
+ {
+ scm_t_inum rr;
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (yy > 0)
+ {
+ rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy);
+ scm_remember_upto_here_1 (x);
+ if (rr < -yy / 2)
+ rr += yy;
+ }
+ else
+ {
+ rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy);
+ scm_remember_upto_here_1 (x);
+ if (rr < yy / 2)
+ rr -= yy;
+ }
+ return SCM_I_MAKINUM (rr);
+ }
+ }
+ else if (SCM_BIGP (y))
+ return scm_i_bigint_centered_remainder (x, y);
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_remainder
+ (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_remainder (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+ s_scm_centered_remainder);
+ }
+ else if (SCM_REALP (x))
+ {
+ if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+ SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_inexact_centered_remainder
+ (SCM_REAL_VALUE (x), scm_to_double (y));
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+ s_scm_centered_remainder);
+ }
+ else if (SCM_FRACTIONP (x))
+ {
+ if (SCM_REALP (y))
+ return scm_i_inexact_centered_remainder
+ (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_remainder (x, y);
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+ s_scm_centered_remainder);
+ }
+ else
+ SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
+ s_scm_centered_remainder);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_remainder (double x, double y)
+{
+ double q;
+
+ /* Although it would be more efficient to use fmod here, we can't
+ because it would in some cases produce results inconsistent with
+ scm_i_inexact_centered_quotient, such that x != r + q * y (not even
+ close). In particular, when x-y/2 is very close to a multiple of
+ y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
+ two cases must correspond to different choices of q. If quotient
+ chooses one and remainder chooses the other, it would be bad. */
+ if (SCM_LIKELY (y > 0))
+ q = floor (x/y + 0.5);
+ else if (SCM_LIKELY (y < 0))
+ q = ceil (x/y - 0.5);
+ else if (y == 0)
+ scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */
+ else
+ return scm_nan ();
+ return scm_from_double (x - q * y);
+}
+
+/* Assumes that both x and y are bigints, though
+ x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_remainder (SCM x, SCM y)
+{
+ SCM r, min_r;
+
+ /* Note that x might be small enough to fit into a
+ fixnum, so we must not let it escape into the wild */
+ r = scm_i_mkbig ();
+
+ /* min_r will eventually become -abs(y)/2 */
+ min_r = scm_i_mkbig ();
+ mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+ SCM_I_BIG_MPZ (y), 1);
+
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ {
+ mpz_cdiv_r (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ mpz_add (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (y));
+ }
+ else
+ {
+ mpz_fdiv_r (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ mpz_sub (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (y));
+ }
+ scm_remember_upto_here_2 (x, y);
+ return scm_i_normbig (r);
+}
+
+static SCM
+scm_i_exact_rational_centered_remainder (SCM x, SCM y)
+{
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
+ SCM r1 = scm_centered_remainder (scm_product (scm_numerator (x), yd),
+ scm_product (scm_numerator (y), xd));
+ return scm_divide (r1, scm_product (xd, yd));
+}
+
+
+static void scm_i_inexact_centered_divide (double x, double y,
+ SCM *qp, SCM *rp);
+static void scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp);
+static void scm_i_exact_rational_centered_divide (SCM x, SCM y,
+ SCM *qp, SCM *rp);
+
+SCM_PRIMITIVE_GENERIC (scm_i_centered_divide, "centered/", 2, 0, 0,
+ (SCM x, SCM y),
+ "Return the integer @var{q} and the real number @var{r}\n"
+ "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+ "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+ "@lisp\n"
+ "(centered/ 123 10) @result{} 12 and 3\n"
+ "(centered/ 123 -10) @result{} -12 and 3\n"
+ "(centered/ -123 10) @result{} -12 and -3\n"
+ "(centered/ -123 -10) @result{} 12 and -3\n"
+ "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+ "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_i_centered_divide
+{
+ SCM q, r;
+
+ scm_centered_divide(x, y, &q, &r);
+ return scm_values (scm_list_2 (q, r));
+}
+#undef FUNC_NAME
+
+#define s_scm_centered_divide s_scm_i_centered_divide
+#define g_scm_centered_divide g_scm_i_centered_divide
+
+void
+scm_centered_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_centered_divide);
+ else
+ {
+ scm_t_inum qq = xx / yy;
+ scm_t_inum rr = xx % yy;
+ if (SCM_LIKELY (xx > 0))
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr >= (yy + 1) / 2)
+ { qq++; rr -= yy; }
+ }
+ else
+ {
+ if (rr >= (1 - yy) / 2)
+ { qq--; rr += yy; }
+ }
+ }
+ else
+ {
+ if (SCM_LIKELY (yy > 0))
+ {
+ if (rr < -yy / 2)
+ { qq--; rr += yy; }
+ }
+ else
+ {
+ if (rr < yy / 2)
+ { qq++; rr -= yy; }
+ }
+ }
+ 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_centered_divide */
+ return scm_i_bigint_centered_divide (scm_i_long2big (xx), y, qp, rp);
+ }
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y), qp, rp);
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_divide (x, y, qp, rp);
+ else
+ return two_valued_wta_dispatch_2
+ (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_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_centered_divide);
+ else
+ {
+ SCM q = scm_i_mkbig ();
+ scm_t_inum rr;
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (yy > 0)
+ {
+ rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (x), yy);
+ scm_remember_upto_here_1 (x);
+ if (rr < -yy / 2)
+ {
+ mpz_sub_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ rr += yy;
+ }
+ }
+ else
+ {
+ rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (x), -yy);
+ scm_remember_upto_here_1 (x);
+ mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+ if (rr < yy / 2)
+ {
+ mpz_add_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ rr -= yy;
+ }
+ }
+ *qp = scm_i_normbig (q);
+ *rp = SCM_I_MAKINUM (rr);
+ }
+ return;
+ }
+ else if (SCM_BIGP (y))
+ return scm_i_bigint_centered_divide (x, y, qp, rp);
+ else if (SCM_REALP (y))
+ return scm_i_inexact_centered_divide
+ (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp);
+ else if (SCM_FRACTIONP (y))
+ return scm_i_exact_rational_centered_divide (x, y, qp, rp);
+ else
+ return two_valued_wta_dispatch_2
+ (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_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_centered_divide
+ (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp);
+ else
+ return two_valued_wta_dispatch_2
+ (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide, qp, rp);
+ }
+ else if (SCM_FRACTIONP (x))
+ {
+ if (SCM_REALP (y))
+ return scm_i_inexact_centered_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_centered_divide (x, y, qp, rp);
+ else
+ return two_valued_wta_dispatch_2
+ (g_scm_centered_divide, x, y, SCM_ARG2,
+ s_scm_centered_divide, qp, rp);
+ }
+ else
+ return two_valued_wta_dispatch_2 (g_scm_centered_divide, x, y, SCM_ARG1,
+ s_scm_centered_divide, qp, rp);
+}
+
+static void
+scm_i_inexact_centered_divide (double x, double y, SCM *qp, SCM *rp)
+{
+ double q, r;
+
+ if (SCM_LIKELY (y > 0))
+ q = floor (x/y + 0.5);
+ else if (SCM_LIKELY (y < 0))
+ q = ceil (x/y - 0.5);
+ else if (y == 0)
+ scm_num_overflow (s_scm_centered_divide); /* or return a NaN? */
+ else
+ q = guile_NaN;
+ 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_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
+{
+ SCM q, r, min_r;
+
+ /* Note that x might be small enough to fit into a
+ fixnum, so we must not let it escape into the wild */
+ q = scm_i_mkbig ();
+ r = scm_i_mkbig ();
+
+ /* min_r will eventually become -abs(y/2) */
+ min_r = scm_i_mkbig ();
+ mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+ SCM_I_BIG_MPZ (y), 1);
+
+ /* Arrange for rr to initially be non-positive,
+ because that simplifies the test to see
+ if it is within the needed bounds. */
+ if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ {
+ mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ {
+ mpz_sub_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ mpz_add (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (y));
+ }
+ }
+ else
+ {
+ mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+ if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+ {
+ mpz_add_ui (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (q), 1);
+ mpz_sub (SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (r),
+ SCM_I_BIG_MPZ (y));
+ }
+ }
+ scm_remember_upto_here_2 (x, y);
+ *qp = scm_i_normbig (q);
+ *rp = scm_i_normbig (r);
+}
+
+static void
+scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
+{
+ SCM r1;
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
+
+ scm_centered_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_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 (SCM_LIKELY (yy > 0))
- {
- if (rr < -yy / 2)
- rr += yy;
- }
- else
- {
- if (rr < yy / 2)
- rr -= yy;
- }
+ if (r2 > ay)
+ qq++;
+ else if (r2 < -ay)
+ qq--;
}
- return SCM_I_MAKINUM (rr);
+ 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_centered_remainder */
- return scm_i_bigint_centered_remainder
- (scm_i_long2big (SCM_I_INUM (x)), y);
+ 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_centered_remainder
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_round_quotient (xx, SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_remainder (x, y);
+ return scm_i_exact_rational_round_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
- s_scm_centered_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2,
+ s_scm_round_quotient);
}
else if (SCM_BIGP (x))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_centered_remainder);
+ 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;
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
+ int needs_adjustment;
+
if (yy > 0)
{
- rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- if (rr < -yy / 2)
- rr += yy;
+ 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_ui (SCM_I_BIG_MPZ (x), -yy);
- scm_remember_upto_here_1 (x);
- if (rr < yy / 2)
- rr -= yy;
+ 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);
}
- return SCM_I_MAKINUM (rr);
+ 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_centered_remainder (x, y);
+ return scm_i_bigint_round_quotient (x, y);
else if (SCM_REALP (y))
- return scm_i_inexact_centered_remainder
+ return scm_i_inexact_round_quotient
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_remainder (x, y);
+ return scm_i_exact_rational_round_quotient (x, y);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
- s_scm_centered_remainder);
+ 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_centered_remainder
+ return scm_i_inexact_round_quotient
(SCM_REAL_VALUE (x), scm_to_double (y));
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
- s_scm_centered_remainder);
+ 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_centered_remainder
+ return scm_i_inexact_round_quotient
(scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+ else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))
+ return scm_i_exact_rational_round_quotient (x, y);
else
- return scm_i_slow_exact_centered_remainder (x, y);
+ SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2,
+ s_scm_round_quotient);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
- s_scm_centered_remainder);
+ SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1,
+ s_scm_round_quotient);
}
#undef FUNC_NAME
static SCM
-scm_i_inexact_centered_remainder (double x, double y)
+scm_i_inexact_round_quotient (double x, double y)
{
- double q;
-
- /* Although it would be more efficient to use fmod here, we can't
- because it would in some cases produce results inconsistent with
- scm_i_inexact_centered_quotient, such that x != r + q * y (not even
- close). In particular, when x-y/2 is very close to a multiple of
- y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
- two cases must correspond to different choices of q. If quotient
- chooses one and remainder chooses the other, it would be bad. */
- if (SCM_LIKELY (y > 0))
- q = floor (x/y + 0.5);
- else if (SCM_LIKELY (y < 0))
- q = ceil (x/y - 0.5);
- else if (y == 0)
- scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */
+ if (SCM_UNLIKELY (y == 0))
+ scm_num_overflow (s_scm_round_quotient); /* or return a NaN? */
else
- return scm_nan ();
- return scm_from_double (x - q * y);
+ 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_centered_remainder (SCM x, SCM y)
+scm_i_bigint_round_quotient (SCM x, SCM y)
{
- SCM r, min_r;
+ 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 ();
- /* min_r will eventually become -abs(y)/2 */
- min_r = scm_i_mkbig ();
- mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
- SCM_I_BIG_MPZ (y), 1);
+ 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);
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
- if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+ 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)))
{
- mpz_cdiv_r (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- mpz_add (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (y));
+ 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
{
- mpz_fdiv_r (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- mpz_sub (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (y));
+ double q = scm_c_round (x / y);
+ return scm_from_double (x - q * y);
}
- scm_remember_upto_here_2 (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_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 centered_remainder the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
static SCM
-scm_i_slow_exact_centered_remainder (SCM x, SCM y)
+scm_i_exact_rational_round_remainder (SCM x, SCM y)
{
- SCM q;
-
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
- s_scm_centered_remainder);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
- s_scm_centered_remainder);
- else if (scm_is_true (scm_positive_p (y)))
- q = scm_floor (scm_sum (scm_divide (x, y), exactly_one_half));
- else if (scm_is_true (scm_negative_p (y)))
- q = scm_ceiling (scm_difference (scm_divide (x, y), exactly_one_half));
- else
- scm_num_overflow (s_scm_centered_remainder);
- return scm_difference (x, scm_product (y, q));
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
+ SCM r1 = scm_round_remainder (scm_product (scm_numerator (x), yd),
+ scm_product (scm_numerator (y), xd));
+ return scm_divide (r1, scm_product (xd, yd));
}
-static SCM scm_i_inexact_centered_divide (double x, double y);
-static SCM scm_i_bigint_centered_divide (SCM x, SCM y);
-static SCM scm_i_slow_exact_centered_divide (SCM x, SCM y);
+static 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_centered_divide, "centered/", 2, 0, 0,
+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 @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\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"
- "(centered/ 123 10) @result{} 12 and 3\n"
- "(centered/ 123 -10) @result{} -12 and 3\n"
- "(centered/ -123 10) @result{} -12 and -3\n"
- "(centered/ -123 -10) @result{} 12 and -3\n"
- "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
- "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
+ "(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_centered_divide
+#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_centered_divide);
+ scm_num_overflow (s_scm_round_divide);
else
{
- scm_t_inum xx = SCM_I_INUM (x);
scm_t_inum qq = xx / yy;
- scm_t_inum rr = xx - qq * yy;
- if (SCM_LIKELY (xx > 0))
+ scm_t_inum rr = xx % yy;
+ scm_t_inum ay = yy;
+ scm_t_inum r2 = 2 * rr;
+
+ if (SCM_LIKELY (yy < 0))
{
- if (SCM_LIKELY (yy > 0))
- {
- if (rr >= (yy + 1) / 2)
- { qq++; rr -= yy; }
- }
- else
- {
- if (rr >= (1 - yy) / 2)
- { qq--; rr += yy; }
- }
+ ay = -ay;
+ r2 = -r2;
+ }
+
+ if (qq & 1L)
+ {
+ if (r2 >= ay)
+ { qq++; rr -= yy; }
+ else if (r2 <= -ay)
+ { qq--; rr += yy; }
}
else
{
- if (SCM_LIKELY (yy > 0))
- {
- if (rr < -yy / 2)
- { qq--; rr += yy; }
- }
- else
- {
- if (rr < yy / 2)
- { qq++; rr -= yy; }
- }
+ if (r2 > ay)
+ { qq++; rr -= yy; }
+ else if (r2 < -ay)
+ { qq--; rr += yy; }
}
- return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
- SCM_I_MAKINUM (rr)));
+ 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_centered_divide */
- return scm_i_bigint_centered_divide
- (scm_i_long2big (SCM_I_INUM (x)), y);
+ 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_centered_divide
- (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+ return scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y), qp, rp);
else if (SCM_FRACTIONP (y))
- return scm_i_slow_exact_centered_divide (x, y);
+ return scm_i_exact_rational_round_divide (x, y, qp, rp);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
- s_scm_centered_divide);
+ 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))
{
{
scm_t_inum yy = SCM_I_INUM (y);
if (SCM_UNLIKELY (yy == 0))
- scm_num_overflow (s_scm_centered_divide);
+ scm_num_overflow (s_scm_round_divide);
else
{
SCM q = scm_i_mkbig ();
scm_t_inum rr;
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
+ int needs_adjustment;
+
if (yy > 0)
{
- rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- if (rr < -yy / 2)
- {
- mpz_sub_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- rr += yy;
- }
+ 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);
- scm_remember_upto_here_1 (x);
mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
- if (rr < yy / 2)
- {
- mpz_add_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- rr -= yy;
- }
+ 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)));
+ *qp = scm_i_normbig (q);
+ *rp = SCM_I_MAKINUM (rr);
}
+ return;
}
else if (SCM_BIGP (y))
- return scm_i_bigint_centered_divide (x, y);
+ return scm_i_bigint_round_divide (x, y, qp, rp);
else if (SCM_REALP (y))
- return scm_i_inexact_centered_divide
- (scm_i_big2dbl (x), SCM_REAL_VALUE (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_slow_exact_centered_divide (x, y);
+ return scm_i_exact_rational_round_divide (x, y, qp, rp);
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
- s_scm_centered_divide);
+ 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_centered_divide
- (SCM_REAL_VALUE (x), scm_to_double (y));
- else
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
- s_scm_centered_divide);
+ 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_centered_divide
- (scm_i_fraction2double (x), SCM_REAL_VALUE (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 scm_i_slow_exact_centered_divide (x, y);
+ return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2,
+ s_scm_round_divide, qp, rp);
}
else
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
- s_scm_centered_divide);
+ return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG1,
+ s_scm_round_divide, qp, rp);
}
-#undef FUNC_NAME
-static SCM
-scm_i_inexact_centered_divide (double x, double y)
+static void
+scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp)
{
- double q, r;
-
- if (SCM_LIKELY (y > 0))
- q = floor (x/y + 0.5);
- else if (SCM_LIKELY (y < 0))
- q = ceil (x/y - 0.5);
- else if (y == 0)
- scm_num_overflow (s_scm_centered_divide); /* or return a NaN? */
+ if (SCM_UNLIKELY (y == 0))
+ scm_num_overflow (s_scm_round_divide); /* or return a NaN? */
else
- q = guile_NaN;
- r = x - q * y;
- return scm_values (scm_list_2 (scm_from_double (q),
- scm_from_double (r)));
+ {
+ 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 SCM
-scm_i_bigint_centered_divide (SCM x, SCM y)
+static void
+scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{
- SCM q, r, min_r;
+ 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 ();
- /* min_r will eventually become -abs(y/2) */
- min_r = scm_i_mkbig ();
- mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
- SCM_I_BIG_MPZ (y), 1);
+ 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 */
- /* Arrange for rr to initially be non-positive,
- because that simplifies the test to see
- if it is within the needed bounds. */
- if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
- {
- mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- {
- mpz_sub_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- mpz_add (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (y));
- }
- }
+ 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_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
- if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
- {
- mpz_add_ui (SCM_I_BIG_MPZ (q),
- SCM_I_BIG_MPZ (q), 1);
- mpz_sub (SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (r),
- SCM_I_BIG_MPZ (y));
- }
+ mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
+ mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y));
}
- scm_remember_upto_here_2 (x, y);
- return scm_values (scm_list_2 (scm_i_normbig (q),
- scm_i_normbig (r)));
+
+ scm_remember_upto_here_2 (r2, y);
+ *qp = scm_i_normbig (q);
+ *rp = scm_i_normbig (r);
}
-/* Compute exact centered quotient and remainder the slow way.
- We use this only if both arguments are exact,
- and at least one of them is a fraction */
-static SCM
-scm_i_slow_exact_centered_divide (SCM x, SCM y)
+static void
+scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{
- SCM q, r;
+ SCM r1;
+ SCM xd = scm_denominator (x);
+ SCM yd = scm_denominator (y);
- if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
- s_scm_centered_divide);
- else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
- SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
- s_scm_centered_divide);
- else if (scm_is_true (scm_positive_p (y)))
- q = scm_floor (scm_sum (scm_divide (x, y),
- exactly_one_half));
- else if (scm_is_true (scm_negative_p (y)))
- q = scm_ceiling (scm_difference (scm_divide (x, y),
- exactly_one_half));
- else
- scm_num_overflow (s_scm_centered_divide);
- r = scm_difference (x, scm_product (q, y));
- return scm_values (scm_list_2 (q, r));
+ 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));
}
else if SCM_BIGP (n2)
{
intbig:
- if (n1 == 0)
+ if (nn1 == 0)
return SCM_INUM0;
{
SCM result_z = scm_i_mkbig ();
else if (SCM_BIGP (n))
{
char *str = mpz_get_str (NULL, base, SCM_I_BIG_MPZ (n));
+ size_t len = strlen (str);
+ void (*freefunc) (void *, size_t);
+ SCM ret;
+ mp_get_memory_functions (NULL, NULL, &freefunc);
scm_remember_upto_here_1 (n);
- return scm_take_locale_string (str);
+ ret = scm_from_latin1_stringn (str, len);
+ freefunc (str, len + 1);
+ return ret;
}
else if (SCM_FRACTIONP (n))
{
scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED)
{
char *str = mpz_get_str (NULL, 10, SCM_I_BIG_MPZ (exp));
+ size_t len = strlen (str);
+ void (*freefunc) (void *, size_t);
+ mp_get_memory_functions (NULL, NULL, &freefunc);
scm_remember_upto_here_1 (exp);
- scm_lfwrite (str, (size_t) strlen (str), port);
- free (str);
+ scm_lfwrite (str, len, port);
+ freefunc (str, len + 1);
return !0;
}
/*** END nums->strs ***/
* in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
* <uinteger R>, ...) that are used to build up numbers in the grammar. Some
* points should be noted about the implementation:
+ *
* * Each function keeps a local index variable 'idx' that points at the
* current position within the parsed string. The global index is only
* updated if the function could parse the corresponding syntactic unit
* successfully.
+ *
* * Similarly, the functions keep track of indicators of inexactness ('#',
- * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
- * global exactness information is only updated after each part has been
- * successfully parsed.
+ * '.' or exponents) using local variables ('hash_seen', 'x').
+ *
* * Sequences of digits are parsed into temporary variables holding fixnums.
* Only if these fixnums would overflow, the result variables are updated
* using the standard functions scm_add, scm_product, scm_divide etc. Then,
* digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
* and the result was computed as 12345 * 100000 + 67890. In other words,
* only every five digits two bignum operations were performed.
+ *
+ * Notes on the handling of exactness specifiers:
+ *
+ * When parsing non-real complex numbers, we apply exactness specifiers on
+ * per-component basis, as is done in PLT Scheme. For complex numbers
+ * written in rectangular form, exactness specifiers are applied to the
+ * real and imaginary parts before calling scm_make_rectangular. For
+ * complex numbers written in polar form, exactness specifiers are applied
+ * to the magnitude and angle before calling scm_make_polar.
+ *
+ * There are two kinds of exactness specifiers: forced and implicit. A
+ * forced exactness specifier is a "#e" or "#i" prefix at the beginning of
+ * the entire number, and applies to both components of a complex number.
+ * "#e" causes each component to be made exact, and "#i" causes each
+ * component to be made inexact. If no forced exactness specifier is
+ * present, then the exactness of each component is determined
+ * independently by the presence or absence of a decimal point or hash mark
+ * within that component. If a decimal point or hash mark is present, the
+ * component is made inexact, otherwise it is made exact.
+ *
+ * After the exactness specifiers have been applied to each component, they
+ * are passed to either scm_make_rectangular or scm_make_polar to produce
+ * the final result. Note that this will result in a real number if the
+ * imaginary part, magnitude, or angle is an exact 0.
+ *
+ * For example, (string->number "#i5.0+0i") does the equivalent of:
+ *
+ * (make-rectangular (exact->inexact 5) (exact->inexact 0))
*/
enum t_exactness {NO_EXACTNESS, INEXACT, EXACT};
return d;
}
+/* Parse the substring of MEM starting at *P_IDX for an unsigned integer
+ in base RADIX. Upon success, return the unsigned integer and update
+ *P_IDX and *P_EXACTNESS accordingly. Return #f on failure. */
static SCM
mem2uinteger (SCM mem, unsigned int *p_idx,
unsigned int radix, enum t_exactness *p_exactness)
if (sign == 1)
result = scm_product (result, e);
else
- result = scm_divide2real (result, e);
+ result = scm_divide (result, e);
/* We've seen an exponent, thus the value is implicitly inexact. */
x = INEXACT;
/* Cobble up the fractional part. We might want to set the
NaN's mantissa from it. */
idx += 4;
- mem2uinteger (mem, &idx, 10, &implicit_x);
+ if (!scm_is_eq (mem2uinteger (mem, &idx, 10, &implicit_x), SCM_INUM0))
+ {
+#if SCM_ENABLE_DEPRECATED == 1
+ scm_c_issue_deprecation_warning
+ ("Non-zero suffixes to `+nan.' are deprecated. Use `+nan.0'.");
+#else
+ return SCM_BOOL_F;
+#endif
+ }
+
*p_idx = idx;
return scm_nan ();
}
double
scm_c_truncate (double x)
{
-#if HAVE_TRUNC
return trunc (x);
-#else
- if (x < 0.0)
- return -floor (-x);
- return floor (x);
-#endif
}
/* scm_c_round is done using floor(x+0.5) to round to nearest and with
: result);
}
-SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0,
- (SCM x),
- "Round the number @var{x} towards zero.")
+SCM_PRIMITIVE_GENERIC (scm_truncate_number, "truncate", 1, 0, 0,
+ (SCM x),
+ "Round the number @var{x} towards zero.")
#define FUNC_NAME s_scm_truncate_number
{
- if (scm_is_false (scm_negative_p (x)))
- return scm_floor (x);
+ if (SCM_I_INUMP (x) || SCM_BIGP (x))
+ return x;
+ else if (SCM_REALP (x))
+ return scm_from_double (trunc (SCM_REAL_VALUE (x)));
+ else if (SCM_FRACTIONP (x))
+ return scm_truncate_quotient (SCM_FRACTION_NUMERATOR (x),
+ SCM_FRACTION_DENOMINATOR (x));
else
- return scm_ceiling (x);
+ SCM_WTA_DISPATCH_1 (g_scm_truncate_number, x, SCM_ARG1,
+ s_scm_truncate_number);
}
#undef FUNC_NAME
-SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
- (SCM x),
- "Round the number @var{x} towards the nearest integer. "
- "When it is exactly halfway between two integers, "
- "round towards the even one.")
+SCM_PRIMITIVE_GENERIC (scm_round_number, "round", 1, 0, 0,
+ (SCM x),
+ "Round the number @var{x} towards the nearest integer. "
+ "When it is exactly halfway between two integers, "
+ "round towards the even one.")
#define FUNC_NAME s_scm_round_number
{
if (SCM_I_INUMP (x) || SCM_BIGP (x))
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. */
- if (scm_is_true (scm_num_eq_p (plus_half, result))
- && scm_is_true (scm_odd_p (result)))
- return scm_difference (result, SCM_INUM1);
- else
- return result;
- }
+ SCM_WTA_DISPATCH_1 (g_scm_round_number, x, SCM_ARG1,
+ s_scm_round_number);
}
#undef FUNC_NAME
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);
}
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);
}
scm_c_make_rectangular (0, 2));
}
else
- SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan);
+ SCM_WTA_DISPATCH_1 (g_scm_atan, z, SCM_ARG1, s_scm_atan);
}
else if (scm_is_real (z))
{
SCM_DEFINE (scm_make_rectangular, "make-rectangular", 2, 0, 0,
(SCM real_part, SCM imaginary_part),
- "Return a complex number constructed of the given @var{real-part} "
- "and @var{imaginary-part} parts.")
+ "Return a complex number constructed of the given @var{real_part} "
+ "and @var{imaginary_part} parts.")
#define FUNC_NAME s_scm_make_rectangular
{
SCM_ASSERT_TYPE (scm_is_real (real_part), real_part,
}
+/* Returns log(x * 2^shift) */
+static SCM
+log_of_shifted_double (double x, long shift)
+{
+ double ans = log (fabs (x)) + shift * M_LN2;
+
+ if (x > 0.0 || double_is_non_negative_zero (x))
+ return scm_from_double (ans);
+ else
+ return scm_c_make_rectangular (ans, M_PI);
+}
+
+/* Returns log(n), for exact integer n of integer-length size */
+static SCM
+log_of_exact_integer_with_size (SCM n, long size)
+{
+ long shift = size - 2 * scm_dblprec[0];
+
+ if (shift > 0)
+ return log_of_shifted_double
+ (scm_to_double (scm_ash (n, scm_from_long(-shift))),
+ shift);
+ else
+ return log_of_shifted_double (scm_to_double (n), 0);
+}
+
+/* Returns log(n), for exact integer n */
+static SCM
+log_of_exact_integer (SCM n)
+{
+ return log_of_exact_integer_with_size
+ (n, scm_to_long (scm_integer_length (n)));
+}
+
+/* Returns log(n/d), for exact non-zero integers n and d */
+static SCM
+log_of_fraction (SCM n, SCM d)
+{
+ long n_size = scm_to_long (scm_integer_length (n));
+ long d_size = scm_to_long (scm_integer_length (d));
+
+ if (abs (n_size - d_size) > 1)
+ return (scm_difference (log_of_exact_integer_with_size (n, n_size),
+ log_of_exact_integer_with_size (d, d_size)));
+ else if (scm_is_false (scm_negative_p (n)))
+ return scm_from_double
+ (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
+ else
+ return scm_c_make_rectangular
+ (log1p (scm_to_double (scm_divide2real
+ (scm_difference (scm_abs (n), d),
+ d))),
+ M_PI);
+}
+
+
/* In the following functions we dispatch to the real-arg funcs like log()
when we know the arg is real, instead of just handing everything to
clog() for instance. This is in case clog() doesn't optimize for a
{
if (SCM_COMPLEXP (z))
{
-#if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
+#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG \
+ && defined (SCM_COMPLEX_VALUE)
return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z)));
#else
double re = SCM_COMPLEX_REAL (z);
atan2 (im, re));
#endif
}
- else if (SCM_NUMBERP (z))
+ else if (SCM_REALP (z))
+ return log_of_shifted_double (SCM_REAL_VALUE (z), 0);
+ else if (SCM_I_INUMP (z))
{
- /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
- although the value itself overflows. */
- double re = scm_to_double (z);
- double l = log (fabs (re));
- if (re >= 0.0)
- return scm_from_double (l);
- else
- return scm_c_make_rectangular (l, M_PI);
+#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
+ if (scm_is_eq (z, SCM_INUM0))
+ scm_num_overflow (s_scm_log);
+#endif
+ return log_of_shifted_double (SCM_I_INUM (z), 0);
}
+ else if (SCM_BIGP (z))
+ return log_of_exact_integer (z);
+ else if (SCM_FRACTIONP (z))
+ return log_of_fraction (SCM_FRACTION_NUMERATOR (z),
+ SCM_FRACTION_DENOMINATOR (z));
else
SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log);
}
M_LOG10E * atan2 (im, re));
#endif
}
- else if (SCM_NUMBERP (z))
+ else if (SCM_REALP (z) || SCM_I_INUMP (z))
{
- /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
- although the value itself overflows. */
- double re = scm_to_double (z);
- double l = log10 (fabs (re));
- if (re >= 0.0)
- return scm_from_double (l);
- else
- return scm_c_make_rectangular (l, M_LOG10E * M_PI);
+#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
+ if (scm_is_eq (z, SCM_INUM0))
+ scm_num_overflow (s_scm_log10);
+#endif
+ {
+ double re = scm_to_double (z);
+ double l = log10 (fabs (re));
+ if (re > 0.0 || double_is_non_negative_zero (re))
+ return scm_from_double (l);
+ else
+ return scm_c_make_rectangular (l, M_LOG10E * M_PI);
+ }
}
+ else if (SCM_BIGP (z))
+ return scm_product (flo_log10e, log_of_exact_integer (z));
+ else if (SCM_FRACTIONP (z))
+ return scm_product (flo_log10e,
+ log_of_fraction (SCM_FRACTION_NUMERATOR (z),
+ SCM_FRACTION_DENOMINATOR (z)));
else
SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10);
}
{
if (SCM_COMPLEXP (z))
{
-#if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
+#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CEXP \
+ && defined (SCM_COMPLEX_VALUE)
return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z)));
#else
return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z)),
#undef FUNC_NAME
+SCM_DEFINE (scm_i_exact_integer_sqrt, "exact-integer-sqrt", 1, 0, 0,
+ (SCM k),
+ "Return two exact non-negative integers @var{s} and @var{r}\n"
+ "such that @math{@var{k} = @var{s}^2 + @var{r}} and\n"
+ "@math{@var{s}^2 <= @var{k} < (@var{s} + 1)^2}.\n"
+ "An error is raised if @var{k} is not an exact non-negative integer.\n"
+ "\n"
+ "@lisp\n"
+ "(exact-integer-sqrt 10) @result{} 3 and 1\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_i_exact_integer_sqrt
+{
+ SCM s, r;
+
+ scm_exact_integer_sqrt (k, &s, &r);
+ return scm_values (scm_list_2 (s, r));
+}
+#undef FUNC_NAME
+
+void
+scm_exact_integer_sqrt (SCM k, SCM *sp, SCM *rp)
+{
+ if (SCM_LIKELY (SCM_I_INUMP (k)))
+ {
+ scm_t_inum kk = SCM_I_INUM (k);
+ scm_t_inum uu = kk;
+ scm_t_inum ss;
+
+ if (SCM_LIKELY (kk > 0))
+ {
+ do
+ {
+ ss = uu;
+ uu = (ss + kk/ss) / 2;
+ } while (uu < ss);
+ *sp = SCM_I_MAKINUM (ss);
+ *rp = SCM_I_MAKINUM (kk - ss*ss);
+ }
+ else if (SCM_LIKELY (kk == 0))
+ *sp = *rp = SCM_INUM0;
+ else
+ scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k,
+ "exact non-negative integer");
+ }
+ else if (SCM_LIKELY (SCM_BIGP (k)))
+ {
+ SCM s, r;
+
+ if (mpz_sgn (SCM_I_BIG_MPZ (k)) < 0)
+ scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k,
+ "exact non-negative integer");
+ s = scm_i_mkbig ();
+ r = scm_i_mkbig ();
+ mpz_sqrtrem (SCM_I_BIG_MPZ (s), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (k));
+ scm_remember_upto_here_1 (k);
+ *sp = scm_i_normbig (s);
+ *rp = scm_i_normbig (r);
+ }
+ else
+ scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k,
+ "exact non-negative integer");
+}
+
+
SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
(SCM z),
"Return the square root of @var{z}. Of the two possible roots\n"
- "(positive and negative), the one with the a positive real part\n"
+ "(positive and negative), the one with positive real part\n"
"is returned, or if that's zero then a positive imaginary part.\n"
"Thus,\n"
"\n"
{
int i;
+ if (scm_install_gmp_memory_functions)
+ mp_set_memory_functions (custom_gmp_malloc,
+ custom_gmp_realloc,
+ custom_gmp_free);
+
mpz_init_set_si (z_negative_one, -1);
/* It may be possible to tune the performance of some algorithms by using
scm_add_feature ("complex");
scm_add_feature ("inexact");
flo0 = scm_from_double (0.0);
+ flo_log10e = scm_from_double (M_LOG10E);
/* determine floating point precision */
for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
HCoop / bpt / guile.git / blobdiff