Merge remote-tracking branch 'origin/stable-2.0'
[bpt/guile.git] / libguile / numbers.c
index fa55b4f..9857e18 100644 (file)
@@ -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 <complex.h>
 #endif
 
+#include <stdarg.h>
+
 #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
+
 \f
 
 /*
@@ -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<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" */
@@ -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 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;
@@ -6645,7 +6877,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_FRACTIONP (x))
     {
@@ -6678,10 +6911,12 @@ 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
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARG1, s_scm_i_num_less_p);
+    return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARG1,
+                               s_scm_i_num_less_p);
 }
 
 
@@ -6710,9 +6945,9 @@ SCM
 scm_gr_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p, x, y, SCM_ARG1, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_gr_p, x, y, SCM_ARG1, FUNC_NAME);
   else if (!SCM_NUMBERP (y))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p, x, y, SCM_ARG2, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_gr_p, x, y, SCM_ARG2, FUNC_NAME);
   else
     return scm_less_p (y, x);
 }
@@ -6744,9 +6979,9 @@ SCM
 scm_leq_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p, x, y, SCM_ARG1, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_leq_p, x, y, SCM_ARG1, FUNC_NAME);
   else if (!SCM_NUMBERP (y))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p, x, y, SCM_ARG2, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_leq_p, x, y, SCM_ARG2, FUNC_NAME);
   else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
     return SCM_BOOL_F;
   else
@@ -6780,9 +7015,9 @@ SCM
 scm_geq_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p, x, y, SCM_ARG1, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_geq_p, x, y, SCM_ARG1, FUNC_NAME);
   else if (!SCM_NUMBERP (y))
-    SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p, x, y, SCM_ARG2, FUNC_NAME);
+    return scm_wta_dispatch_2 (g_scm_i_num_geq_p, x, y, SCM_ARG2, FUNC_NAME);
   else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
     return SCM_BOOL_F;
   else
@@ -6809,7 +7044,7 @@ SCM_PRIMITIVE_GENERIC (scm_zero_p, "zero?", 1, 0, 0,
   else if (SCM_FRACTIONP (z))
     return SCM_BOOL_F;
   else
-    SCM_WTA_DISPATCH_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p);
+    return scm_wta_dispatch_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p);
 }
 #undef FUNC_NAME
 
@@ -6833,7 +7068,7 @@ SCM_PRIMITIVE_GENERIC (scm_positive_p, "positive?", 1, 0, 0,
   else if (SCM_FRACTIONP (x))
     return scm_positive_p (SCM_FRACTION_NUMERATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p);
+    return scm_wta_dispatch_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p);
 }
 #undef FUNC_NAME
 
@@ -6857,7 +7092,7 @@ SCM_PRIMITIVE_GENERIC (scm_negative_p, "negative?", 1, 0, 0,
   else if (SCM_FRACTIONP (x))
     return scm_negative_p (SCM_FRACTION_NUMERATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p);
+    return scm_wta_dispatch_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p);
 }
 #undef FUNC_NAME
 
@@ -6891,11 +7126,11 @@ scm_max (SCM x, SCM y)
   if (SCM_UNBNDP (y))
     {
       if (SCM_UNBNDP (x))
-       SCM_WTA_DISPATCH_0 (g_max, s_max);
+       return scm_wta_dispatch_0 (g_max, s_max);
       else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
        return x;
       else
-       SCM_WTA_DISPATCH_1 (g_max, x, SCM_ARG1, s_max);
+       return scm_wta_dispatch_1 (g_max, x, SCM_ARG1, s_max);
     }
   
   if (SCM_I_INUMP (x))
@@ -6934,7 +7169,7 @@ scm_max (SCM x, SCM y)
           return (scm_is_false (scm_less_p (x, y)) ? x : y);
        }
       else
-       SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
+       return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
     }
   else if (SCM_BIGP (x))
     {
@@ -6964,7 +7199,7 @@ scm_max (SCM x, SCM y)
           goto use_less;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
+       return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
     }
   else if (SCM_REALP (x))
     {
@@ -7019,7 +7254,7 @@ scm_max (SCM x, SCM y)
          return (xx < yy) ? scm_from_double (yy) : x;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
+       return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -7042,10 +7277,10 @@ scm_max (SCM x, SCM y)
           goto use_less;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
+       return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARG1, s_max);
+    return scm_wta_dispatch_2 (g_max, x, y, SCM_ARG1, s_max);
 }
 
 
@@ -7072,11 +7307,11 @@ scm_min (SCM x, SCM y)
   if (SCM_UNBNDP (y))
     {
       if (SCM_UNBNDP (x))
-       SCM_WTA_DISPATCH_0 (g_min, s_min);
+       return scm_wta_dispatch_0 (g_min, s_min);
       else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
        return x;
       else
-       SCM_WTA_DISPATCH_1 (g_min, x, SCM_ARG1, s_min);
+       return scm_wta_dispatch_1 (g_min, x, SCM_ARG1, s_min);
     }
   
   if (SCM_I_INUMP (x))
@@ -7105,7 +7340,7 @@ scm_min (SCM x, SCM y)
           return (scm_is_false (scm_less_p (x, y)) ? y : x);
        }
       else
-       SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
+       return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
     }
   else if (SCM_BIGP (x))
     {
@@ -7135,7 +7370,7 @@ scm_min (SCM x, SCM y)
           goto use_less;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
+       return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
     }
   else if (SCM_REALP (x))
     {
@@ -7179,7 +7414,7 @@ scm_min (SCM x, SCM y)
          return (yy < xx) ? scm_from_double (yy) : x;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
+       return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -7202,10 +7437,10 @@ scm_min (SCM x, SCM y)
           goto use_less;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
+       return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARG1, s_min);
+    return scm_wta_dispatch_2 (g_min, x, y, SCM_ARG1, s_min);
 }
 
 
@@ -7234,7 +7469,7 @@ scm_sum (SCM x, SCM y)
     {
       if (SCM_NUMBERP (x)) return x;
       if (SCM_UNBNDP (x)) return SCM_INUM0;
-      SCM_WTA_DISPATCH_1 (g_sum, x, SCM_ARG1, s_sum);
+      return scm_wta_dispatch_1 (g_sum, x, SCM_ARG1, s_sum);
     }
 
   if (SCM_LIKELY (SCM_I_INUMP (x)))
@@ -7267,7 +7502,7 @@ scm_sum (SCM x, SCM y)
                                        scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
                               SCM_FRACTION_DENOMINATOR (y));
       else
-        SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
+        return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
     } else if (SCM_BIGP (x))
       {
        if (SCM_I_INUMP (y))
@@ -7332,7 +7567,7 @@ scm_sum (SCM x, SCM y)
                                          scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
                                 SCM_FRACTION_DENOMINATOR (y));
        else
-         SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
+         return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
       }
   else if (SCM_REALP (x))
     {
@@ -7352,7 +7587,7 @@ scm_sum (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
       else
-       SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
+       return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
     }
   else if (SCM_COMPLEXP (x))
     {
@@ -7376,7 +7611,7 @@ scm_sum (SCM x, SCM y)
        return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + scm_i_fraction2double (y),
                                 SCM_COMPLEX_IMAG (x));
       else
-       SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
+       return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -7399,10 +7634,10 @@ scm_sum (SCM x, SCM y)
                                        scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
                               scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
       else
-       SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
+       return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARG1, s_sum);
+    return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARG1, s_sum);
 }
 
 
@@ -7442,7 +7677,7 @@ scm_difference (SCM x, SCM y)
   if (SCM_UNLIKELY (SCM_UNBNDP (y)))
     {
       if (SCM_UNBNDP (x))
-        SCM_WTA_DISPATCH_0 (g_difference, s_difference);
+        return scm_wta_dispatch_0 (g_difference, s_difference);
       else 
         if (SCM_I_INUMP (x))
           {
@@ -7466,7 +7701,7 @@ scm_difference (SCM x, SCM y)
            (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
             SCM_FRACTION_DENOMINATOR (x));
         else
-          SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference);
+          return scm_wta_dispatch_1 (g_difference, x, SCM_ARG1, s_difference);
     }
   
   if (SCM_LIKELY (SCM_I_INUMP (x)))
@@ -7553,7 +7788,7 @@ scm_difference (SCM x, SCM y)
                                               SCM_FRACTION_NUMERATOR (y)),
                               SCM_FRACTION_DENOMINATOR (y));
       else
-       SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
+       return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
   else if (SCM_BIGP (x))
     {
@@ -7617,7 +7852,8 @@ scm_difference (SCM x, SCM y)
        return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
                                               SCM_FRACTION_NUMERATOR (y)),
                               SCM_FRACTION_DENOMINATOR (y));
-      else SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
+      else
+        return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
   else if (SCM_REALP (x))
     {
@@ -7637,7 +7873,7 @@ scm_difference (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
       else
-       SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
+       return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
   else if (SCM_COMPLEXP (x))
     {
@@ -7661,7 +7897,7 @@ scm_difference (SCM x, SCM y)
        return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - scm_i_fraction2double (y),
                                 SCM_COMPLEX_IMAG (x));
       else
-       SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
+       return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -7685,10 +7921,10 @@ scm_difference (SCM x, SCM y)
                                               scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
                               scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
       else
-       SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
+       return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARG1, s_difference);
+    return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARG1, s_difference);
 }
 #undef FUNC_NAME
 
@@ -7731,7 +7967,7 @@ scm_product (SCM x, SCM y)
       else if (SCM_NUMBERP (x))
        return x;
       else
-       SCM_WTA_DISPATCH_1 (g_product, x, SCM_ARG1, s_product);
+       return scm_wta_dispatch_1 (g_product, x, SCM_ARG1, s_product);
     }
   
   if (SCM_LIKELY (SCM_I_INUMP (x)))
@@ -7764,7 +8000,7 @@ scm_product (SCM x, SCM y)
          else if (SCM_NUMP (y))
            return SCM_INUM0;
          else
-           SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+           return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
          break;
         case -1:
          /*
@@ -7815,7 +8051,7 @@ scm_product (SCM x, SCM y)
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
                               SCM_FRACTION_DENOMINATOR (y));
       else
-       SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
   else if (SCM_BIGP (x))
     {
@@ -7850,7 +8086,7 @@ scm_product (SCM x, SCM y)
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
                               SCM_FRACTION_DENOMINATOR (y));
       else
-       SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
   else if (SCM_REALP (x))
     {
@@ -7873,7 +8109,7 @@ scm_product (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
       else
-       SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
   else if (SCM_COMPLEXP (x))
     {
@@ -7906,7 +8142,7 @@ scm_product (SCM x, SCM y)
                                   yy * SCM_COMPLEX_IMAG (x));
        }
       else
-       SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -7931,10 +8167,10 @@ scm_product (SCM x, SCM y)
                               scm_product (SCM_FRACTION_DENOMINATOR (x),
                                            SCM_FRACTION_DENOMINATOR (y)));
       else
-       SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
+       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARG1, s_product);
+    return scm_wta_dispatch_2 (g_product, x, y, SCM_ARG1, s_product);
 }
 
 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
@@ -7989,8 +8225,8 @@ SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1,
 #define s_divide s_scm_i_divide
 #define g_divide g_scm_i_divide
 
-static SCM
-do_divide (SCM x, SCM y, int inexact)
+SCM
+scm_divide (SCM x, SCM y)
 #define FUNC_NAME s_divide
 {
   double a;
@@ -7998,7 +8234,7 @@ do_divide (SCM x, SCM y, int inexact)
   if (SCM_UNLIKELY (SCM_UNBNDP (y)))
     {
       if (SCM_UNBNDP (x))
-       SCM_WTA_DISPATCH_0 (g_divide, s_divide);
+       return scm_wta_dispatch_0 (g_divide, s_divide);
       else if (SCM_I_INUMP (x))
        {
          scm_t_inum xx = SCM_I_INUM (x);
@@ -8009,18 +8245,10 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
 #endif
          else
-           {
-             if (inexact)
-               return scm_from_double (1.0 / (double) xx);
-             else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
-           }
+           return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
        }
       else if (SCM_BIGP (x))
-       {
-         if (inexact)
-           return scm_from_double (1.0 / scm_i_big2dbl (x));
-         else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
-       }
+       return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
       else if (SCM_REALP (x))
        {
          double xx = SCM_REAL_VALUE (x);
@@ -8052,7 +8280,7 @@ do_divide (SCM x, SCM y, int inexact)
        return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x),
                                                 SCM_FRACTION_NUMERATOR (x));
       else
-       SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
+       return scm_wta_dispatch_1 (g_divide, x, SCM_ARG1, s_divide);
     }
 
   if (SCM_LIKELY (SCM_I_INUMP (x)))
@@ -8070,11 +8298,7 @@ do_divide (SCM x, SCM y, int inexact)
 #endif
            }
          else if (xx % yy != 0)
-           {
-             if (inexact)
-               return scm_from_double ((double) xx / (double) yy);
-             else return scm_i_make_ratio (x, y);
-           }
+           return scm_i_make_ratio (x, y);
          else
            {
              scm_t_inum z = xx / yy;
@@ -8085,11 +8309,7 @@ do_divide (SCM x, SCM y, int inexact)
            }
        }
       else if (SCM_BIGP (y))
-       {
-         if (inexact)
-           return scm_from_double ((double) xx / scm_i_big2dbl (y));
-         else return scm_i_make_ratio (x, y);
-       }
+       return scm_i_make_ratio (x, y);
       else if (SCM_REALP (y))
        {
          double yy = SCM_REAL_VALUE (y);
@@ -8098,6 +8318,9 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* 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);
        }
       else if (SCM_COMPLEXP (y))
@@ -8124,9 +8347,9 @@ do_divide (SCM x, SCM y, int inexact)
       else if (SCM_FRACTIONP (y))
        /* a / b/c = ac / b */
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
-                              SCM_FRACTION_NUMERATOR (y));
+                                 SCM_FRACTION_NUMERATOR (y));
       else
