diff --git a/libguile/eq.c b/libguile/eq.c index 5d8c19a71..6870613c1 100644 --- a/libguile/eq.c +++ b/libguile/eq.c @@ -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) { diff --git a/libguile/integers.c b/libguile/integers.c index 4b26ea3e0..2e35bc2d5 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -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 diff --git a/libguile/integers.h b/libguile/integers.h index 98c19b90f..bda575774 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -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); diff --git a/libguile/numbers.c b/libguile/numbers.c index 17131b9ba..549b730ec 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -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 diff --git a/libguile/numbers.h b/libguile/numbers.h index df5c9110c..f5cfe221e 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -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);