X-Git-Url: http://git.hcoop.net/bpt/guile.git/blobdiff_plain/7f34acd8a48198c7fec2daf8d2f4161eaa9963ec..28d5d2537c0321643c3b511a2195cd491204e7f2:/libguile/numbers.c diff --git a/libguile/numbers.c b/libguile/numbers.c index fa55b4f9e..9857e182d 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -1,4 +1,6 @@ -/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc. +/* 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. * * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories * and Bellcore. See scm_divide. @@ -56,6 +58,8 @@ #include #endif +#include + #include "libguile/_scm.h" #include "libguile/feature.h" #include "libguile/ports.h" @@ -96,6 +100,38 @@ typedef scm_t_signed_bits scm_t_inum; #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 + +/* GMP < 5.0.0 lacks `mpz_inits' and `mpz_clears'. Provide them. */ + +#define VARARG_MPZ_ITERATOR(func) \ + static void \ + func ## s (mpz_t x, ...) \ + { \ + va_list ap; \ + \ + va_start (ap, x); \ + while (x != NULL) \ + { \ + func (x); \ + x = va_arg (ap, mpz_ptr); \ + } \ + va_end (ap); \ + } + +VARARG_MPZ_ITERATOR (mpz_init) +VARARG_MPZ_ITERATOR (mpz_clear) + +#endif + /* @@ -186,7 +222,7 @@ finalize_bignum (void *ptr, void *data) { SCM bignum; - bignum = PTR2SCM (ptr); + bignum = SCM_PACK_POINTER (ptr); mpz_clear (SCM_I_BIG_MPZ (bignum)); } @@ -410,9 +446,6 @@ scm_i_mpz2num (mpz_t b) } } -/* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */ -static SCM scm_divide2real (SCM x, SCM y); - /* Make the ratio NUMERATOR/DENOMINATOR, where: 1. NUMERATOR and DENOMINATOR are exact integers 2. NUMERATOR and DENOMINATOR are reduced to lowest terms: gcd(n,d) == 1 */ @@ -466,11 +499,159 @@ scm_i_make_ratio (SCM numerator, SCM denominator) } #undef FUNC_NAME +static mpz_t scm_i_divide2double_lo2b; + +/* Return the double that is closest to the exact rational N/D, with + ties rounded toward even mantissas. N and D must be exact + integers. */ +static double +scm_i_divide2double (SCM n, SCM d) +{ + int neg; + mpz_t nn, dd, lo, hi, x; + ssize_t e; + + if (SCM_LIKELY (SCM_I_INUMP (d))) + { + 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. */ + return ((double) SCM_I_INUM (n)) / ((double) SCM_I_INUM (d)); + + if (SCM_UNLIKELY (scm_is_eq (d, SCM_INUM0))) + { + if (scm_is_true (scm_positive_p (n))) + return 1.0 / 0.0; + else if (scm_is_true (scm_negative_p (n))) + return -1.0 / 0.0; + else + return 0.0 / 0.0; + } + + mpz_init_set_si (dd, SCM_I_INUM (d)); + } + else + mpz_init_set (dd, SCM_I_BIG_MPZ (d)); + + if (SCM_I_INUMP (n)) + mpz_init_set_si (nn, SCM_I_INUM (n)); + else + mpz_init_set (nn, SCM_I_BIG_MPZ (n)); + + neg = (mpz_sgn (nn) < 0) ^ (mpz_sgn (dd) < 0); + mpz_abs (nn, nn); + mpz_abs (dd, dd); + + /* Now we need to find the value of e such that: + + For e <= 0: + b^{p-1} - 1/2b <= b^-e n / d < b^p - 1/2 [1A] + (2 b^p - 1) <= 2 b b^-e n / d < (2 b^p - 1) b [2A] + (2 b^p - 1) d <= 2 b b^-e n < (2 b^p - 1) d b [3A] + + For e >= 0: + b^{p-1} - 1/2b <= n / b^e d < b^p - 1/2 [1B] + (2 b^p - 1) <= 2 b n / b^e d < (2 b^p - 1) b [2B] + (2 b^p - 1) d b^e <= 2 b n < (2 b^p - 1) d b b^e [3B] + + where: p = DBL_MANT_DIG + b = FLT_RADIX (here assumed to be 2) + + After rounding, the mantissa must be an integer between b^{p-1} and + (b^p - 1), except for subnormal numbers. In the inequations [1A] + and [1B], the middle expression represents the mantissa *before* + rounding, and therefore is bounded by the range of values that will + round to a floating-point number with the exponent e. The upper + bound is (b^p - 1 + 1/2) = (b^p - 1/2), and is exclusive because + ties will round up to the next power of b. The lower bound is + (b^{p-1} - 1/2b), and is inclusive because ties will round toward + this power of b. Here we subtract 1/2b instead of 1/2 because it + is in the range of the next smaller exponent, where the + representable numbers are closer together by a factor of b. + + Inequations [2A] and [2B] are derived from [1A] and [1B] by + multiplying by 2b, and in [3A] and [3B] we multiply by the + denominator of the middle value to obtain integer expressions. + + In the code below, we refer to the three expressions in [3A] or + [3B] as lo, x, and hi. If the number is normalizable, we will + achieve the goal: lo <= x < hi */ + + /* Make an initial guess for e */ + e = mpz_sizeinbase (nn, 2) - mpz_sizeinbase (dd, 2) - (DBL_MANT_DIG-1); + if (e < DBL_MIN_EXP - DBL_MANT_DIG) + e = DBL_MIN_EXP - DBL_MANT_DIG; + + /* Compute the initial values of lo, x, and hi + based on the initial guess of e */ + mpz_inits (lo, hi, x, NULL); + mpz_mul_2exp (x, nn, 2 + ((e < 0) ? -e : 0)); + mpz_mul (lo, dd, scm_i_divide2double_lo2b); + if (e > 0) + mpz_mul_2exp (lo, lo, e); + mpz_mul_2exp (hi, lo, 1); + + /* Adjust e as needed to satisfy the inequality lo <= x < hi, + (but without making e less then the minimum exponent) */ + while (mpz_cmp (x, lo) < 0 && e > DBL_MIN_EXP - DBL_MANT_DIG) + { + mpz_mul_2exp (x, x, 1); + e--; + } + while (mpz_cmp (x, hi) >= 0) + { + /* If we ever used lo's value again, + we would need to double lo here. */ + mpz_mul_2exp (hi, hi, 1); + e++; + } + + /* Now compute the rounded mantissa: + n / b^e d (if e >= 0) + n b^-e / d (if e <= 0) */ + { + int cmp; + double result; + + if (e < 0) + mpz_mul_2exp (nn, nn, -e); + else + mpz_mul_2exp (dd, dd, e); + + /* mpz does not directly support rounded right + shifts, so we have to do it the hard way. + For efficiency, we reuse lo and hi. + hi == quotient, lo == remainder */ + mpz_fdiv_qr (hi, lo, nn, dd); + + /* The fractional part of the unrounded mantissa would be + remainder/dividend, i.e. lo/dd. So we have a tie if + lo/dd = 1/2. Multiplying both sides by 2*dd yields the + integer expression 2*lo = dd. Here we do that comparison + to decide whether to round up or down. */ + mpz_mul_2exp (lo, lo, 1); + cmp = mpz_cmp (lo, dd); + if (cmp > 0 || (cmp == 0 && mpz_odd_p (hi))) + mpz_add_ui (hi, hi, 1); + + result = ldexp (mpz_get_d (hi), e); + if (neg) + result = -result; + + mpz_clears (nn, dd, lo, hi, x, NULL); + return result; + } +} + double scm_i_fraction2double (SCM z) { - return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z), - SCM_FRACTION_DENOMINATOR (z))); + return scm_i_divide2double (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z)); } static int @@ -492,7 +673,7 @@ SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0, else if (SCM_NUMBERP (x)) return SCM_BOOL_T; else - SCM_WTA_DISPATCH_1 (g_scm_exact_p, x, 1, s_scm_exact_p); + return scm_wta_dispatch_1 (g_scm_exact_p, x, 1, s_scm_exact_p); } #undef FUNC_NAME @@ -513,7 +694,7 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_p, "inexact?", 1, 0, 0, else if (SCM_NUMBERP (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p); + return scm_wta_dispatch_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p); } #undef FUNC_NAME @@ -552,7 +733,7 @@ SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0, return SCM_BOOL_F; } } - SCM_WTA_DISPATCH_1 (g_scm_odd_p, n, 1, s_scm_odd_p); + return scm_wta_dispatch_1 (g_scm_odd_p, n, 1, s_scm_odd_p); } #undef FUNC_NAME @@ -586,7 +767,7 @@ SCM_PRIMITIVE_GENERIC (scm_even_p, "even?", 1, 0, 0, return SCM_BOOL_T; } } - SCM_WTA_DISPATCH_1 (g_scm_even_p, n, 1, s_scm_even_p); + return scm_wta_dispatch_1 (g_scm_even_p, n, 1, s_scm_even_p); } #undef FUNC_NAME @@ -601,7 +782,7 @@ SCM_PRIMITIVE_GENERIC (scm_finite_p, "finite?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_T; else - SCM_WTA_DISPATCH_1 (g_scm_finite_p, x, 1, s_scm_finite_p); + return scm_wta_dispatch_1 (g_scm_finite_p, x, 1, s_scm_finite_p); } #undef FUNC_NAME @@ -616,7 +797,7 @@ SCM_PRIMITIVE_GENERIC (scm_inf_p, "inf?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_inf_p, x, 1, s_scm_inf_p); + return scm_wta_dispatch_1 (g_scm_inf_p, x, 1, s_scm_inf_p); } #undef FUNC_NAME @@ -631,7 +812,7 @@ SCM_PRIMITIVE_GENERIC (scm_nan_p, "nan?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_nan_p, x, 1, s_scm_nan_p); + return scm_wta_dispatch_1 (g_scm_nan_p, x, 1, s_scm_nan_p); } #undef FUNC_NAME @@ -759,7 +940,7 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0, SCM_FRACTION_DENOMINATOR (x)); } else - SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs); + return scm_wta_dispatch_1 (g_scm_abs, x, 1, s_scm_abs); } #undef FUNC_NAME @@ -774,10 +955,10 @@ SCM_PRIMITIVE_GENERIC (scm_quotient, "quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient); + return scm_wta_dispatch_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient); } #undef FUNC_NAME @@ -795,10 +976,10 @@ SCM_PRIMITIVE_GENERIC (scm_remainder, "remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder); + return scm_wta_dispatch_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder); } #undef FUNC_NAME @@ -817,10 +998,10 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo); } else - SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG1, s_scm_modulo); + return scm_wta_dispatch_2 (g_scm_modulo, x, y, SCM_ARG1, s_scm_modulo); } #undef FUNC_NAME @@ -919,10 +1100,9 @@ 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 vals = scm_wta_dispatch_2 (gf, a1, a2, pos, subr); + + scm_i_extract_values_2 (vals, rp1, rp2); } SCM_DEFINE (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, @@ -1054,8 +1234,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, - s_scm_floor_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_BIGP (x)) { @@ -1095,8 +1275,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, - s_scm_floor_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_REALP (x)) { @@ -1105,8 +1285,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, return scm_i_inexact_floor_quotient (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, - s_scm_floor_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_FRACTIONP (x)) { @@ -1116,12 +1296,12 @@ SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG2, - s_scm_floor_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_floor_quotient, x, y, SCM_ARG1, - s_scm_floor_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG1, + s_scm_floor_quotient); } #undef FUNC_NAME @@ -1214,8 +1394,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, - s_scm_floor_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_BIGP (x)) { @@ -1250,8 +1430,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, - s_scm_floor_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_REALP (x)) { @@ -1260,8 +1440,8 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, return scm_i_inexact_floor_remainder (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, - s_scm_floor_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_FRACTIONP (x)) { @@ -1271,12 +1451,12 @@ SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_exact_rational_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG2, - s_scm_floor_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_floor_remainder, x, y, SCM_ARG1, - s_scm_floor_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG1, + s_scm_floor_remainder); } #undef FUNC_NAME @@ -1587,8 +1767,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_ceiling_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, - s_scm_ceiling_quotient); + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); } else if (SCM_BIGP (x)) { @@ -1628,8 +1808,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_ceiling_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, - s_scm_ceiling_quotient); + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); } else if (SCM_REALP (x)) { @@ -1638,8 +1818,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, return scm_i_inexact_ceiling_quotient (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, - s_scm_ceiling_quotient); + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); } else if (SCM_FRACTIONP (x)) { @@ -1649,12 +1829,12 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_exact_rational_ceiling_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, - s_scm_ceiling_quotient); + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, - s_scm_ceiling_quotient); + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, + s_scm_ceiling_quotient); } #undef FUNC_NAME @@ -1757,8 +1937,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else if (SCM_BIGP (x)) { @@ -1793,8 +1973,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else if (SCM_REALP (x)) { @@ -1803,8 +1983,8 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else if (SCM_FRACTIONP (x)) { @@ -1814,12 +1994,12 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, + s_scm_ceiling_remainder); } #undef FUNC_NAME @@ -2119,8 +2299,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); } else if (SCM_BIGP (x)) { @@ -2160,8 +2340,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); } else if (SCM_REALP (x)) { @@ -2170,8 +2350,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); } else if (SCM_FRACTIONP (x)) { @@ -2181,12 +2361,12 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, + s_scm_truncate_quotient); } #undef FUNC_NAME @@ -2254,8 +2434,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); } else if (SCM_BIGP (x)) { @@ -2288,8 +2468,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); } else if (SCM_REALP (x)) { @@ -2298,8 +2478,8 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); } else if (SCM_FRACTIONP (x)) { @@ -2309,12 +2489,12 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, + s_scm_truncate_remainder); } #undef FUNC_NAME @@ -2601,8 +2781,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); } else if (SCM_BIGP (x)) { @@ -2650,8 +2830,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); } else if (SCM_REALP (x)) { @@ -2660,8 +2840,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); } else if (SCM_FRACTIONP (x)) { @@ -2671,12 +2851,12 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG1, + s_scm_centered_quotient); } #undef FUNC_NAME @@ -2815,8 +2995,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); } else if (SCM_BIGP (x)) { @@ -2856,8 +3036,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); } else if (SCM_REALP (x)) { @@ -2866,8 +3046,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); } else if (SCM_FRACTIONP (x)) { @@ -2877,12 +3057,12 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG1, + s_scm_centered_remainder); } #undef FUNC_NAME @@ -3297,8 +3477,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_round_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, - s_scm_round_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_BIGP (x)) { @@ -3348,8 +3528,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_round_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, - s_scm_round_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_REALP (x)) { @@ -3358,8 +3538,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, return scm_i_inexact_round_quotient (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, - s_scm_round_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_FRACTIONP (x)) { @@ -3369,12 +3549,12 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) return scm_i_exact_rational_round_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG2, - s_scm_round_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_round_quotient, x, y, SCM_ARG1, - s_scm_round_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG1, + s_scm_round_quotient); } #undef FUNC_NAME @@ -3501,8 +3681,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_BIGP (x)) { @@ -3549,8 +3729,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_REALP (x)) { @@ -3559,8 +3739,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, 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); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_FRACTIONP (x)) { @@ -3570,12 +3750,12 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, 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); + return 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); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG1, + s_scm_round_remainder); } #undef FUNC_NAME @@ -3965,7 +4145,7 @@ scm_gcd (SCM x, SCM y) goto big_inum; } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd); } else if (SCM_BIGP (x)) { @@ -3995,10 +4175,10 @@ scm_gcd (SCM x, SCM y) return scm_i_normbig (result); } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd); } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG1, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG1, s_gcd); } SCM_PRIMITIVE_GENERIC (scm_i_lcm, "lcm", 0, 2, 1, @@ -4029,10 +4209,11 @@ scm_lcm (SCM n1, SCM n2) n2 = SCM_I_MAKINUM (1L); } - SCM_GASSERT2 (SCM_I_INUMP (n1) || SCM_BIGP (n1), - g_lcm, n1, n2, SCM_ARG1, s_lcm); - SCM_GASSERT2 (SCM_I_INUMP (n2) || SCM_BIGP (n2), - g_lcm, n1, n2, SCM_ARGn, s_lcm); + 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)) { @@ -5115,229 +5296,230 @@ SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0, #undef FUNC_NAME /*** NUMBERS -> STRINGS ***/ -#define SCM_MAX_DBL_PREC 60 #define SCM_MAX_DBL_RADIX 36 -/* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */ -static int scm_dblprec[SCM_MAX_DBL_RADIX - 1]; -static double fx_per_radix[SCM_MAX_DBL_RADIX - 1][SCM_MAX_DBL_PREC]; - -static -void init_dblprec(int *prec, int radix) { - /* determine floating point precision by adding successively - smaller increments to 1.0 until it is considered == 1.0 */ - double f = ((double)1.0)/radix; - double fsum = 1.0 + f; - - *prec = 0; - while (fsum != 1.0) - { - if (++(*prec) > SCM_MAX_DBL_PREC) - fsum = 1.0; - else - { - f /= radix; - fsum = f + 1.0; - } - } - (*prec) -= 1; -} - -static -void init_fx_radix(double *fx_list, int radix) -{ - /* initialize a per-radix list of tolerances. When added - to a number < 1.0, we can determine if we should raund - up and quit converting a number to a string. */ - int i; - fx_list[0] = 0.0; - fx_list[1] = 0.5; - for( i = 2 ; i < SCM_MAX_DBL_PREC; ++i ) - fx_list[i] = (fx_list[i-1] / radix); -} - /* use this array as a way to generate a single digit */ static const char number_chars[] = "0123456789abcdefghijklmnopqrstuvwxyz"; +static mpz_t dbl_minimum_normal_mantissa; + static size_t -idbl2str (double f, char *a, int radix) +idbl2str (double dbl, char *a, int radix) { - int efmt, dpt, d, i, wp; - double *fx; -#ifdef DBL_MIN_10_EXP - double f_cpy; - int exp_cpy; -#endif /* DBL_MIN_10_EXP */ - size_t ch = 0; - int exp = 0; - - if(radix < 2 || - radix > SCM_MAX_DBL_RADIX) - { - /* revert to existing behavior */ - radix = 10; - } + int ch = 0; - wp = scm_dblprec[radix-2]; - fx = fx_per_radix[radix-2]; + if (radix < 2 || radix > SCM_MAX_DBL_RADIX) + /* revert to existing behavior */ + radix = 10; - if (f == 0.0) + if (isinf (dbl)) { -#ifdef HAVE_COPYSIGN - double sgn = copysign (1.0, f); - - if (sgn < 0.0) - a[ch++] = '-'; -#endif - goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */ + strcpy (a, (dbl > 0.0) ? "+inf.0" : "-inf.0"); + return 6; } - - if (isinf (f)) + else if (dbl > 0.0) + ; + else if (dbl < 0.0) { - if (f < 0) - strcpy (a, "-inf.0"); - else - strcpy (a, "+inf.0"); - return ch+6; + dbl = -dbl; + a[ch++] = '-'; } - else if (isnan (f)) + else if (dbl == 0.0) { - strcpy (a, "+nan.0"); - return ch+6; + if (!double_is_non_negative_zero (dbl)) + a[ch++] = '-'; + strcpy (a + ch, "0.0"); + return ch + 3; } - - if (f < 0.0) + else if (isnan (dbl)) { - f = -f; - a[ch++] = '-'; + strcpy (a, "+nan.0"); + return 6; } -#ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from - make-uniform-vector, from causing infinite loops. */ - /* just do the checking...if it passes, we do the conversion for our - radix again below */ - f_cpy = f; - exp_cpy = exp; + /* Algorithm taken from "Printing Floating-Point Numbers Quickly and + Accurately" by Robert G. Burger and R. Kent Dybvig */ + { + int e, k; + mpz_t f, r, s, mplus, mminus, hi, digit; + int f_is_even, f_is_odd; + int expon; + int show_exp = 0; + + mpz_inits (f, r, s, mplus, mminus, hi, digit, NULL); + mpz_set_d (f, ldexp (frexp (dbl, &e), DBL_MANT_DIG)); + if (e < DBL_MIN_EXP) + { + mpz_tdiv_q_2exp (f, f, DBL_MIN_EXP - e); + e = DBL_MIN_EXP; + } + e -= DBL_MANT_DIG; - while (f_cpy < 1.0) - { - f_cpy *= 10.0; - if (exp_cpy-- < DBL_MIN_10_EXP) - { - a[ch++] = '#'; - a[ch++] = '.'; - a[ch++] = '#'; - return ch; - } - } - while (f_cpy > 10.0) - { - f_cpy *= 0.10; - if (exp_cpy++ > DBL_MAX_10_EXP) - { - a[ch++] = '#'; - a[ch++] = '.'; - a[ch++] = '#'; - return ch; - } - } -#endif + f_is_even = !mpz_odd_p (f); + f_is_odd = !f_is_even; - while (f < 1.0) - { - f *= radix; - exp--; - } - while (f > radix) - { - f /= radix; - exp++; - } + /* Initialize r, s, mplus, and mminus according + to Table 1 from the paper. */ + if (e < 0) + { + mpz_set_ui (mminus, 1); + if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0 + || e == DBL_MIN_EXP - DBL_MANT_DIG) + { + mpz_set_ui (mplus, 1); + mpz_mul_2exp (r, f, 1); + mpz_mul_2exp (s, mminus, 1 - e); + } + else + { + mpz_set_ui (mplus, 2); + mpz_mul_2exp (r, f, 2); + mpz_mul_2exp (s, mminus, 2 - e); + } + } + else + { + mpz_set_ui (mminus, 1); + mpz_mul_2exp (mminus, mminus, e); + if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0) + { + mpz_set (mplus, mminus); + mpz_mul_2exp (r, f, 1 + e); + mpz_set_ui (s, 2); + } + else + { + mpz_mul_2exp (mplus, mminus, 1); + mpz_mul_2exp (r, f, 2 + e); + mpz_set_ui (s, 4); + } + } - if (f + fx[wp] >= radix) - { - f = 1.0; - exp++; - } - zero: -#ifdef ENGNOT - /* adding 9999 makes this equivalent to abs(x) % 3 */ - dpt = (exp + 9999) % 3; - exp -= dpt++; - efmt = 1; -#else - efmt = (exp < -3) || (exp > wp + 2); - if (!efmt) + /* Find the smallest k such that: + (r + mplus) / s < radix^k (if f is even) + (r + mplus) / s <= radix^k (if f is odd) */ { - if (exp < 0) - { - a[ch++] = '0'; - a[ch++] = '.'; - dpt = exp; - while (++dpt) - a[ch++] = '0'; - } - else - dpt = exp + 1; + /* IMPROVE-ME: Make an initial guess to speed this up */ + mpz_add (hi, r, mplus); + k = 0; + while (mpz_cmp (hi, s) >= f_is_odd) + { + mpz_mul_ui (s, s, radix); + k++; + } + if (k == 0) + { + mpz_mul_ui (hi, hi, radix); + while (mpz_cmp (hi, s) < f_is_odd) + { + mpz_mul_ui (r, r, radix); + mpz_mul_ui (mplus, mplus, radix); + mpz_mul_ui (mminus, mminus, radix); + mpz_mul_ui (hi, hi, radix); + k--; + } + } } - else - dpt = 1; -#endif - do - { - d = f; - f -= d; - a[ch++] = number_chars[d]; - if (f < fx[wp]) - break; - if (f + fx[wp] >= 1.0) - { - a[ch - 1] = number_chars[d+1]; - break; - } - f *= radix; - if (!(--dpt)) - a[ch++] = '.'; - } - while (wp--); + expon = k - 1; + if (k <= 0) + { + if (k <= -3) + { + /* Use scientific notation */ + show_exp = 1; + k = 1; + } + else + { + int i; - if (dpt > 0) - { -#ifndef ENGNOT - if ((dpt > 4) && (exp > 6)) - { - d = (a[0] == '-' ? 2 : 1); - for (i = ch++; i > d; i--) - a[i] = a[i - 1]; - a[d] = '.'; - efmt = 1; - } - else -#endif - { - while (--dpt) - a[ch++] = '0'; - a[ch++] = '.'; - } - } - if (a[ch - 1] == '.') - a[ch++] = '0'; /* trailing zero */ - if (efmt && exp) - { - a[ch++] = 'e'; - if (exp < 0) - { - exp = -exp; - a[ch++] = '-'; - } - for (i = radix; i <= exp; i *= radix); - for (i /= radix; i; i /= radix) - { - a[ch++] = number_chars[exp / i]; - exp %= i; - } - } + /* Print leading zeroes */ + a[ch++] = '0'; + a[ch++] = '.'; + for (i = 0; i > k; i--) + a[ch++] = '0'; + } + } + + for (;;) + { + int end_1_p, end_2_p; + int d; + + mpz_mul_ui (mplus, mplus, radix); + mpz_mul_ui (mminus, mminus, radix); + mpz_mul_ui (r, r, radix); + mpz_fdiv_qr (digit, r, r, s); + d = mpz_get_ui (digit); + + mpz_add (hi, r, mplus); + end_1_p = (mpz_cmp (r, mminus) < f_is_even); + end_2_p = (mpz_cmp (s, hi) < f_is_even); + if (end_1_p || end_2_p) + { + mpz_mul_2exp (r, r, 1); + if (!end_2_p) + ; + else if (!end_1_p) + d++; + else if (mpz_cmp (r, s) >= !(d & 1)) + d++; + a[ch++] = number_chars[d]; + if (--k == 0) + a[ch++] = '.'; + break; + } + else + { + a[ch++] = number_chars[d]; + if (--k == 0) + a[ch++] = '.'; + } + } + + if (k > 0) + { + if (expon >= 7 && k >= 4 && expon >= k) + { + /* Here we would have to print more than three zeroes + followed by a decimal point and another zero. It + makes more sense to use scientific notation. */ + + /* Adjust k to what it would have been if we had chosen + scientific notation from the beginning. */ + k -= expon; + + /* k will now be <= 0, with magnitude equal to the number of + digits that we printed which should now be put after the + decimal point. */ + + /* Insert a decimal point */ + memmove (a + ch + k + 1, a + ch + k, -k); + a[ch + k] = '.'; + ch++; + + show_exp = 1; + } + else + { + for (; k > 0; k--) + a[ch++] = '0'; + a[ch++] = '.'; + } + } + + if (k == 0) + a[ch++] = '0'; + + if (show_exp) + { + a[ch++] = 'e'; + ch += scm_iint2str (expon, radix, a + ch); + } + + mpz_clears (f, r, s, mplus, mminus, hi, digit, NULL); + } return ch; } @@ -5476,7 +5658,7 @@ int scm_print_real (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED) { char num_buf[FLOBUFLEN]; - scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port); + scm_lfwrite_unlocked (num_buf, iflo2str (sexp, num_buf, 10), port); return !0; } @@ -5484,7 +5666,7 @@ void scm_i_print_double (double val, SCM port) { char num_buf[FLOBUFLEN]; - scm_lfwrite (num_buf, idbl2str (val, num_buf, 10), port); + scm_lfwrite_unlocked (num_buf, idbl2str (val, num_buf, 10), port); } int @@ -5492,7 +5674,7 @@ scm_print_complex (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED) { char num_buf[FLOBUFLEN]; - scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port); + scm_lfwrite_unlocked (num_buf, iflo2str (sexp, num_buf, 10), port); return !0; } @@ -5500,7 +5682,7 @@ void scm_i_print_complex (double real, double imag, SCM port) { char num_buf[FLOBUFLEN]; - scm_lfwrite (num_buf, icmplx2str (real, imag, num_buf, 10), port); + scm_lfwrite_unlocked (num_buf, icmplx2str (real, imag, num_buf, 10), port); } int @@ -5521,7 +5703,7 @@ scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED) void (*freefunc) (void *, size_t); mp_get_memory_functions (NULL, NULL, &freefunc); scm_remember_upto_here_1 (exp); - scm_lfwrite (str, len, port); + scm_lfwrite_unlocked (str, len, port); freefunc (str, len + 1); return !0; } @@ -5821,7 +6003,7 @@ mem2decimal_from_point (SCM result, SCM mem, break; } - if (exponent > SCM_MAXEXP) + if (exponent > ((sign == 1) ? SCM_MAXEXP : SCM_MAXEXP + DBL_DIG + 1)) { size_t exp_len = idx - start; SCM exp_string = scm_i_substring_copy (mem, start, start + exp_len); @@ -6360,9 +6542,11 @@ scm_num_eq_p (SCM x, SCM y) 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 @@ -6377,12 +6561,19 @@ scm_num_eq_p (SCM x, SCM y) || 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 - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_BIGP (x)) { @@ -6417,7 +6608,8 @@ scm_num_eq_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_REALP (x)) { @@ -6433,35 +6625,40 @@ scm_num_eq_p (SCM x, SCM y) 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 - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } 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; @@ -6475,25 +6672,24 @@ scm_num_eq_p (SCM x, SCM y) } 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 - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_FRACTIONP (x)) { @@ -6504,10 +6700,8 @@ scm_num_eq_p (SCM x, SCM y) 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; } @@ -6517,20 +6711,20 @@ scm_num_eq_p (SCM x, SCM y) 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; } else if (SCM_FRACTIONP (y)) return scm_i_fraction_equalp (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARG1, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARG1, + s_scm_i_num_eq_p); } @@ -6579,7 +6773,25 @@ scm_less_p (SCM x, SCM y) 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= (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" */ @@ -6589,7 +6801,8 @@ scm_less_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } else if (SCM_BIGP (x)) { @@ -6617,12 +6830,31 @@ scm_less_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) goto int_frac; else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } 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 xinexact", 1, 0, 0, else if (SCM_INEXACTP (z)) return z; else - SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact, z, 1, s_scm_exact_to_inexact); + return scm_wta_dispatch_1 (g_scm_exact_to_inexact, z, 1, + s_scm_exact_to_inexact); } #undef FUNC_NAME @@ -9105,26 +9324,40 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, else if (SCM_COMPLEXP (z) && SCM_COMPLEX_IMAG (z) == 0.0) val = SCM_COMPLEX_REAL (z); else - SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact, z, 1, s_scm_inexact_to_exact); + return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1, + s_scm_inexact_to_exact); if (!SCM_LIKELY (DOUBLE_IS_FINITE (val))) SCM_OUT_OF_RANGE (1, z); + else if (val == 0.0) + return SCM_INUM0; else { - mpq_t frac; - SCM q; - - mpq_init (frac); - mpq_set_d (frac, val); - q = scm_i_make_ratio_already_reduced - (scm_i_mpz2num (mpq_numref (frac)), - scm_i_mpz2num (mpq_denref (frac))); - - /* When scm_i_make_ratio throws, we leak the memory allocated - for frac... - */ - mpq_clear (frac); - return q; + int expon; + SCM numerator; + + numerator = scm_i_dbl2big (ldexp (frexp (val, &expon), + DBL_MANT_DIG)); + expon -= DBL_MANT_DIG; + if (expon < 0) + { + int shift = mpz_scan1 (SCM_I_BIG_MPZ (numerator), 0); + + if (shift > -expon) + shift = -expon; + mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (numerator), + SCM_I_BIG_MPZ (numerator), + shift); + expon += shift; + } + numerator = scm_i_normbig (numerator); + if (expon < 0) + return scm_i_make_ratio_already_reduced + (numerator, left_shift_exact_integer (SCM_INUM1, -expon)); + else if (expon > 0) + return left_shift_exact_integer (numerator, expon); + else + return numerator; } } } @@ -9486,7 +9719,7 @@ scm_from_double (double val) { SCM z; - z = PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real")); + 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; @@ -9494,46 +9727,6 @@ scm_from_double (double val) return z; } -#if SCM_ENABLE_DEPRECATED == 1 - -float -scm_num2float (SCM num, unsigned long pos, const char *s_caller) -{ - scm_c_issue_deprecation_warning - ("`scm_num2float' is deprecated. Use scm_to_double instead."); - - if (SCM_BIGP (num)) - { - float res = mpz_get_d (SCM_I_BIG_MPZ (num)); - if (!isinf (res)) - return res; - else - scm_out_of_range (NULL, num); - } - else - return scm_to_double (num); -} - -double -scm_num2double (SCM num, unsigned long pos, const char *s_caller) -{ - scm_c_issue_deprecation_warning - ("`scm_num2double' is deprecated. Use scm_to_double instead."); - - if (SCM_BIGP (num)) - { - double res = mpz_get_d (SCM_I_BIG_MPZ (num)); - if (!isinf (res)) - return res; - else - scm_out_of_range (NULL, num); - } - else - return scm_to_double (num); -} - -#endif - int scm_is_complex (SCM val) { @@ -9628,12 +9821,11 @@ log_of_fraction (SCM n, SCM d) log_of_exact_integer (d))); else if (scm_is_false (scm_negative_p (n))) return scm_from_double - (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d)))); + (log1p (scm_i_divide2double (scm_difference (n, d), d))); else return scm_c_make_rectangular - (log1p (scm_to_double (scm_divide2real - (scm_difference (scm_abs (n), d), - d))), + (log1p (scm_i_divide2double (scm_difference (scm_abs (n), d), + d)), M_PI); } @@ -9677,7 +9869,7 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0, 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); + return scm_wta_dispatch_1 (g_scm_log, z, 1, s_scm_log); } #undef FUNC_NAME @@ -9724,7 +9916,7 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0, log_of_fraction (SCM_FRACTION_NUMERATOR (z), SCM_FRACTION_DENOMINATOR (z))); else - SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10); + return scm_wta_dispatch_1 (g_scm_log10, z, 1, s_scm_log10); } #undef FUNC_NAME @@ -9752,7 +9944,7 @@ SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0, return scm_from_double (exp (scm_to_double (z))); } else - SCM_WTA_DISPATCH_1 (g_scm_exp, z, 1, s_scm_exp); + return scm_wta_dispatch_1 (g_scm_exp, z, 1, s_scm_exp); } #undef FUNC_NAME @@ -9781,25 +9973,17 @@ 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; + mpz_t kk, ss, rr; - 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 + if (SCM_I_INUM (k) < 0) scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k, "exact non-negative integer"); + mpz_init_set_ui (kk, SCM_I_INUM (k)); + mpz_inits (ss, rr, NULL); + mpz_sqrtrem (ss, rr, kk); + *sp = SCM_I_MAKINUM (mpz_get_ui (ss)); + *rp = SCM_I_MAKINUM (mpz_get_ui (rr)); + mpz_clears (kk, ss, rr, NULL); } else if (SCM_LIKELY (SCM_BIGP (k))) { @@ -9820,6 +10004,56 @@ scm_exact_integer_sqrt (SCM k, SCM *sp, SCM *rp) "exact non-negative integer"); } +/* Return true iff K is a perfect square. + K must be an exact integer. */ +static int +exact_integer_is_perfect_square (SCM k) +{ + int result; + + if (SCM_LIKELY (SCM_I_INUMP (k))) + { + mpz_t kk; + + mpz_init_set_si (kk, SCM_I_INUM (k)); + result = mpz_perfect_square_p (kk); + mpz_clear (kk); + } + else + { + result = mpz_perfect_square_p (SCM_I_BIG_MPZ (k)); + scm_remember_upto_here_1 (k); + } + return result; +} + +/* Return the floor of the square root of K. + K must be an exact integer. */ +static SCM +exact_integer_floor_square_root (SCM k) +{ + if (SCM_LIKELY (SCM_I_INUMP (k))) + { + mpz_t kk; + scm_t_inum ss; + + mpz_init_set_ui (kk, SCM_I_INUM (k)); + mpz_sqrt (kk, kk); + ss = mpz_get_ui (kk); + mpz_clear (kk); + return SCM_I_MAKINUM (ss); + } + else + { + SCM s; + + s = scm_i_mkbig (); + mpz_sqrt (SCM_I_BIG_MPZ (s), SCM_I_BIG_MPZ (k)); + scm_remember_upto_here_1 (k); + return scm_i_normbig (s); + } +} + SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, (SCM z), @@ -9850,14 +10084,114 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, } else if (SCM_NUMBERP (z)) { - double xx = scm_to_double (z); - if (xx < 0) - return scm_c_make_rectangular (0.0, sqrt (-xx)); - else - return scm_from_double (sqrt (xx)); + if (SCM_I_INUMP (z)) + { + scm_t_inum x = SCM_I_INUM (z); + + if (SCM_LIKELY (x >= 0)) + { + if (SCM_LIKELY (SCM_I_FIXNUM_BIT < DBL_MANT_DIG + || x < (1L << (DBL_MANT_DIG - 1)))) + { + double root = sqrt (x); + + /* If 0 <= x < 2^(DBL_MANT_DIG-1) and sqrt(x) is an + integer, then the result is exact. */ + if (root == floor (root)) + return SCM_I_MAKINUM ((scm_t_inum) root); + else + return scm_from_double (root); + } + else + { + mpz_t xx; + scm_t_inum root; + + mpz_init_set_ui (xx, x); + if (mpz_perfect_square_p (xx)) + { + mpz_sqrt (xx, xx); + root = mpz_get_ui (xx); + mpz_clear (xx); + return SCM_I_MAKINUM (root); + } + else + mpz_clear (xx); + } + } + } + else if (SCM_BIGP (z)) + { + if (mpz_perfect_square_p (SCM_I_BIG_MPZ (z))) + { + SCM root = scm_i_mkbig (); + + mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (z)); + scm_remember_upto_here_1 (z); + return scm_i_normbig (root); + } + else + { + long expon; + double signif = scm_i_big2dbl_2exp (z, &expon); + + if (expon & 1) + { + signif *= 2; + expon--; + } + if (signif < 0) + return scm_c_make_rectangular + (0.0, ldexp (sqrt (-signif), expon / 2)); + else + return scm_from_double (ldexp (sqrt (signif), expon / 2)); + } + } + else if (SCM_FRACTIONP (z)) + { + SCM n = SCM_FRACTION_NUMERATOR (z); + SCM d = SCM_FRACTION_DENOMINATOR (z); + + if (exact_integer_is_perfect_square (n) + && exact_integer_is_perfect_square (d)) + return scm_i_make_ratio_already_reduced + (exact_integer_floor_square_root (n), + exact_integer_floor_square_root (d)); + else + { + double xx = scm_i_divide2double (n, d); + double abs_xx = fabs (xx); + long shift = 0; + + if (SCM_UNLIKELY (abs_xx > DBL_MAX || abs_xx < DBL_MIN)) + { + shift = (scm_to_long (scm_integer_length (n)) + - scm_to_long (scm_integer_length (d))) / 2; + if (shift > 0) + d = left_shift_exact_integer (d, 2 * shift); + else + n = left_shift_exact_integer (n, -2 * shift); + xx = scm_i_divide2double (n, d); + } + + if (xx < 0) + return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift)); + else + return scm_from_double (ldexp (sqrt (xx), shift)); + } + } + + /* Fallback method, when the cases above do not apply. */ + { + double xx = scm_to_double (z); + if (xx < 0) + return scm_c_make_rectangular (0.0, sqrt (-xx)); + else + return scm_from_double (sqrt (xx)); + } } else - SCM_WTA_DISPATCH_1 (g_scm_sqrt, z, 1, s_scm_sqrt); + return scm_wta_dispatch_1 (g_scm_sqrt, z, 1, s_scm_sqrt); } #undef FUNC_NAME @@ -9866,8 +10200,6 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, void scm_init_numbers () { - int i; - if (scm_install_gmp_memory_functions) mp_set_memory_functions (custom_gmp_malloc, custom_gmp_realloc, @@ -9889,18 +10221,25 @@ scm_init_numbers () 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) - { - init_dblprec(&scm_dblprec[i-2],i); - init_fx_radix(fx_per_radix[i-2],i); - } -#ifdef DBL_DIG - /* hard code precision for base 10 if the preprocessor tells us to... */ - scm_dblprec[10-2] = (DBL_DIG > 20) ? 20 : DBL_DIG; -#endif - exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2)); + + { + /* Set scm_i_divide2double_lo2b to (2 b^p - 1) */ + mpz_init_set_ui (scm_i_divide2double_lo2b, 1); + mpz_mul_2exp (scm_i_divide2double_lo2b, + scm_i_divide2double_lo2b, + DBL_MANT_DIG + 1); /* 2 b^p */ + mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1); + } + + { + /* Set dbl_minimum_normal_mantissa to b^{p-1} */ + mpz_init_set_ui (dbl_minimum_normal_mantissa, 1); + mpz_mul_2exp (dbl_minimum_normal_mantissa, + dbl_minimum_normal_mantissa, + DBL_MANT_DIG - 1); + } + #include "libguile/numbers.x" }