-       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
+       return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else if (SCM_BIGP (x))
     {
@@ -8168,43 +8391,24 @@ do_divide (SCM x, SCM y, int inexact)
                  return scm_i_normbig (result);
                }
              else
-               {
-                 if (inexact)
-                   return scm_from_double (scm_i_big2dbl (x) / (double) yy);
-                 else return scm_i_make_ratio (x, y);
-               }
+               return scm_i_make_ratio (x, y);
            }
        }
       else if (SCM_BIGP (y))
        {
-         /* big_x / big_y */
-         if (inexact)
-           {
-             /* It's easily possible for the ratio x/y to fit a double
-                but one or both x and y be too big to fit a double,
-                hence the use of mpq_get_d rather than converting and
-                dividing.  */
-             mpq_t q;
-             *mpq_numref(q) = *SCM_I_BIG_MPZ (x);
-             *mpq_denref(q) = *SCM_I_BIG_MPZ (y);
-             return scm_from_double (mpq_get_d (q));
-           }
-         else
-           {
-             int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
-                                                SCM_I_BIG_MPZ (y));
-             if (divisible_p)
-               {
-                 SCM result = scm_i_mkbig ();
-                 mpz_divexact (SCM_I_BIG_MPZ (result),
-                               SCM_I_BIG_MPZ (x),
-                               SCM_I_BIG_MPZ (y));
-                 scm_remember_upto_here_2 (x, y);
-                 return scm_i_normbig (result);
-               }
-             else
-               return scm_i_make_ratio (x, y);
-           }
+          int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
+                                             SCM_I_BIG_MPZ (y));
+          if (divisible_p)
+            {
+              SCM result = scm_i_mkbig ();
+              mpz_divexact (SCM_I_BIG_MPZ (result),
+                            SCM_I_BIG_MPZ (x),
+                            SCM_I_BIG_MPZ (y));
+              scm_remember_upto_here_2 (x, y);
+              return scm_i_normbig (result);
+            }
+          else
+            return scm_i_make_ratio (x, y);
        }
       else if (SCM_REALP (y))
        {
@@ -8214,6 +8418,8 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #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);
        }
       else if (SCM_COMPLEXP (y))
@@ -8223,9 +8429,9 @@ do_divide (SCM x, SCM y, int inexact)
        }
       else if (SCM_FRACTIONP (y))
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
-                              SCM_FRACTION_NUMERATOR (y));
+                                 SCM_FRACTION_NUMERATOR (y));
       else
-       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
+       return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else if (SCM_REALP (x))
     {
@@ -8238,10 +8444,16 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* 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);
        }
       else if (SCM_BIGP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from bignum to double
+             (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);
@@ -8264,7 +8476,7 @@ do_divide (SCM x, SCM y, int inexact)
       else if (SCM_FRACTIONP (y))
        return scm_from_double (rx / scm_i_fraction2double (y));
       else
-       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
+       return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else if (SCM_COMPLEXP (x))
     {
@@ -8279,12 +8491,18 @@ do_divide (SCM x, SCM y, int inexact)
          else
 #endif
            {
+              /* FIXME: Precision may be lost here due to:
+                 (1) The conversion from 'scm_t_inum' to double
+                 (2) Double rounding */
              double d = yy;
              return scm_c_make_rectangular (rx / d, ix / d);
            }
        }
       else if (SCM_BIGP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from bignum to double
+             (2) Double rounding */
          double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
          scm_remember_upto_here_1 (y);
          return scm_c_make_rectangular (rx / dby, ix / dby);
@@ -8318,11 +8536,14 @@ do_divide (SCM x, SCM y, int inexact)
        }
       else if (SCM_FRACTIONP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from fraction to double
+             (2) Double rounding */
          double yy = scm_i_fraction2double (y);
          return scm_c_make_rectangular (rx / yy, ix / yy);
        }
       else
-       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
+       return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else if (SCM_FRACTIONP (x))
     {
@@ -8335,12 +8556,12 @@ do_divide (SCM x, SCM y, int inexact)
          else
 #endif
            return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
-                                  scm_product (SCM_FRACTION_DENOMINATOR (x), y));
+                                     scm_product (SCM_FRACTION_DENOMINATOR (x), y));
        } 
       else if (SCM_BIGP (y)) 
        {
          return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
-                                scm_product (SCM_FRACTION_DENOMINATOR (x), y));
+                                   scm_product (SCM_FRACTION_DENOMINATOR (x), y));
        } 
       else if (SCM_REALP (y)) 
        {
@@ -8350,32 +8571,27 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* 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);
        }
       else if (SCM_COMPLEXP (y)) 
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from fraction to double
+             (2) Double rounding */
          a = scm_i_fraction2double (x);
          goto complex_div;
        } 
       else if (SCM_FRACTIONP (y))
        return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
-                              scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
+                                 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
       else 
-       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
+       return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide);
-}
-
-SCM
-scm_divide (SCM x, SCM y)
-{
-  return do_divide (x, y, 0);
-}
-
-static SCM scm_divide2real (SCM x, SCM y)
-{
-  return do_divide (x, y, 1);
+    return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARG1, s_divide);
 }
 #undef FUNC_NAME
 
@@ -8441,7 +8657,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_number, "truncate", 1, 0, 0,
     return scm_truncate_quotient (SCM_FRACTION_NUMERATOR (x),
                                  SCM_FRACTION_DENOMINATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_truncate_number, x, SCM_ARG1,
+    return scm_wta_dispatch_1 (g_scm_truncate_number, x, SCM_ARG1,
                        s_scm_truncate_number);
 }
 #undef FUNC_NAME
@@ -8461,8 +8677,8 @@ SCM_PRIMITIVE_GENERIC (scm_round_number, "round", 1, 0, 0,
     return scm_round_quotient (SCM_FRACTION_NUMERATOR (x),
                               SCM_FRACTION_DENOMINATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_round_number, x, SCM_ARG1,
-                       s_scm_round_number);
+    return scm_wta_dispatch_1 (g_scm_round_number, x, SCM_ARG1,
+                               s_scm_round_number);
 }
 #undef FUNC_NAME
 
@@ -8479,7 +8695,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0,
     return scm_floor_quotient (SCM_FRACTION_NUMERATOR (x),
                               SCM_FRACTION_DENOMINATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_floor, x, 1, s_scm_floor);
+    return scm_wta_dispatch_1 (g_scm_floor, x, 1, s_scm_floor);
 }  
 #undef FUNC_NAME
 
@@ -8496,7 +8712,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
     return scm_ceiling_quotient (SCM_FRACTION_NUMERATOR (x),
                                 SCM_FRACTION_DENOMINATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_ceiling, x, 1, s_scm_ceiling);
+    return scm_wta_dispatch_1 (g_scm_ceiling, x, 1, s_scm_ceiling);
 }
 #undef FUNC_NAME
 
@@ -8535,9 +8751,9 @@ SCM_PRIMITIVE_GENERIC (scm_expt, "expt", 2, 0, 0,
   else if (scm_is_complex (x) && scm_is_complex (y))
     return scm_exp (scm_product (scm_log (x), y));
   else if (scm_is_complex (x))
-    SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG2, s_scm_expt);
+    return scm_wta_dispatch_2 (g_scm_expt, x, y, SCM_ARG2, s_scm_expt);
   else
-    SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt);
+    return scm_wta_dispatch_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt);
 }
 #undef FUNC_NAME
 
@@ -8564,7 +8780,7 @@ SCM_PRIMITIVE_GENERIC (scm_sin, "sin", 1, 0, 0,
                                      cos (x) * sinh (y));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_sin, z, 1, s_scm_sin);
+    return scm_wta_dispatch_1 (g_scm_sin, z, 1, s_scm_sin);
 }
 #undef FUNC_NAME
 
@@ -8585,7 +8801,7 @@ SCM_PRIMITIVE_GENERIC (scm_cos, "cos", 1, 0, 0,
                                      -sin (x) * sinh (y));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_cos, z, 1, s_scm_cos);
+    return scm_wta_dispatch_1 (g_scm_cos, z, 1, s_scm_cos);
 }
 #undef FUNC_NAME
 
@@ -8610,7 +8826,7 @@ SCM_PRIMITIVE_GENERIC (scm_tan, "tan", 1, 0, 0,
       return scm_c_make_rectangular (sin (x) / w, sinh (y) / w);
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_tan, z, 1, s_scm_tan);
+    return scm_wta_dispatch_1 (g_scm_tan, z, 1, s_scm_tan);
 }
 #undef FUNC_NAME
 
@@ -8631,7 +8847,7 @@ SCM_PRIMITIVE_GENERIC (scm_sinh, "sinh", 1, 0, 0,
                                      cosh (x) * sin (y));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_sinh, z, 1, s_scm_sinh);
+    return scm_wta_dispatch_1 (g_scm_sinh, z, 1, s_scm_sinh);
 }
 #undef FUNC_NAME
 
@@ -8652,7 +8868,7 @@ SCM_PRIMITIVE_GENERIC (scm_cosh, "cosh", 1, 0, 0,
                                      sinh (x) * sin (y));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_cosh, z, 1, s_scm_cosh);
+    return scm_wta_dispatch_1 (g_scm_cosh, z, 1, s_scm_cosh);
 }
 #undef FUNC_NAME
 
@@ -8677,7 +8893,7 @@ SCM_PRIMITIVE_GENERIC (scm_tanh, "tanh", 1, 0, 0,
       return scm_c_make_rectangular (sinh (x) / w, sin (y) / w);
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_tanh, z, 1, s_scm_tanh);
+    return scm_wta_dispatch_1 (g_scm_tanh, z, 1, s_scm_tanh);
 }
 #undef FUNC_NAME
 
@@ -8705,7 +8921,7 @@ SCM_PRIMITIVE_GENERIC (scm_asin, "asin", 1, 0, 0,
                           scm_sys_asinh (scm_c_make_rectangular (-y, x)));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_asin, z, 1, s_scm_asin);
+    return scm_wta_dispatch_1 (g_scm_asin, z, 1, s_scm_asin);
 }
 #undef FUNC_NAME
 
@@ -8735,7 +8951,7 @@ SCM_PRIMITIVE_GENERIC (scm_acos, "acos", 1, 0, 0,
                                    scm_sys_asinh (scm_c_make_rectangular (-y, x))));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_acos, z, 1, s_scm_acos);
+    return scm_wta_dispatch_1 (g_scm_acos, z, 1, s_scm_acos);
 }
 #undef FUNC_NAME
 
