mirror of
https://git.savannah.gnu.org/git/guile.git
synced 2025-04-29 19:30:36 +02:00
Fix scm_integer_to_double_z to always round; clean ups
* libguile/integers.c (scm_integer_to_double_z): Doubles that can't be exactly represented as integers should round. (bignum_frexp_z): New helper. (scm_integer_from_mpz, scm_integer_from_double): New internal functions. * libguile/numbers.h: * libguile/numbers.c (scm_i_bigcmp, scm_i_dbl2big, scm_i_dbl2num): Remove unused internal functions. (scm_inexact_to_exact): Rework to avoid scm_i_dbl2big. (scm_bigequal): Move here, from eq.c. (scm_integer_to_double_z): Use the new helper. (scm_i_big2dbl): Use scm_integer_to_double_z.
This commit is contained in:
Andy Wingo
parent
f2390e510f
commit
7c53325c31
5 changed files with 111 additions and 104 deletions
|
|
@ -1,4 +1,4 @@
|
|||
/* Copyright 1995-1998,2000-2001,2003-2004,2006,2009-2011,2017-2018
|
||||
/* Copyright 1995-1998,2000-2001,2003-2004,2006,2009-2011,2017-2018,2022
|
||||
Free Software Foundation, Inc.
|
||||
|
||||
This file is part of Guile.
|
||||
|
|
@ -131,12 +131,6 @@ scm_real_equalp (SCM x, SCM y)
|
|||
SCM_REAL_VALUE (y)));
|
||||
}
|
||||
|
||||
SCM
|
||||
scm_bigequal (SCM x, SCM y)
|
||||
{
|
||||
return scm_from_bool (scm_i_bigcmp (x, y) == 0);
|
||||
}
|
||||
|
||||
SCM
|
||||
scm_complex_equalp (SCM x, SCM y)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -286,6 +286,12 @@ bignum_cmp_long (struct scm_bignum *z, long l)
|
|||
}
|
||||
}
|
||||
|
||||
SCM
|
||||
scm_integer_from_mpz (mpz_srcptr mpz)
|
||||
{
|
||||
return normalize_bignum (make_bignum_from_mpz (mpz));
|
||||
}
|
||||
|
||||
int
|
||||
scm_is_integer_odd_i (scm_t_inum i)
|
||||
{
|
||||
|
|
@ -2503,14 +2509,68 @@ scm_is_integer_negative_z (struct scm_bignum *x)
|
|||
return bignum_is_negative (x);
|
||||
}
|
||||
|
||||
double
|
||||
scm_integer_to_double_z (struct scm_bignum *x)
|
||||
#if SCM_ENABLE_MINI_GMP
|
||||
static double
|
||||
mpz_get_d_2exp (long *exp, mpz_srcptr z)
|
||||
{
|
||||
double signif = mpz_get_d (z);
|
||||
int iexp;
|
||||
signif = frexp (signif, &iexp);
|
||||
*exp = iexp;
|
||||
return signif;
|
||||
}
|
||||
#endif
|
||||
|
||||
static double
|
||||
bignum_frexp (struct scm_bignum *x, long *exp)
|
||||
{
|
||||
mpz_t zx;
|
||||
alias_bignum_to_mpz (x, zx);
|
||||
double result = mpz_get_d (zx);
|
||||
|
||||
size_t bits = mpz_sizeinbase (zx, 2);
|
||||
ASSERT (bits != 0);
|
||||
size_t shift = 0;
|
||||
if (bits > DBL_MANT_DIG)
|
||||
{
|
||||
shift = bits - DBL_MANT_DIG;
|
||||
SCM xx = scm_integer_round_rsh_zu (x, shift);
|
||||
if (SCM_I_INUMP (xx))
|
||||
{
|
||||
int expon;
|
||||
double signif = frexp (SCM_I_INUM (xx), &expon);
|
||||
*exp = expon + shift;
|
||||
return signif;
|
||||
}
|
||||
x = scm_bignum (xx);
|
||||
alias_bignum_to_mpz (x, zx);
|
||||
}
|
||||
|
||||
double significand = mpz_get_d_2exp (exp, zx);
|
||||
scm_remember_upto_here_1 (x);
|
||||
return result;
|
||||
*exp += shift;
|
||||
return significand;
|
||||
}
|
||||
|
||||
double
|
||||
scm_integer_to_double_z (struct scm_bignum *x)
|
||||
{
|
||||
long exponent;
|
||||
double significand = bignum_frexp (x, &exponent);
|
||||
return ldexp (significand, exponent);
|
||||
}
|
||||
|
||||
SCM
|
||||
scm_integer_from_double (double val)
|
||||
{
|
||||
if (!isfinite (val))
|
||||
scm_out_of_range ("inexact->exact", scm_from_double (val));
|
||||
|
||||
if (((double) INT64_MIN) <= val && val <= ((double) INT64_MAX))
|
||||
return scm_from_int64 (val);
|
||||
|
||||
mpz_t result;
|
||||
mpz_init_set_d (result, val);
|
||||
return take_mpz (result);
|
||||
}
|
||||
|
||||
SCM
|
||||
|
|
|
|||
|
|
@ -32,6 +32,8 @@ scm_bignum (SCM x)
|
|||
return (struct scm_bignum *) SCM_UNPACK (x);
|
||||
}
|
||||
|
||||
SCM_INTERNAL SCM scm_integer_from_mpz (mpz_srcptr mpz);
|
||||
|
||||
SCM_INTERNAL int scm_is_integer_odd_i (scm_t_inum i);
|
||||
SCM_INTERNAL int scm_is_integer_odd_z (struct scm_bignum *z);
|
||||
|
||||
|
|
@ -167,6 +169,7 @@ SCM_INTERNAL int scm_is_integer_positive_z (struct scm_bignum *x);
|
|||
SCM_INTERNAL int scm_is_integer_negative_z (struct scm_bignum *x);
|
||||
|
||||
SCM_INTERNAL double scm_integer_to_double_z (struct scm_bignum *x);
|
||||
SCM_INTERNAL SCM scm_integer_from_double (double val);
|
||||
|
||||
SCM_INTERNAL SCM scm_integer_add_ii (scm_t_inum x, scm_t_inum y);
|
||||
SCM_INTERNAL SCM scm_integer_add_zi (struct scm_bignum *x, scm_t_inum y);
|
||||
|
|
|
|||
|
|
@ -338,51 +338,6 @@ scm_i_clonebig (SCM src_big, int same_sign_p)
|
|||
return z;
|
||||
}
|
||||
|
||||
int
|
||||
scm_i_bigcmp (SCM x, SCM y)
|
||||
{
|
||||
/* Return neg if x < y, pos if x > y, and 0 if x == y */
|
||||
/* presume we already know x and y are bignums */
|
||||
int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
|
||||
scm_remember_upto_here_2 (x, y);
|
||||
return result;
|
||||
}
|
||||
|
||||
SCM
|
||||
scm_i_dbl2big (double d)
|
||||
{
|
||||
/* results are only defined if d is an integer */
|
||||
SCM z = make_bignum ();
|
||||
mpz_init_set_d (SCM_I_BIG_MPZ (z), d);
|
||||
return z;
|
||||
}
|
||||
|
||||
/* Convert a integer in double representation to a SCM number. */
|
||||
|
||||
SCM
|
||||
scm_i_dbl2num (double u)
|
||||
{
|
||||
/* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
|
||||
powers of 2, so there's no rounding when making "double" values
|
||||
from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
|
||||
get rounded on a 64-bit machine, hence the "+1".
|
||||
|
||||
The use of floor() to force to an integer value ensures we get a
|
||||
"numerically closest" value without depending on how a
|
||||
double->long cast or how mpz_set_d will round. For reference,
|
||||
double->long probably follows the hardware rounding mode,
|
||||
mpz_set_d truncates towards zero. */
|
||||
|
||||
/* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
|
||||
representable as a double? */
|
||||
|
||||
if (u < (double) (SCM_MOST_POSITIVE_FIXNUM+1)
|
||||
&& u >= (double) SCM_MOST_NEGATIVE_FIXNUM)
|
||||
return SCM_I_MAKINUM ((scm_t_inum) u);
|
||||
else
|
||||
return scm_i_dbl2big (u);
|
||||
}
|
||||
|
||||
static SCM round_rsh (SCM n, SCM count);
|
||||
|
||||
/* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the
|
||||
|
|
@ -434,9 +389,7 @@ scm_i_big2dbl_2exp (SCM b, long *expon_p)
|
|||
double
|
||||
scm_i_big2dbl (SCM b)
|
||||
{
|
||||
long expon;
|
||||
double signif = scm_i_big2dbl_2exp (b, &expon);
|
||||
return ldexp (signif, expon);
|
||||
return scm_integer_to_double_z (scm_bignum (b));
|
||||
}
|
||||
|
||||
SCM
|
||||
|
|
@ -4619,6 +4572,12 @@ SCM_DEFINE (scm_exact_integer_p, "exact-integer?", 1, 0, 0,
|
|||
}
|
||||
#undef FUNC_NAME
|
||||
|
||||
SCM
|
||||
scm_bigequal (SCM x, SCM y)
|
||||
{
|
||||
return scm_from_bool
|
||||
(scm_is_integer_equal_zz (scm_bignum (x), scm_bignum (y)));
|
||||
}
|
||||
|
||||
SCM scm_i_num_eq_p (SCM, SCM, SCM);
|
||||
SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p, "=", 0, 2, 1,
|
||||
|
|
@ -6561,51 +6520,45 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
|
|||
{
|
||||
if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
|
||||
return z;
|
||||
|
||||
double val;
|
||||
|
||||
if (SCM_REALP (z))
|
||||
val = SCM_REAL_VALUE (z);
|
||||
else if (SCM_COMPLEXP (z) && SCM_COMPLEX_IMAG (z) == 0.0)
|
||||
val = SCM_COMPLEX_REAL (z);
|
||||
else
|
||||
return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1,
|
||||
s_scm_inexact_to_exact);
|
||||
|
||||
if (!SCM_LIKELY (isfinite (val)))
|
||||
SCM_OUT_OF_RANGE (1, z);
|
||||
if (val == 0)
|
||||
return SCM_INUM0;
|
||||
|
||||
int expon;
|
||||
mpz_t zn;
|
||||
|
||||
mpz_init_set_d (zn, ldexp (frexp (val, &expon), DBL_MANT_DIG));
|
||||
expon -= DBL_MANT_DIG;
|
||||
if (expon < 0)
|
||||
{
|
||||
double val;
|
||||
int shift = mpz_scan1 (zn, 0);
|
||||
|
||||
if (SCM_REALP (z))
|
||||
val = SCM_REAL_VALUE (z);
|
||||
else if (SCM_COMPLEXP (z) && SCM_COMPLEX_IMAG (z) == 0.0)
|
||||
val = SCM_COMPLEX_REAL (z);
|
||||
else
|
||||
return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1,
|
||||
s_scm_inexact_to_exact);
|
||||
|
||||
if (!SCM_LIKELY (isfinite (val)))
|
||||
SCM_OUT_OF_RANGE (1, z);
|
||||
else if (val == 0.0)
|
||||
return SCM_INUM0;
|
||||
else
|
||||
{
|
||||
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, scm_integer_lsh_iu (1, -expon));
|
||||
else if (expon > 0)
|
||||
return lsh (numerator, scm_from_int (expon), FUNC_NAME);
|
||||
else
|
||||
return numerator;
|
||||
}
|
||||
if (shift > -expon)
|
||||
shift = -expon;
|
||||
mpz_fdiv_q_2exp (zn, zn, shift);
|
||||
expon += shift;
|
||||
}
|
||||
SCM numerator = scm_integer_from_mpz (zn);
|
||||
mpz_clear (zn);
|
||||
if (expon < 0)
|
||||
return scm_i_make_ratio_already_reduced
|
||||
(numerator, scm_integer_lsh_iu (1, -expon));
|
||||
else if (expon > 0)
|
||||
return lsh (numerator, scm_from_int (expon), FUNC_NAME);
|
||||
else
|
||||
return numerator;
|
||||
}
|
||||
#undef FUNC_NAME
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
#ifndef SCM_NUMBERS_H
|
||||
#define SCM_NUMBERS_H
|
||||
|
||||
/* Copyright 1995-1996,1998,2000-2006,2008-2011,2013-2014,2016-2018,2021
|
||||
/* Copyright 1995-1996,1998,2000-2006,2008-2011,2013-2014,2016-2018,2021,2022
|
||||
Free Software Foundation, Inc.
|
||||
|
||||
This file is part of Guile.
|
||||
|
|
@ -349,9 +349,6 @@ SCM_INTERNAL SCM scm_i_exact_integer_sqrt (SCM k);
|
|||
/* bignum internal functions */
|
||||
SCM_INTERNAL SCM scm_i_mkbig (void);
|
||||
SCM_API /* FIXME: not internal */ SCM scm_i_normbig (SCM x);
|
||||
SCM_INTERNAL int scm_i_bigcmp (SCM a, SCM b);
|
||||
SCM_INTERNAL SCM scm_i_dbl2big (double d);
|
||||
SCM_INTERNAL SCM scm_i_dbl2num (double d);
|
||||
SCM_API /* FIXME: not internal */ double scm_i_big2dbl (SCM b);
|
||||
SCM_API /* FIXME: not internal */ SCM scm_i_long2big (long n);
|
||||
SCM_API /* FIXME: not internal */ SCM scm_i_ulong2big (unsigned long n);
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue