/* Copyright (C) 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003,
* 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012,
- * 2013 Free Software Foundation, Inc.
+ * 2013, 2014 Free Software Foundation, Inc.
*
* Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
* and Bellcore. See scm_divide.
#endif
#include <verify.h>
+#include <assert.h>
#include <math.h>
#include <string.h>
typedef scm_t_signed_bits scm_t_inum;
#define scm_from_inum(x) (scm_from_signed_integer (x))
-/* Tests to see if a C double is neither infinite nor a NaN.
- TODO: if it's available, use C99's isfinite(x) instead */
-#define DOUBLE_IS_FINITE(x) (!isinf(x) && !isnan(x))
-
-/* On some platforms, isinf(x) returns 0, 1 or -1, indicating the sign
- of the infinity, but other platforms return a boolean only. */
-#define DOUBLE_IS_POSITIVE_INFINITY(x) (isinf(x) && ((x) > 0))
-#define DOUBLE_IS_NEGATIVE_INFINITY(x) (isinf(x) && ((x) < 0))
-
+/* Test an inum to see if it can be converted to a double without loss
+ of precision. Note that this will sometimes return 0 even when 1
+ could have been returned, e.g. for large powers of 2. It is designed
+ to be a fast check to optimize common cases. */
+#define INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE(n) \
+ (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG \
+ || ((n) ^ ((n) >> (SCM_I_FIXNUM_BIT-1))) < (1L << DBL_MANT_DIG))
#if ! HAVE_DECL_MPZ_INITS
if (SCM_LIKELY (SCM_I_INUMP (d)))
{
- if (SCM_LIKELY (SCM_I_INUMP (n)
- && (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG
- || (SCM_I_INUM (n) < (1L << DBL_MANT_DIG)
- && SCM_I_INUM (d) < (1L << DBL_MANT_DIG)))))
+ if (SCM_LIKELY
+ (SCM_I_INUMP (n)
+ && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (n))
+ && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (d))))
/* If both N and D can be losslessly converted to doubles, then
we can rely on IEEE floating point to do proper rounding much
faster than we can. */
SCM_FRACTION_DENOMINATOR (z));
}
-static int
-double_is_non_negative_zero (double x)
+static SCM
+scm_i_from_double (double val)
{
- static double zero = 0.0;
+ SCM z;
+
+ z = SCM_PACK_POINTER (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real"));
+
+ SCM_SET_CELL_TYPE (z, scm_tc16_real);
+ SCM_REAL_VALUE (z) = val;
- return !memcmp (&x, &zero, sizeof(double));
+ return z;
}
SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0,
else if (SCM_REALP (n))
{
double val = SCM_REAL_VALUE (n);
- if (DOUBLE_IS_FINITE (val))
+ if (isfinite (val))
{
double rem = fabs (fmod (val, 2.0));
if (rem == 1.0)
else if (SCM_REALP (n))
{
double val = SCM_REAL_VALUE (n);
- if (DOUBLE_IS_FINITE (val))
+ if (isfinite (val))
{
double rem = fabs (fmod (val, 2.0));
if (rem == 1.0)
#define FUNC_NAME s_scm_finite_p
{
if (SCM_REALP (x))
- return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x)));
+ return scm_from_bool (isfinite (SCM_REAL_VALUE (x)));
else if (scm_is_real (x))
return SCM_BOOL_T;
else
guile_ieee_init ();
initialized = 1;
}
- return scm_from_double (guile_Inf);
+ return scm_i_from_double (guile_Inf);
}
#undef FUNC_NAME
guile_ieee_init ();
initialized = 1;
}
- return scm_from_double (guile_NaN);
+ return scm_i_from_double (guile_NaN);
}
#undef FUNC_NAME
double xx = SCM_REAL_VALUE (x);
/* If x is a NaN then xx<0 is false so we return x unchanged */
if (xx < 0.0)
- return scm_from_double (-xx);
+ return scm_i_from_double (-xx);
/* Handle signed zeroes properly */
else if (SCM_UNLIKELY (xx == 0.0))
return flo0;
if (SCM_UNLIKELY (y == 0))
scm_num_overflow (s_scm_floor_quotient); /* or return a NaN? */
else
- return scm_from_double (floor (x / y));
+ return scm_i_from_double (floor (x / y));
}
static SCM
if (SCM_UNLIKELY (y == 0))
scm_num_overflow (s_scm_floor_remainder); /* or return a NaN? */
else
- return scm_from_double (x - y * floor (x / y));
+ return scm_i_from_double (x - y * floor (x / y));
}
static SCM
{
double q = floor (x / y);
double r = x - q * y;
- *qp = scm_from_double (q);
- *rp = scm_from_double (r);
+ *qp = scm_i_from_double (q);
+ *rp = scm_i_from_double (r);
}
}
if (SCM_UNLIKELY (y == 0))
scm_num_overflow (s_scm_ceiling_quotient); /* or return a NaN? */
else
- return scm_from_double (ceil (x / y));
+ return scm_i_from_double (ceil (x / y));
}
static SCM
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));
+ return scm_i_from_double (x - y * ceil (x / y));
}
static SCM
{
double q = ceil (x / y);
double r = x - q * y;
- *qp = scm_from_double (q);
- *rp = scm_from_double (r);
+ *qp = scm_i_from_double (q);
+ *rp = scm_i_from_double (r);
}
}
if (SCM_UNLIKELY (y == 0))
scm_num_overflow (s_scm_truncate_quotient); /* or return a NaN? */
else
- return scm_from_double (trunc (x / y));
+ return scm_i_from_double (trunc (x / y));
}
static SCM
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));
+ return scm_i_from_double (x - y * trunc (x / y));
}
static SCM
{
double q = trunc (x / y);
double r = x - q * y;
- *qp = scm_from_double (q);
- *rp = scm_from_double (r);
+ *qp = scm_i_from_double (q);
+ *rp = scm_i_from_double (r);
}
}
scm_i_inexact_centered_quotient (double x, double y)
{
if (SCM_LIKELY (y > 0))
- return scm_from_double (floor (x/y + 0.5));
+ return scm_i_from_double (floor (x/y + 0.5));
else if (SCM_LIKELY (y < 0))
- return scm_from_double (ceil (x/y - 0.5));
+ return scm_i_from_double (ceil (x/y - 0.5));
else if (y == 0)
scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */
else
scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */
else
return scm_nan ();
- return scm_from_double (x - q * y);
+ return scm_i_from_double (x - q * y);
}
/* Assumes that both x and y are bigints, though
else
q = guile_NaN;
r = x - q * y;
- *qp = scm_from_double (q);
- *rp = scm_from_double (r);
+ *qp = scm_i_from_double (q);
+ *rp = scm_i_from_double (r);
}
/* Assumes that both x and y are bigints, though
if (SCM_UNLIKELY (y == 0))
scm_num_overflow (s_scm_round_quotient); /* or return a NaN? */
else
- return scm_from_double (scm_c_round (x / y));
+ return scm_i_from_double (scm_c_round (x / y));
}
/* Assumes that both x and y are bigints, though
else
{
double q = scm_c_round (x / y);
- return scm_from_double (x - q * y);
+ return scm_i_from_double (x - q * y);
}
}
{
double q = scm_c_round (x / y);
double r = x - q * y;
- *qp = scm_from_double (q);
- *rp = scm_from_double (r);
+ *qp = scm_i_from_double (q);
+ *rp = scm_i_from_double (r);
}
}
SCM_SWAP (x, y);
goto big_inum;
}
+ else if (SCM_REALP (y) && scm_is_integer (y))
+ goto handle_inexacts;
else
return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
}
scm_remember_upto_here_2 (x, y);
return scm_i_normbig (result);
}
+ else if (SCM_REALP (y) && scm_is_integer (y))
+ goto handle_inexacts;
+ else
+ return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
+ }
+ else if (SCM_REALP (x) && scm_is_integer (x))
+ {
+ if (SCM_I_INUMP (y) || SCM_BIGP (y)
+ || (SCM_REALP (y) && scm_is_integer (y)))
+ {
+ handle_inexacts:
+ return scm_exact_to_inexact (scm_gcd (scm_inexact_to_exact (x),
+ scm_inexact_to_exact (y)));
+ }
else
return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
}
SCM
scm_lcm (SCM n1, SCM n2)
{
- if (SCM_UNBNDP (n2))
- {
- if (SCM_UNBNDP (n1))
- return SCM_I_MAKINUM (1L);
- n2 = SCM_I_MAKINUM (1L);
- }
+ if (SCM_UNLIKELY (SCM_UNBNDP (n2)))
+ return SCM_UNBNDP (n1) ? SCM_INUM1 : scm_abs (n1);
- if (SCM_UNLIKELY (!(SCM_I_INUMP (n1) || SCM_BIGP (n1))))
- return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG1, s_lcm);
-
- if (SCM_UNLIKELY (!(SCM_I_INUMP (n2) || SCM_BIGP (n2))))
- return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG2, s_lcm);
-
- if (SCM_I_INUMP (n1))
+ if (SCM_LIKELY (SCM_I_INUMP (n1)))
{
- if (SCM_I_INUMP (n2))
+ if (SCM_LIKELY (SCM_I_INUMP (n2)))
{
SCM d = scm_gcd (n1, n2);
if (scm_is_eq (d, SCM_INUM0))
else
return scm_abs (scm_product (n1, scm_quotient (n2, d)));
}
- else
+ else if (SCM_LIKELY (SCM_BIGP (n2)))
{
/* inum n1, big n2 */
inumbig:
return result;
}
}
+ else if (SCM_REALP (n2) && scm_is_integer (n2))
+ goto handle_inexacts;
+ else
+ return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG2, s_lcm);
}
- else
+ else if (SCM_LIKELY (SCM_BIGP (n1)))
{
/* big n1 */
if (SCM_I_INUMP (n2))
SCM_SWAP (n1, n2);
goto inumbig;
}
- else
+ else if (SCM_LIKELY (SCM_BIGP (n2)))
{
SCM result = scm_i_mkbig ();
mpz_lcm(SCM_I_BIG_MPZ (result),
/* shouldn't need to normalize b/c lcm of 2 bigs should be big */
return result;
}
+ else if (SCM_REALP (n2) && scm_is_integer (n2))
+ goto handle_inexacts;
+ else
+ return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG2, s_lcm);
}
+ else if (SCM_REALP (n1) && scm_is_integer (n1))
+ {
+ if (SCM_I_INUMP (n2) || SCM_BIGP (n2)
+ || (SCM_REALP (n2) && scm_is_integer (n2)))
+ {
+ handle_inexacts:
+ return scm_exact_to_inexact (scm_lcm (scm_inexact_to_exact (n1),
+ scm_inexact_to_exact (n2)));
+ }
+ else
+ return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG2, s_lcm);
+ }
+ else
+ return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG1, s_lcm);
}
/* Emulating 2's complement bignums with sign magnitude arithmetic:
if (SCM_I_INUMP (j))
{
- /* bits above what's in an inum follow the sign bit */
- iindex = min (iindex, SCM_LONG_BIT - 1);
- return scm_from_bool ((1L << iindex) & SCM_I_INUM (j));
+ if (iindex < SCM_LONG_BIT - 1)
+ /* Arrange for the number to be converted to unsigned before
+ checking the bit, to ensure that we're testing the bit in a
+ two's complement representation (regardless of the native
+ representation. */
+ return scm_from_bool ((1UL << iindex) & SCM_I_INUM (j));
+ else
+ /* Portably check the sign. */
+ return scm_from_bool (SCM_I_INUM (j) < 0);
}
else if (SCM_BIGP (j))
{
{
scm_t_inum nn = SCM_I_INUM (n);
- /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will always
+ /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will almost[*] always
overflow a non-zero fixnum. For smaller shifts we check the
bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
- Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)". */
+ Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)".
+
+ [*] There's one exception:
+ (-1) << SCM_I_FIXNUM_BIT-1 == SCM_MOST_NEGATIVE_FIXNUM */
if (nn == 0)
return n;
else if (count < SCM_I_FIXNUM_BIT-1 &&
((scm_t_bits) (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - count)) + 1)
<= 1))
- return SCM_I_MAKINUM (nn << count);
+ return SCM_I_MAKINUM (nn < 0 ? -(-nn << count) : (nn << count));
else
{
SCM result = scm_i_inum2big (nn);
mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
count);
- return result;
+ return scm_i_normbig (result);
}
}
else if (SCM_BIGP (n))
return result;
}
else
- scm_syserror ("left_shift_exact_integer");
+ assert (0);
}
/* Efficiently compute floor (N / 2^COUNT),
return scm_i_normbig (result);
}
else
- scm_syserror ("floor_right_shift_exact_integer");
+ assert (0);
}
/* Efficiently compute round (N / 2^COUNT),
return scm_i_normbig (q);
}
else
- scm_syserror ("round_right_shift_exact_integer");
+ assert (0);
}
SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
}
else if (dbl == 0.0)
{
- if (!double_is_non_negative_zero (dbl))
+ if (copysign (1.0, dbl) < 0.0)
a[ch++] = '-';
strcpy (a + ch, "0.0");
return ch + 3;
#endif
/* Don't output a '+' for negative numbers or for Inf and
NaN. They will provide their own sign. */
- if (sgn >= 0 && DOUBLE_IS_FINITE (imag))
+ if (sgn >= 0 && isfinite (imag))
str[i++] = '+';
i += idbl2str (imag, &str[i], radix);
str[i++] = 'i';
static unsigned int
char_decimal_value (scm_t_uint32 c)
{
- /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
- that's certainly above any valid decimal, so we take advantage of
- that to elide some tests. */
- unsigned int d = (unsigned int) uc_decimal_value (c);
-
- /* If that failed, try extended hexadecimals, then. Only accept ascii
- hexadecimals. */
- if (d >= 10U)
+ if (c >= (scm_t_uint32) '0' && c <= (scm_t_uint32) '9')
+ return c - (scm_t_uint32) '0';
+ else
{
- c = uc_tolower (c);
- if (c >= (scm_t_uint32) 'a')
- d = c - (scm_t_uint32)'a' + 10U;
+ /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
+ that's certainly above any valid decimal, so we take advantage of
+ that to elide some tests. */
+ unsigned int d = (unsigned int) uc_decimal_value (c);
+
+ /* If that failed, try extended hexadecimals, then. Only accept ascii
+ hexadecimals. */
+ if (d >= 10U)
+ {
+ c = uc_tolower (c);
+ if (c >= (scm_t_uint32) 'a')
+ d = c - (scm_t_uint32)'a' + 10U;
+ }
+ return d;
}
- return d;
}
/* Parse the substring of MEM starting at *P_IDX for an unsigned integer
}
/* We should never get here */
- scm_syserror ("mem2ureal");
+ assert (0);
}
else if (SCM_REALP (x))
/* due to their limited precision, finite floating point numbers are
rational as well. (finite means neither infinity nor a NaN) */
- return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x)));
+ return scm_from_bool (isfinite (SCM_REAL_VALUE (x)));
else
return SCM_BOOL_F;
}
SCM_DEFINE (scm_integer_p, "integer?", 1, 0, 0,
(SCM x),
- "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
- "else.")
+ "Return @code{#t} if @var{x} is an integer number,\n"
+ "else return @code{#f}.")
#define FUNC_NAME s_scm_integer_p
{
if (SCM_I_INUMP (x) || SCM_BIGP (x))
}
#undef FUNC_NAME
+SCM_DEFINE (scm_exact_integer_p, "exact-integer?", 1, 0, 0,
+ (SCM x),
+ "Return @code{#t} if @var{x} is an exact integer number,\n"
+ "else return @code{#f}.")
+#define FUNC_NAME s_scm_exact_integer_p
+{
+ if (SCM_I_INUMP (x) || SCM_BIGP (x))
+ return SCM_BOOL_T;
+ else
+ return SCM_BOOL_F;
+}
+#undef FUNC_NAME
+
SCM scm_i_num_eq_p (SCM, SCM, SCM);
SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p, "=", 0, 2, 1,
to a double and compare.
But on a 64-bit system an inum is bigger than a double and
- casting it to a double (call that dxx) will round. dxx is at
- worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
- an integer and fits a long. So we cast yy to a long and
+ casting it to a double (call that dxx) will round.
+ Although dxx will not in general be equal to xx, dxx will
+ always be an integer and within a factor of 2 of xx, so if
+ dxx==yy, we know that yy is an integer and fits in
+ scm_t_signed_bits. So we cast yy to scm_t_signed_bits and
compare with plain xx.
An alternative (for any size system actually) would be to check
|| xx == (scm_t_signed_bits) yy));
}
else if (SCM_COMPLEXP (y))
- return scm_from_bool (((double) xx == SCM_COMPLEX_REAL (y))
- && (0.0 == SCM_COMPLEX_IMAG (y)));
+ {
+ /* see comments with inum/real above */
+ double ry = SCM_COMPLEX_REAL (y);
+ return scm_from_bool ((double) xx == ry
+ && 0.0 == SCM_COMPLEX_IMAG (y)
+ && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+ || xx == (scm_t_signed_bits) ry));
+ }
else if (SCM_FRACTIONP (y))
return SCM_BOOL_F;
else
else if (SCM_BIGP (y))
{
int cmp;
- if (isnan (SCM_REAL_VALUE (x)))
+ if (isnan (xx))
return SCM_BOOL_F;
- cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
+ cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), xx);
scm_remember_upto_here_1 (y);
return scm_from_bool (0 == cmp);
}
else if (SCM_REALP (y))
- return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
+ return scm_from_bool (xx == SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
- return scm_from_bool ((SCM_REAL_VALUE (x) == SCM_COMPLEX_REAL (y))
- && (0.0 == SCM_COMPLEX_IMAG (y)));
+ return scm_from_bool ((xx == SCM_COMPLEX_REAL (y))
+ && (0.0 == SCM_COMPLEX_IMAG (y)));
else if (SCM_FRACTIONP (y))
{
- double xx = SCM_REAL_VALUE (x);
- if (isnan (xx))
+ if (isnan (xx) || isinf (xx))
return SCM_BOOL_F;
- if (isinf (xx))
- return scm_from_bool (xx < 0.0);
x = scm_inexact_to_exact (x); /* with x as frac or int */
goto again;
}
else if (SCM_COMPLEXP (x))
{
if (SCM_I_INUMP (y))
- return scm_from_bool ((SCM_COMPLEX_REAL (x) == (double) SCM_I_INUM (y))
- && (SCM_COMPLEX_IMAG (x) == 0.0));
+ {
+ /* see comments with inum/real above */
+ double rx = SCM_COMPLEX_REAL (x);
+ scm_t_signed_bits yy = SCM_I_INUM (y);
+ return scm_from_bool (rx == (double) yy
+ && 0.0 == SCM_COMPLEX_IMAG (x)
+ && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+ || (scm_t_signed_bits) rx == yy));
+ }
else if (SCM_BIGP (y))
{
int cmp;
}
else if (SCM_REALP (y))
return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_REAL_VALUE (y))
- && (SCM_COMPLEX_IMAG (x) == 0.0));
+ && (SCM_COMPLEX_IMAG (x) == 0.0));
else if (SCM_COMPLEXP (y))
return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y))
- && (SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y)));
+ && (SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y)));
else if (SCM_FRACTIONP (y))
{
double xx;
if (SCM_COMPLEX_IMAG (x) != 0.0)
return SCM_BOOL_F;
xx = SCM_COMPLEX_REAL (x);
- if (isnan (xx))
+ if (isnan (xx) || isinf (xx))
return SCM_BOOL_F;
- if (isinf (xx))
- return scm_from_bool (xx < 0.0);
x = scm_inexact_to_exact (x); /* with x as frac or int */
goto again;
}
else if (SCM_REALP (y))
{
double yy = SCM_REAL_VALUE (y);
- if (isnan (yy))
+ if (isnan (yy) || isinf (yy))
return SCM_BOOL_F;
- if (isinf (yy))
- return scm_from_bool (0.0 < yy);
y = scm_inexact_to_exact (y); /* with y as frac or int */
goto again;
}
if (SCM_COMPLEX_IMAG (y) != 0.0)
return SCM_BOOL_F;
yy = SCM_COMPLEX_REAL (y);
- if (isnan (yy))
+ if (isnan (yy) || isinf(yy))
return SCM_BOOL_F;
- if (isinf (yy))
- return scm_from_bool (0.0 < yy);
y = scm_inexact_to_exact (y); /* with y as frac or int */
goto again;
}
return scm_from_bool (sgn > 0);
}
else if (SCM_REALP (y))
- return scm_from_bool ((double) xx < SCM_REAL_VALUE (y));
+ {
+ /* We can safely take the ceiling of y without changing the
+ result of x<y, given that x is an integer. */
+ double yy = ceil (SCM_REAL_VALUE (y));
+
+ /* In the following comparisons, it's important that the right
+ hand side always be a power of 2, so that it can be
+ losslessly converted to a double even on 64-bit
+ machines. */
+ if (yy >= (double) (SCM_MOST_POSITIVE_FIXNUM+1))
+ return SCM_BOOL_T;
+ else if (!(yy > (double) SCM_MOST_NEGATIVE_FIXNUM))
+ /* The condition above is carefully written to include the
+ case where yy==NaN. */
+ return SCM_BOOL_F;
+ else
+ /* yy is a finite integer that fits in an inum. */
+ return scm_from_bool (xx < (scm_t_inum) yy);
+ }
else if (SCM_FRACTIONP (y))
{
/* "x < a/b" becomes "x*b < a" */
else if (SCM_REALP (x))
{
if (SCM_I_INUMP (y))
- return scm_from_bool (SCM_REAL_VALUE (x) < (double) SCM_I_INUM (y));
+ {
+ /* We can safely take the floor of x without changing the
+ result of x<y, given that y is an integer. */
+ double xx = floor (SCM_REAL_VALUE (x));
+
+ /* In the following comparisons, it's important that the right
+ hand side always be a power of 2, so that it can be
+ losslessly converted to a double even on 64-bit
+ machines. */
+ if (xx < (double) SCM_MOST_NEGATIVE_FIXNUM)
+ return SCM_BOOL_T;
+ else if (!(xx < (double) (SCM_MOST_POSITIVE_FIXNUM+1)))
+ /* The condition above is carefully written to include the
+ case where xx==NaN. */
+ return SCM_BOOL_F;
+ else
+ /* xx is a finite integer that fits in an inum. */
+ return scm_from_bool ((scm_t_inum) xx < SCM_I_INUM (y));
+ }
else if (SCM_BIGP (y))
{
int cmp;
double yyd = SCM_REAL_VALUE (y);
if (xxd > yyd)
- return scm_from_double (xxd);
+ return scm_i_from_double (xxd);
/* If y is a NaN, then "==" is false and we return the NaN */
else if (SCM_LIKELY (!(xxd == yyd)))
return y;
big_real:
xx = scm_i_big2dbl (x);
yy = SCM_REAL_VALUE (y);
- return (xx > yy ? scm_from_double (xx) : y);
+ return (xx > yy ? scm_i_from_double (xx) : y);
}
else if (SCM_FRACTIONP (y))
{
double yyd = yy;
if (yyd > xxd)
- return scm_from_double (yyd);
+ return scm_i_from_double (yyd);
/* If x is a NaN, then "==" is false and we return the NaN */
else if (SCM_LIKELY (!(xxd == yyd)))
return x;
double xx = SCM_REAL_VALUE (x);
double yy = SCM_REAL_VALUE (y);
- /* For purposes of max: +inf.0 > nan > everything else, per R6RS */
+ /* For purposes of max: nan > +inf.0 > everything else,
+ per the R6RS errata */
if (xx > yy)
return x;
else if (SCM_LIKELY (xx < yy))
return y;
/* If neither (xx > yy) nor (xx < yy), then
either they're equal or one is a NaN */
- else if (SCM_UNLIKELY (isnan (xx)))
- return DOUBLE_IS_POSITIVE_INFINITY (yy) ? y : x;
- else if (SCM_UNLIKELY (isnan (yy)))
- return DOUBLE_IS_POSITIVE_INFINITY (xx) ? x : y;
+ else if (SCM_UNLIKELY (xx != yy))
+ return (xx != xx) ? x : y; /* Return the NaN */
/* xx == yy, but handle signed zeroes properly */
- else if (double_is_non_negative_zero (yy))
- return y;
- else
+ else if (copysign (1.0, yy) < 0.0)
return x;
+ else
+ return y;
}
else if (SCM_FRACTIONP (y))
{
double yy = scm_i_fraction2double (y);
double xx = SCM_REAL_VALUE (x);
- return (xx < yy) ? scm_from_double (yy) : x;
+ return (xx < yy) ? scm_i_from_double (yy) : x;
}
else
return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
{
double xx = scm_i_fraction2double (x);
/* if y==NaN then ">" is false, so we return the NaN y */
- return (xx > SCM_REAL_VALUE (y)) ? scm_from_double (xx) : y;
+ return (xx > SCM_REAL_VALUE (y)) ? scm_i_from_double (xx) : y;
}
else if (SCM_FRACTIONP (y))
{
{
double z = xx;
/* if y==NaN then "<" is false and we return NaN */
- return (z < SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
+ return (z < SCM_REAL_VALUE (y)) ? scm_i_from_double (z) : y;
}
else if (SCM_FRACTIONP (y))
{
big_real:
xx = scm_i_big2dbl (x);
yy = SCM_REAL_VALUE (y);
- return (xx < yy ? scm_from_double (xx) : y);
+ return (xx < yy ? scm_i_from_double (xx) : y);
}
else if (SCM_FRACTIONP (y))
{
{
double z = SCM_I_INUM (y);
/* if x==NaN then "<" is false and we return NaN */
- return (z < SCM_REAL_VALUE (x)) ? scm_from_double (z) : x;
+ return (z < SCM_REAL_VALUE (x)) ? scm_i_from_double (z) : x;
}
else if (SCM_BIGP (y))
{
double xx = SCM_REAL_VALUE (x);
double yy = SCM_REAL_VALUE (y);
- /* For purposes of min: -inf.0 < nan < everything else, per R6RS */
+ /* For purposes of min: nan < -inf.0 < everything else,
+ per the R6RS errata */
if (xx < yy)
return x;
else if (SCM_LIKELY (xx > yy))
return y;
/* If neither (xx < yy) nor (xx > yy), then
either they're equal or one is a NaN */
- else if (SCM_UNLIKELY (isnan (xx)))
- return DOUBLE_IS_NEGATIVE_INFINITY (yy) ? y : x;
- else if (SCM_UNLIKELY (isnan (yy)))
- return DOUBLE_IS_NEGATIVE_INFINITY (xx) ? x : y;
+ else if (SCM_UNLIKELY (xx != yy))
+ return (xx != xx) ? x : y; /* Return the NaN */
/* xx == yy, but handle signed zeroes properly */
- else if (double_is_non_negative_zero (xx))
- return y;
- else
+ else if (copysign (1.0, xx) < 0.0)
return x;
+ else
+ return y;
}
else if (SCM_FRACTIONP (y))
{
double yy = scm_i_fraction2double (y);
double xx = SCM_REAL_VALUE (x);
- return (yy < xx) ? scm_from_double (yy) : x;
+ return (yy < xx) ? scm_i_from_double (yy) : x;
}
else
return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
{
double xx = scm_i_fraction2double (x);
/* if y==NaN then "<" is false, so we return the NaN y */
- return (xx < SCM_REAL_VALUE (y)) ? scm_from_double (xx) : y;
+ return (xx < SCM_REAL_VALUE (y)) ? scm_i_from_double (xx) : y;
}
else if (SCM_FRACTIONP (y))
{
else if (SCM_REALP (y))
{
scm_t_inum xx = SCM_I_INUM (x);
- return scm_from_double (xx + SCM_REAL_VALUE (y));
+ return scm_i_from_double (xx + SCM_REAL_VALUE (y));
}
else if (SCM_COMPLEXP (y))
{
{
double result = mpz_get_d (SCM_I_BIG_MPZ (x)) + SCM_REAL_VALUE (y);
scm_remember_upto_here_1 (x);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_COMPLEXP (y))
{
else if (SCM_REALP (x))
{
if (SCM_I_INUMP (y))
- return scm_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
else if (SCM_BIGP (y))
{
double result = mpz_get_d (SCM_I_BIG_MPZ (y)) + SCM_REAL_VALUE (x);
scm_remember_upto_here_1 (y);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_REALP (y))
- return scm_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (SCM_REAL_VALUE (x) + SCM_COMPLEX_REAL (y),
SCM_COMPLEX_IMAG (y));
else if (SCM_FRACTIONP (y))
- return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
else
return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
}
scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
SCM_FRACTION_DENOMINATOR (x));
else if (SCM_REALP (y))
- return scm_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
+ return scm_i_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (SCM_COMPLEX_REAL (y) + scm_i_fraction2double (x),
SCM_COMPLEX_IMAG (y));
bignum, but negating that gives a fixnum. */
return scm_i_normbig (scm_i_clonebig (x, 0));
else if (SCM_REALP (x))
- return scm_from_double (-SCM_REAL_VALUE (x));
+ return scm_i_from_double (-SCM_REAL_VALUE (x));
else if (SCM_COMPLEXP (x))
return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
-SCM_COMPLEX_IMAG (x));
* (0.0 - 0.0) ==> 0.0, but (- 0.0) ==> -0.0.
*/
if (xx == 0)
- return scm_from_double (- SCM_REAL_VALUE (y));
+ return scm_i_from_double (- SCM_REAL_VALUE (y));
else
- return scm_from_double (xx - SCM_REAL_VALUE (y));
+ return scm_i_from_double (xx - SCM_REAL_VALUE (y));
}
else if (SCM_COMPLEXP (y))
{
{
double result = mpz_get_d (SCM_I_BIG_MPZ (x)) - SCM_REAL_VALUE (y);
scm_remember_upto_here_1 (x);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_COMPLEXP (y))
{
else if (SCM_REALP (x))
{
if (SCM_I_INUMP (y))
- return scm_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
else if (SCM_BIGP (y))
{
double result = SCM_REAL_VALUE (x) - mpz_get_d (SCM_I_BIG_MPZ (y));
scm_remember_upto_here_1 (x);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_REALP (y))
- return scm_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (SCM_REAL_VALUE (x) - SCM_COMPLEX_REAL (y),
-SCM_COMPLEX_IMAG (y));
else if (SCM_FRACTIONP (y))
- return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
else
return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
}
scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
SCM_FRACTION_DENOMINATOR (x));
else if (SCM_REALP (y))
- return scm_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
+ return scm_i_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (scm_i_fraction2double (x) - SCM_COMPLEX_REAL (y),
-SCM_COMPLEX_IMAG (y));
and we must do the multiplication in order to handle
infinities and NaNs properly. */
else if (SCM_REALP (y))
- return scm_from_double (0.0 * SCM_REAL_VALUE (y));
+ return scm_i_from_double (0.0 * SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y),
0.0 * SCM_COMPLEX_IMAG (y));
return result;
}
else if (SCM_REALP (y))
- return scm_from_double (xx * SCM_REAL_VALUE (y));
+ return scm_i_from_double (xx * SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
xx * SCM_COMPLEX_IMAG (y));
{
double result = mpz_get_d (SCM_I_BIG_MPZ (x)) * SCM_REAL_VALUE (y);
scm_remember_upto_here_1 (x);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_COMPLEXP (y))
{
{
double result = mpz_get_d (SCM_I_BIG_MPZ (y)) * SCM_REAL_VALUE (x);
scm_remember_upto_here_1 (y);
- return scm_from_double (result);
+ return scm_i_from_double (result);
}
else if (SCM_REALP (y))
- return scm_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
return scm_c_make_rectangular (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL (y),
SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
else if (SCM_FRACTIONP (y))
- return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
+ return scm_i_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
else
return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
}
return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
SCM_FRACTION_DENOMINATOR (x));
else if (SCM_REALP (y))
- return scm_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
+ return scm_i_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
else if (SCM_COMPLEXP (y))
{
double xx = scm_i_fraction2double (x);
scm_num_overflow (s_divide);
else
#endif
- return scm_from_double (1.0 / xx);
+ return scm_i_from_double (1.0 / xx);
}
else if (SCM_COMPLEXP (x))
{
#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
scm_num_overflow (s_divide);
#else
- return scm_from_double ((double) xx / (double) yy);
+ return scm_i_from_double ((double) xx / (double) yy);
#endif
}
else if (xx % yy != 0)
/* FIXME: Precision may be lost here due to:
(1) The cast from 'scm_t_inum' to 'double'
(2) Double rounding */
- return scm_from_double ((double) xx / yy);
+ return scm_i_from_double ((double) xx / yy);
}
else if (SCM_COMPLEXP (y))
{
#endif
/* FIXME: Precision may be lost here due to:
(1) scm_i_big2dbl (2) Double rounding */
- return scm_from_double (scm_i_big2dbl (x) / yy);
+ return scm_i_from_double (scm_i_big2dbl (x) / yy);
}
else if (SCM_COMPLEXP (y))
{
/* FIXME: Precision may be lost here due to:
(1) The cast from 'scm_t_inum' to 'double'
(2) Double rounding */
- return scm_from_double (rx / (double) yy);
+ return scm_i_from_double (rx / (double) yy);
}
else if (SCM_BIGP (y))
{
(2) Double rounding */
double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
scm_remember_upto_here_1 (y);
- return scm_from_double (rx / dby);
+ return scm_i_from_double (rx / dby);
}
else if (SCM_REALP (y))
{
scm_num_overflow (s_divide);
else
#endif
- return scm_from_double (rx / yy);
+ return scm_i_from_double (rx / yy);
}
else if (SCM_COMPLEXP (y))
{
goto complex_div;
}
else if (SCM_FRACTIONP (y))
- return scm_from_double (rx / scm_i_fraction2double (y));
+ return scm_i_from_double (rx / scm_i_fraction2double (y));
else
return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
}
/* FIXME: Precision may be lost here due to:
(1) The conversion from fraction to double
(2) Double rounding */
- return scm_from_double (scm_i_fraction2double (x) / yy);
+ return scm_i_from_double (scm_i_fraction2double (x) / yy);
}
else if (SCM_COMPLEXP (y))
{
if (SCM_I_INUMP (x) || SCM_BIGP (x))
return x;
else if (SCM_REALP (x))
- return scm_from_double (trunc (SCM_REAL_VALUE (x)));
+ return scm_i_from_double (trunc (SCM_REAL_VALUE (x)));
else if (SCM_FRACTIONP (x))
return scm_truncate_quotient (SCM_FRACTION_NUMERATOR (x),
SCM_FRACTION_DENOMINATOR (x));
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)));
+ return scm_i_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));
if (SCM_I_INUMP (x) || SCM_BIGP (x))
return x;
else if (SCM_REALP (x))
- return scm_from_double (floor (SCM_REAL_VALUE (x)));
+ return scm_i_from_double (floor (SCM_REAL_VALUE (x)));
else if (SCM_FRACTIONP (x))
return scm_floor_quotient (SCM_FRACTION_NUMERATOR (x),
SCM_FRACTION_DENOMINATOR (x));
if (SCM_I_INUMP (x) || SCM_BIGP (x))
return x;
else if (SCM_REALP (x))
- return scm_from_double (ceil (SCM_REAL_VALUE (x)));
+ return scm_i_from_double (ceil (SCM_REAL_VALUE (x)));
else if (SCM_FRACTIONP (x))
return scm_ceiling_quotient (SCM_FRACTION_NUMERATOR (x),
SCM_FRACTION_DENOMINATOR (x));
}
else if (scm_is_real (x) && scm_is_real (y) && scm_to_double (x) >= 0.0)
{
- return scm_from_double (pow (scm_to_double (x), scm_to_double (y)));
+ return scm_i_from_double (pow (scm_to_double (x), scm_to_double (y)));
}
else if (scm_is_complex (x) && scm_is_complex (y))
return scm_exp (scm_product (scm_log (x), y));
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* sin(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (sin (scm_to_double (z)));
+ return scm_i_from_double (sin (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y;
x = SCM_COMPLEX_REAL (z);
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return SCM_INUM1; /* cos(exact0) = exact1 */
else if (scm_is_real (z))
- return scm_from_double (cos (scm_to_double (z)));
+ return scm_i_from_double (cos (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y;
x = SCM_COMPLEX_REAL (z);
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* tan(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (tan (scm_to_double (z)));
+ return scm_i_from_double (tan (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y, w;
x = 2.0 * SCM_COMPLEX_REAL (z);
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* sinh(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (sinh (scm_to_double (z)));
+ return scm_i_from_double (sinh (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y;
x = SCM_COMPLEX_REAL (z);
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return SCM_INUM1; /* cosh(exact0) = exact1 */
else if (scm_is_real (z))
- return scm_from_double (cosh (scm_to_double (z)));
+ return scm_i_from_double (cosh (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y;
x = SCM_COMPLEX_REAL (z);
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* tanh(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (tanh (scm_to_double (z)));
+ return scm_i_from_double (tanh (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{ double x, y, w;
x = 2.0 * SCM_COMPLEX_REAL (z);
{
double w = scm_to_double (z);
if (w >= -1.0 && w <= 1.0)
- return scm_from_double (asin (w));
+ return scm_i_from_double (asin (w));
else
return scm_product (scm_c_make_rectangular (0, -1),
scm_sys_asinh (scm_c_make_rectangular (0, w)));
{
double w = scm_to_double (z);
if (w >= -1.0 && w <= 1.0)
- return scm_from_double (acos (w));
+ return scm_i_from_double (acos (w));
else
- return scm_sum (scm_from_double (acos (0.0)),
+ return scm_sum (scm_i_from_double (acos (0.0)),
scm_product (scm_c_make_rectangular (0, 1),
scm_sys_asinh (scm_c_make_rectangular (0, w))));
}
{ double x, y;
x = SCM_COMPLEX_REAL (z);
y = SCM_COMPLEX_IMAG (z);
- return scm_sum (scm_from_double (acos (0.0)),
+ return scm_sum (scm_i_from_double (acos (0.0)),
scm_product (scm_c_make_rectangular (0, 1),
scm_sys_asinh (scm_c_make_rectangular (-y, x))));
}
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* atan(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (atan (scm_to_double (z)));
+ return scm_i_from_double (atan (scm_to_double (z)));
else if (SCM_COMPLEXP (z))
{
double v, w;
else if (scm_is_real (z))
{
if (scm_is_real (y))
- return scm_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
+ return scm_i_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
else
return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan);
}
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* asinh(exact0) = exact0 */
else if (scm_is_real (z))
- return scm_from_double (asinh (scm_to_double (z)));
+ return scm_i_from_double (asinh (scm_to_double (z)));
else if (scm_is_number (z))
return scm_log (scm_sum (z,
scm_sqrt (scm_sum (scm_product (z, z),
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM1)))
return SCM_INUM0; /* acosh(exact1) = exact0 */
else if (scm_is_real (z) && scm_to_double (z) >= 1.0)
- return scm_from_double (acosh (scm_to_double (z)));
+ return scm_i_from_double (acosh (scm_to_double (z)));
else if (scm_is_number (z))
return scm_log (scm_sum (z,
scm_sqrt (scm_difference (scm_product (z, z),
if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
return z; /* atanh(exact0) = exact0 */
else if (scm_is_real (z) && scm_to_double (z) >= -1.0 && scm_to_double (z) <= 1.0)
- return scm_from_double (atanh (scm_to_double (z)));
+ return scm_i_from_double (atanh (scm_to_double (z)));
else if (scm_is_number (z))
return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1, z),
scm_difference (SCM_INUM1, z))),
#define FUNC_NAME s_scm_real_part
{
if (SCM_COMPLEXP (z))
- return scm_from_double (SCM_COMPLEX_REAL (z));
+ return scm_i_from_double (SCM_COMPLEX_REAL (z));
else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z))
return z;
else
#define FUNC_NAME s_scm_imag_part
{
if (SCM_COMPLEXP (z))
- return scm_from_double (SCM_COMPLEX_IMAG (z));
+ return scm_i_from_double (SCM_COMPLEX_IMAG (z));
else if (SCM_I_INUMP (z) || SCM_REALP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
return SCM_INUM0;
else
else if (SCM_FRACTIONP (z))
return SCM_FRACTION_NUMERATOR (z);
else if (SCM_REALP (z))
- return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
+ {
+ double zz = SCM_REAL_VALUE (z);
+ if (zz == floor (zz))
+ /* Handle -0.0 and infinities in accordance with R6RS
+ flnumerator, and optimize handling of integers. */
+ return z;
+ else
+ return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
+ }
else
return scm_wta_dispatch_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator);
}
else if (SCM_FRACTIONP (z))
return SCM_FRACTION_DENOMINATOR (z);
else if (SCM_REALP (z))
- return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
+ {
+ double zz = SCM_REAL_VALUE (z);
+ if (zz == floor (zz))
+ /* Handle infinities in accordance with R6RS fldenominator, and
+ optimize handling of integers. */
+ return scm_i_from_double (1.0);
+ else
+ return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
+ }
else
return scm_wta_dispatch_1 (g_scm_denominator, z, SCM_ARG1,
s_scm_denominator);
return z;
}
else if (SCM_REALP (z))
- return scm_from_double (fabs (SCM_REAL_VALUE (z)));
+ return scm_i_from_double (fabs (SCM_REAL_VALUE (z)));
else if (SCM_COMPLEXP (z))
- return scm_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
+ return scm_i_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
else if (SCM_FRACTIONP (z))
{
if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
#define FUNC_NAME s_scm_angle
{
/* atan(0,-1) is pi and it'd be possible to have that as a constant like
- flo0 to save allocating a new flonum with scm_from_double each time.
+ flo0 to save allocating a new flonum with scm_i_from_double each time.
But if atan2 follows the floating point rounding mode, then the value
is not a constant. Maybe it'd be close enough though. */
if (SCM_I_INUMP (z))
if (SCM_I_INUM (z) >= 0)
return flo0;
else
- return scm_from_double (atan2 (0.0, -1.0));
+ return scm_i_from_double (atan2 (0.0, -1.0));
}
else if (SCM_BIGP (z))
{
int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
scm_remember_upto_here_1 (z);
if (sgn < 0)
- return scm_from_double (atan2 (0.0, -1.0));
+ return scm_i_from_double (atan2 (0.0, -1.0));
else
return flo0;
}
else if (SCM_REALP (z))
{
double x = SCM_REAL_VALUE (z);
- if (x > 0.0 || double_is_non_negative_zero (x))
+ if (copysign (1.0, x) > 0.0)
return flo0;
else
- return scm_from_double (atan2 (0.0, -1.0));
+ return scm_i_from_double (atan2 (0.0, -1.0));
}
else if (SCM_COMPLEXP (z))
- return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
+ return scm_i_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
else if (SCM_FRACTIONP (z))
{
if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
return flo0;
- else return scm_from_double (atan2 (0.0, -1.0));
+ else return scm_i_from_double (atan2 (0.0, -1.0));
}
else
return scm_wta_dispatch_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
#define FUNC_NAME s_scm_exact_to_inexact
{
if (SCM_I_INUMP (z))
- return scm_from_double ((double) SCM_I_INUM (z));
+ return scm_i_from_double ((double) SCM_I_INUM (z));
else if (SCM_BIGP (z))
- return scm_from_double (scm_i_big2dbl (z));
+ return scm_i_from_double (scm_i_big2dbl (z));
else if (SCM_FRACTIONP (z))
- return scm_from_double (scm_i_fraction2double (z));
+ return scm_i_from_double (scm_i_fraction2double (z));
else if (SCM_INEXACTP (z))
return z;
else
return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1,
s_scm_inexact_to_exact);
- if (!SCM_LIKELY (DOUBLE_IS_FINITE (val)))
+ if (!SCM_LIKELY (isfinite (val)))
SCM_OUT_OF_RANGE (1, z);
else if (val == 0.0)
return SCM_INUM0;
{
SCM_ASSERT_TYPE (scm_is_real (x), x, SCM_ARG1, FUNC_NAME, "real");
SCM_ASSERT_TYPE (scm_is_real (eps), eps, SCM_ARG2, FUNC_NAME, "real");
- eps = scm_abs (eps);
- if (scm_is_false (scm_positive_p (eps)))
- {
- /* eps is either zero or a NaN */
- if (scm_is_true (scm_nan_p (eps)))
- return scm_nan ();
- else if (SCM_INEXACTP (eps))
- return scm_exact_to_inexact (x);
- else
- return x;
- }
- else if (scm_is_false (scm_finite_p (eps)))
- {
- if (scm_is_true (scm_finite_p (x)))
- return flo0;
- else
- return scm_nan ();
- }
- else if (scm_is_false (scm_finite_p (x))) /* checks for both inf and nan */
- return x;
- else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x, eps)),
- scm_ceiling (scm_difference (x, eps)))))
+
+ if (SCM_UNLIKELY (!scm_is_exact (eps) || !scm_is_exact (x)))
{
- /* There's an integer within range; we want the one closest to zero */
- if (scm_is_false (scm_less_p (eps, scm_abs (x))))
- {
- /* zero is within range */
- if (SCM_INEXACTP (x) || SCM_INEXACTP (eps))
- return flo0;
- else
- return SCM_INUM0;
- }
- else if (scm_is_true (scm_positive_p (x)))
- return scm_ceiling (scm_difference (x, eps));
+ if (SCM_UNLIKELY (scm_is_false (scm_finite_p (eps))))
+ {
+ if (scm_is_false (scm_nan_p (eps)) && scm_is_true (scm_finite_p (x)))
+ return flo0;
+ else
+ return scm_nan ();
+ }
+ else if (SCM_UNLIKELY (scm_is_false (scm_finite_p (x))))
+ return x;
else
- return scm_floor (scm_sum (x, eps));
- }
- else
- {
- /* Use continued fractions to find closest ratio. All
- arithmetic is done with exact numbers.
+ return scm_exact_to_inexact
+ (scm_rationalize (scm_inexact_to_exact (x),
+ scm_inexact_to_exact (eps)));
+ }
+ else
+ {
+ /* X and EPS are exact rationals.
+
+ The code that follows is equivalent to the following Scheme code:
+
+ (define (exact-rationalize x eps)
+ (let ((n1 (if (negative? x) -1 1))
+ (x (abs x))
+ (eps (abs eps)))
+ (let ((lo (- x eps))
+ (hi (+ x eps)))
+ (if (<= lo 0)
+ 0
+ (let loop ((nlo (numerator lo)) (dlo (denominator lo))
+ (nhi (numerator hi)) (dhi (denominator hi))
+ (n1 n1) (d1 0) (n2 0) (d2 1))
+ (let-values (((qlo rlo) (floor/ nlo dlo))
+ ((qhi rhi) (floor/ nhi dhi)))
+ (let ((n0 (+ n2 (* n1 qlo)))
+ (d0 (+ d2 (* d1 qlo))))
+ (cond ((zero? rlo) (/ n0 d0))
+ ((< qlo qhi) (/ (+ n0 n1) (+ d0 d1)))
+ (else (loop dhi rhi dlo rlo n0 d0 n1 d1))))))))))
*/
- SCM ex = scm_inexact_to_exact (x);
- SCM int_part = scm_floor (ex);
- SCM tt = SCM_INUM1;
- SCM a1 = SCM_INUM0, a2 = SCM_INUM1, a = SCM_INUM0;
- SCM b1 = SCM_INUM1, b2 = SCM_INUM0, b = SCM_INUM0;
- SCM rx;
- int i = 0;
+ int n1_init = 1;
+ SCM lo, hi;
- ex = scm_difference (ex, int_part); /* x = x-int_part */
- rx = scm_divide (ex, SCM_UNDEFINED); /* rx = 1/x */
+ eps = scm_abs (eps);
+ if (scm_is_true (scm_negative_p (x)))
+ {
+ n1_init = -1;
+ x = scm_difference (x, SCM_UNDEFINED);
+ }
- /* We stop after a million iterations just to be absolutely sure
- that we don't go into an infinite loop. The process normally
- converges after less than a dozen iterations.
- */
+ /* X and EPS are non-negative exact rationals. */
- while (++i < 1000000)
- {
- a = scm_sum (scm_product (a1, tt), a2); /* a = a1*tt + a2 */
- b = scm_sum (scm_product (b1, tt), b2); /* b = b1*tt + b2 */
- if (scm_is_false (scm_zero_p (b)) && /* b != 0 */
- scm_is_false
- (scm_gr_p (scm_abs (scm_difference (ex, scm_divide (a, b))),
- eps))) /* abs(x-a/b) <= eps */
- {
- SCM res = scm_sum (int_part, scm_divide (a, b));
- if (SCM_INEXACTP (x) || SCM_INEXACTP (eps))
- return scm_exact_to_inexact (res);
- else
- return res;
- }
- rx = scm_divide (scm_difference (rx, tt), /* rx = 1/(rx - tt) */
- SCM_UNDEFINED);
- tt = scm_floor (rx); /* tt = floor (rx) */
- a2 = a1;
- b2 = b1;
- a1 = a;
- b1 = b;
- }
- scm_num_overflow (s_scm_rationalize);
+ lo = scm_difference (x, eps);
+ hi = scm_sum (x, eps);
+
+ if (scm_is_false (scm_positive_p (lo)))
+ /* If zero is included in the interval, return it.
+ It is the simplest rational of all. */
+ return SCM_INUM0;
+ else
+ {
+ SCM result;
+ mpz_t n0, d0, n1, d1, n2, d2;
+ mpz_t nlo, dlo, nhi, dhi;
+ mpz_t qlo, rlo, qhi, rhi;
+
+ /* LO and HI are positive exact rationals. */
+
+ /* Our approach here follows the method described by Alan
+ Bawden in a message entitled "(rationalize x y)" on the
+ rrrs-authors mailing list, dated 16 Feb 1988 14:08:28 EST:
+
+ http://groups.csail.mit.edu/mac/ftpdir/scheme-mail/HTML/rrrs-1988/msg00063.html
+
+ In brief, we compute the continued fractions of the two
+ endpoints of the interval (LO and HI). The continued
+ fraction of the result consists of the common prefix of the
+ continued fractions of LO and HI, plus one final term. The
+ final term of the result is the smallest integer contained
+ in the interval between the remainders of LO and HI after
+ the common prefix has been removed.
+
+ The following code lazily computes the continued fraction
+ representations of LO and HI, and simultaneously converts
+ the continued fraction of the result into a rational
+ number. We use MPZ functions directly to avoid type
+ dispatch and GC allocation during the loop. */
+
+ mpz_inits (n0, d0, n1, d1, n2, d2,
+ nlo, dlo, nhi, dhi,
+ qlo, rlo, qhi, rhi,
+ NULL);
+
+ /* The variables N1, D1, N2 and D2 are used to compute the
+ resulting rational from its continued fraction. At each
+ step, N2/D2 and N1/D1 are the last two convergents. They
+ are normally initialized to 0/1 and 1/0, respectively.
+ However, if we negated X then we must negate the result as
+ well, and we do that by initializing N1/D1 to -1/0. */
+ mpz_set_si (n1, n1_init);
+ mpz_set_ui (d1, 0);
+ mpz_set_ui (n2, 0);
+ mpz_set_ui (d2, 1);
+
+ /* The variables NLO, DLO, NHI, and DHI are used to lazily
+ compute the continued fraction representations of LO and HI
+ using Euclid's algorithm. Initially, NLO/DLO == LO and
+ NHI/DHI == HI. */
+ scm_to_mpz (scm_numerator (lo), nlo);
+ scm_to_mpz (scm_denominator (lo), dlo);
+ scm_to_mpz (scm_numerator (hi), nhi);
+ scm_to_mpz (scm_denominator (hi), dhi);
+
+ /* As long as we're using exact arithmetic, the following loop
+ is guaranteed to terminate. */
+ for (;;)
+ {
+ /* Compute the next terms (QLO and QHI) of the continued
+ fractions of LO and HI. */
+ mpz_fdiv_qr (qlo, rlo, nlo, dlo); /* QLO <-- floor (NLO/DLO), RLO <-- NLO - QLO * DLO */
+ mpz_fdiv_qr (qhi, rhi, nhi, dhi); /* QHI <-- floor (NHI/DHI), RHI <-- NHI - QHI * DHI */
+
+ /* The next term of the result will be either QLO or
+ QLO+1. Here we compute the next convergent of the
+ result based on the assumption that QLO is the next
+ term. If that turns out to be wrong, we'll adjust
+ these later by adding N1 to N0 and D1 to D0. */
+ mpz_set (n0, n2); mpz_addmul (n0, n1, qlo); /* N0 <-- N2 + (QLO * N1) */
+ mpz_set (d0, d2); mpz_addmul (d0, d1, qlo); /* D0 <-- D2 + (QLO * D1) */
+
+ /* We stop iterating when an integer is contained in the
+ interval between the remainders NLO/DLO and NHI/DHI.
+ There are two cases to consider: either NLO/DLO == QLO
+ is an integer (indicated by RLO == 0), or QLO < QHI. */
+ if (mpz_sgn (rlo) == 0 || mpz_cmp (qlo, qhi) != 0)
+ break;
+
+ /* Efficiently shuffle variables around for the next
+ iteration. First we shift the recent convergents. */
+ mpz_swap (n2, n1); mpz_swap (n1, n0); /* N2 <-- N1 <-- N0 */
+ mpz_swap (d2, d1); mpz_swap (d1, d0); /* D2 <-- D1 <-- D0 */
+
+ /* The following shuffling is a bit confusing, so some
+ explanation is in order. Conceptually, we're doing a
+ couple of things here. After substracting the floor of
+ NLO/DLO, the remainder is RLO/DLO. The rest of the
+ continued fraction will represent the remainder's
+ reciprocal DLO/RLO. Similarly for the HI endpoint.
+ So in the next iteration, the new endpoints will be
+ DLO/RLO and DHI/RHI. However, when we take the
+ reciprocals of these endpoints, their order is
+ switched. So in summary, we want NLO/DLO <-- DHI/RHI
+ and NHI/DHI <-- DLO/RLO. */
+ mpz_swap (nlo, dhi); mpz_swap (dhi, rlo); /* NLO <-- DHI <-- RLO */
+ mpz_swap (nhi, dlo); mpz_swap (dlo, rhi); /* NHI <-- DLO <-- RHI */
+ }
+
+ /* There is now an integer in the interval [NLO/DLO NHI/DHI].
+ The last term of the result will be the smallest integer in
+ that interval, which is ceiling(NLO/DLO). We have already
+ computed floor(NLO/DLO) in QLO, so now we adjust QLO to be
+ equal to the ceiling. */
+ if (mpz_sgn (rlo) != 0)
+ {
+ /* If RLO is non-zero, then NLO/DLO is not an integer and
+ the next term will be QLO+1. QLO was used in the
+ computation of N0 and D0 above. Here we adjust N0 and
+ D0 to be based on QLO+1 instead of QLO. */
+ mpz_add (n0, n0, n1); /* N0 <-- N0 + N1 */
+ mpz_add (d0, d0, d1); /* D0 <-- D0 + D1 */
+ }
+
+ /* The simplest rational in the interval is N0/D0 */
+ result = scm_i_make_ratio_already_reduced (scm_from_mpz (n0),
+ scm_from_mpz (d0));
+ mpz_clears (n0, d0, n1, d1, n2, d2,
+ nlo, dlo, nhi, dhi,
+ qlo, rlo, qhi, rhi,
+ NULL);
+ return result;
+ }
}
}
#undef FUNC_NAME
return scm_is_true (scm_integer_p (val));
}
+int
+scm_is_exact_integer (SCM val)
+{
+ return scm_is_true (scm_exact_integer_p (val));
+}
+
int
scm_is_signed_integer (SCM val, scm_t_intmax min, scm_t_intmax max)
{
SCM
scm_from_double (double val)
{
- SCM z;
-
- z = SCM_PACK_POINTER (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real"));
-
- SCM_SET_CELL_TYPE (z, scm_tc16_real);
- SCM_REAL_VALUE (z) = val;
-
- return z;
+ return scm_i_from_double (val);
}
int
{
double ans = log (fabs (x)) + shift * M_LN2;
- if (x > 0.0 || double_is_non_negative_zero (x))
- return scm_from_double (ans);
+ if (copysign (1.0, x) > 0.0)
+ return scm_i_from_double (ans);
else
return scm_c_make_rectangular (ans, M_PI);
}
return (scm_difference (log_of_exact_integer (n),
log_of_exact_integer (d)));
else if (scm_is_false (scm_negative_p (n)))
- return scm_from_double
+ return scm_i_from_double
(log1p (scm_i_divide2double (scm_difference (n, d), d)));
else
return scm_c_make_rectangular
{
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);
+ if (copysign (1.0, re) > 0.0)
+ return scm_i_from_double (l);
else
return scm_c_make_rectangular (l, M_LOG10E * M_PI);
}
{
/* When z is a negative bignum the conversion to double overflows,
giving -infinity, but that's ok, the exp is still 0.0. */
- return scm_from_double (exp (scm_to_double (z)));
+ return scm_i_from_double (exp (scm_to_double (z)));
}
else
return scm_wta_dispatch_1 (g_scm_exp, z, 1, s_scm_exp);
if (root == floor (root))
return SCM_I_MAKINUM ((scm_t_inum) root);
else
- return scm_from_double (root);
+ return scm_i_from_double (root);
}
else
{
return scm_c_make_rectangular
(0.0, ldexp (sqrt (-signif), expon / 2));
else
- return scm_from_double (ldexp (sqrt (signif), expon / 2));
+ return scm_i_from_double (ldexp (sqrt (signif), expon / 2));
}
}
else if (SCM_FRACTIONP (z))
if (xx < 0)
return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift));
else
- return scm_from_double (ldexp (sqrt (xx), shift));
+ return scm_i_from_double (ldexp (sqrt (xx), shift));
}
}
if (xx < 0)
return scm_c_make_rectangular (0.0, sqrt (-xx));
else
- return scm_from_double (sqrt (xx));
+ return scm_i_from_double (sqrt (xx));
}
}
else
scm_add_feature ("complex");
scm_add_feature ("inexact");
- flo0 = scm_from_double (0.0);
- flo_log10e = scm_from_double (M_LOG10E);
+ flo0 = scm_i_from_double (0.0);
+ flo_log10e = scm_i_from_double (M_LOG10E);
exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
HCoop / bpt / guile.git / blobdiff