[Guile-commits] 48/69: Fix scm_integer_to_double_z to always round; clean ups

[Andy Wingo]

wingo pushed a commit to branch wip-inline-digits
in repository guile.

commit e79687477708aa373e85dd04fafbbd344ae238b3
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Wed Jan 5 21:12:01 2022 +0100

    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.
---
 libguile/eq.c       |   8 +---
 libguile/integers.c |  68 +++++++++++++++++++++++++--
 libguile/integers.h |   3 ++
 libguile/numbers.c  | 129 +++++++++++++++++-----------------------------------
 libguile/numbers.h  |   5 +-
 5 files changed, 110 insertions(+), 103 deletions(-)

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
-    {
-      double val;
+    return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1,
+                               s_scm_inexact_to_exact);
 
-      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;
 
-      if (!SCM_LIKELY (isfinite (val)))
-       SCM_OUT_OF_RANGE (1, z);
-      else if (val == 0.0)
-        return SCM_INUM0;
-      else
-       {
-          int expon;
-          SCM numerator;
+  int expon;
+  mpz_t zn;
 
-          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;
-       }
+  mpz_init_set_d (zn, ldexp (frexp (val, &expon), DBL_MANT_DIG));
+  expon -= DBL_MANT_DIG;
+  if (expon < 0)
+    {
+      int shift = mpz_scan1 (zn, 0);
+
+      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);