@@ -8762,17 +8978,17 @@ SCM_PRIMITIVE_GENERIC (scm_atan, "atan", 1, 1, 0,
                              scm_c_make_rectangular (0, 2));
         }
       else
-        SCM_WTA_DISPATCH_1 (g_scm_atan, z, SCM_ARG1, s_scm_atan);
+        return scm_wta_dispatch_1 (g_scm_atan, z, SCM_ARG1, s_scm_atan);
     }
   else if (scm_is_real (z))
     {
       if (scm_is_real (y))
         return scm_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
       else
-        SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan);
+        return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan);
+    return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan);
 }
 #undef FUNC_NAME
 
@@ -8790,7 +9006,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_asinh, "asinh", 1, 0, 0,
                              scm_sqrt (scm_sum (scm_product (z, z),
                                                 SCM_INUM1))));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh);
+    return scm_wta_dispatch_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh);
 }
 #undef FUNC_NAME
 
@@ -8808,7 +9024,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_acosh, "acosh", 1, 0, 0,
                              scm_sqrt (scm_difference (scm_product (z, z),
                                                        SCM_INUM1))));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh);
+    return scm_wta_dispatch_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh);
 }
 #undef FUNC_NAME
 
@@ -8826,7 +9042,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_atanh, "atanh", 1, 0, 0,
                                             scm_difference (SCM_INUM1, z))),
                        SCM_I_MAKINUM (2));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_sys_atanh, z, 1, s_scm_sys_atanh);
+    return scm_wta_dispatch_1 (g_scm_sys_atanh, z, 1, s_scm_sys_atanh);
 }
 #undef FUNC_NAME
 
@@ -8835,7 +9051,7 @@ scm_c_make_rectangular (double re, double im)
 {
   SCM z;
 
-  z = PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex),
+  z = SCM_PACK_POINTER (scm_gc_malloc_pointerless (sizeof (scm_t_complex),
                                          "complex"));
   SCM_SET_CELL_TYPE (z, scm_tc16_complex);
   SCM_COMPLEX_REAL (z) = re;
@@ -8927,7 +9143,7 @@ SCM_PRIMITIVE_GENERIC (scm_real_part, "real-part", 1, 0, 0,
   else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z))
     return z;
   else
-    SCM_WTA_DISPATCH_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part);
+    return scm_wta_dispatch_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part);
 }
 #undef FUNC_NAME
 
@@ -8942,7 +9158,7 @@ SCM_PRIMITIVE_GENERIC (scm_imag_part, "imag-part", 1, 0, 0,
   else if (SCM_I_INUMP (z) || SCM_REALP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
     return SCM_INUM0;
   else
-    SCM_WTA_DISPATCH_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part);
+    return scm_wta_dispatch_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part);
 }
 #undef FUNC_NAME
 
@@ -8958,7 +9174,7 @@ SCM_PRIMITIVE_GENERIC (scm_numerator, "numerator", 1, 0, 0,
   else if (SCM_REALP (z))
     return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator);
+    return scm_wta_dispatch_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator);
 }
 #undef FUNC_NAME
 
@@ -8975,7 +9191,8 @@ SCM_PRIMITIVE_GENERIC (scm_denominator, "denominator", 1, 0, 0,
   else if (SCM_REALP (z))
     return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
   else
-    SCM_WTA_DISPATCH_1 (g_scm_denominator, z, SCM_ARG1, s_scm_denominator);
+    return scm_wta_dispatch_1 (g_scm_denominator, z, SCM_ARG1,
+                               s_scm_denominator);
 }
 #undef FUNC_NAME
 
@@ -9018,7 +9235,8 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0,
         SCM_FRACTION_DENOMINATOR (z));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude);
+    return scm_wta_dispatch_1 (g_scm_magnitude, z, SCM_ARG1,
+                               s_scm_magnitude);
 }
 #undef FUNC_NAME
 
@@ -9065,7 +9283,7 @@ SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0,
       else return scm_from_double (atan2 (0.0, -1.0));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
+    return scm_wta_dispatch_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
 }
 #undef FUNC_NAME
 
@@ -9084,7 +9302,8 @@ SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact, "exact->inexact", 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"
 }