Rename {euclidean,centered}_quo_rem to {euclidean,centered}_divide
[bpt/guile.git] / libguile / numbers.c
index 6583103..41d178b 100644 (file)
@@ -1,4 +1,4 @@
-/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
  *
  * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
  * and Bellcore.  See scm_divide.
@@ -60,6 +60,7 @@
 #include "libguile/root.h"
 #include "libguile/smob.h"
 #include "libguile/strings.h"
+#include "libguile/bdw-gc.h"
 
 #include "libguile/validate.h"
 #include "libguile/numbers.h"
@@ -67,8 +68,6 @@
 
 #include "libguile/eq.h"
 
-#include "libguile/discouraged.h"
-
 /* values per glibc, if not already defined */
 #ifndef M_LOG10E
 #define M_LOG10E   0.43429448190325182765
 #define M_PI       3.14159265358979323846
 #endif
 
+typedef scm_t_signed_bits scm_t_inum;
+#define scm_from_inum(x) (scm_from_signed_integer (x))
+
+/* Tests to see if a C double is neither infinite nor a NaN.
+   TODO: if it's available, use C99's isfinite(x) instead */
+#define DOUBLE_IS_FINITE(x) (!isinf(x) && !isnan(x))
+
 \f
 
 /*
 /* the macro above will not work as is with fractions */
 
 
+static SCM flo0;
+static SCM exactly_one_half;
+
 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
 
 /* FLOBUFLEN is the maximum number of characters neccessary for the
  */
 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
 
-#if defined (SCO)
-#if ! defined (HAVE_ISNAN)
-#define HAVE_ISNAN
-static int
-isnan (double x)
-{
-  return (IsNANorINF (x) && NaN (x) && ! IsINF (x)) ? 1 : 0;
-}
-#endif
-#if ! defined (HAVE_ISINF)
-#define HAVE_ISINF
-static int
-isinf (double x)
-{
-  return (IsNANorINF (x) && IsINF (x)) ? 1 : 0;
-}
-
-#endif
-#endif
-
 
 #if !defined (HAVE_ASINH)
 static double asinh (double x) { return log (x + sqrt (x * x + 1)); }
@@ -141,35 +130,11 @@ static double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); }
    mpz_cmp_d is supposed to do this itself.  */
 #if 1
 #define xmpz_cmp_d(z, d)                                \
-  (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
+  (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
 #else
 #define xmpz_cmp_d(z, d)  mpz_cmp_d (z, d)
 #endif
 
-/* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
-   isinf.  It does have finite and isnan though, hence the use of those.
-   fpclass would be a possibility on that system too.  */
-static int
-xisinf (double x)
-{
-#if defined (HAVE_ISINF)
-  return isinf (x);
-#elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
-  return (! (finite (x) || isnan (x)));
-#else
-  return 0;
-#endif
-}
-
-static int
-xisnan (double x)
-{
-#if defined (HAVE_ISNAN)
-  return isnan (x);
-#else
-  return 0;
-#endif
-}
 
 #if defined (GUILE_I)
 #if HAVE_COMPLEX_DOUBLE
@@ -196,21 +161,66 @@ scm_from_complex_double (complex double z)
 static mpz_t z_negative_one;
 
 \f
+/* Clear the `mpz_t' embedded in bignum PTR.  */
+static void
+finalize_bignum (GC_PTR ptr, GC_PTR data)
+{
+  SCM bignum;
+
+  bignum = PTR2SCM (ptr);
+  mpz_clear (SCM_I_BIG_MPZ (bignum));
+}
+
+/* Return a new uninitialized bignum.  */
+static inline SCM
+make_bignum (void)
+{
+  scm_t_bits *p;
+  GC_finalization_proc prev_finalizer;
+  GC_PTR prev_finalizer_data;
+
+  /* Allocate one word for the type tag and enough room for an `mpz_t'.  */
+  p = scm_gc_malloc_pointerless (sizeof (scm_t_bits) + sizeof (mpz_t),
+                                "bignum");
+  p[0] = scm_tc16_big;
+
+  GC_REGISTER_FINALIZER_NO_ORDER (p, finalize_bignum, NULL,
+                                 &prev_finalizer,
+                                 &prev_finalizer_data);
+
+  return SCM_PACK (p);
+}
+
 
 SCM
 scm_i_mkbig ()
 {
   /* Return a newly created bignum. */
-  SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+  SCM z = make_bignum ();
   mpz_init (SCM_I_BIG_MPZ (z));
   return z;
 }
 
+static SCM
+scm_i_inum2big (scm_t_inum x)
+{
+  /* Return a newly created bignum initialized to X. */
+  SCM z = make_bignum ();
+#if SIZEOF_VOID_P == SIZEOF_LONG
+  mpz_init_set_si (SCM_I_BIG_MPZ (z), x);
+#else
+  /* Note that in this case, you'll also have to check all mpz_*_ui and
+     mpz_*_si invocations in Guile. */
+#error creation of mpz not implemented for this inum size
+#endif
+  return z;
+}
+
 SCM
 scm_i_long2big (long x)
 {
   /* Return a newly created bignum initialized to X. */
-  SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+  SCM z = make_bignum ();
   mpz_init_set_si (SCM_I_BIG_MPZ (z), x);
   return z;
 }
@@ -219,7 +229,7 @@ SCM
 scm_i_ulong2big (unsigned long x)
 {
   /* Return a newly created bignum initialized to X. */
-  SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+  SCM z = make_bignum ();
   mpz_init_set_ui (SCM_I_BIG_MPZ (z), x);
   return z;
 }
@@ -228,7 +238,7 @@ SCM
 scm_i_clonebig (SCM src_big, int same_sign_p)
 {
   /* Copy src_big's value, negate it if same_sign_p is false, and return. */
-  SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+  SCM z = make_bignum ();
   mpz_init_set (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (src_big));
   if (!same_sign_p)
     mpz_neg (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (z));
@@ -249,7 +259,7 @@ SCM
 scm_i_dbl2big (double d)
 {
   /* results are only defined if d is an integer */
-  SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+  SCM z = make_bignum ();
   mpz_init_set_d (SCM_I_BIG_MPZ (z), d);
   return z;
 }
@@ -275,7 +285,7 @@ scm_i_dbl2num (double u)
 
   if (u < (double) (SCM_MOST_POSITIVE_FIXNUM+1)
       && u >= (double) SCM_MOST_NEGATIVE_FIXNUM)
-    return SCM_I_MAKINUM ((long) u);
+    return SCM_I_MAKINUM ((scm_t_inum) u);
   else
     return scm_i_dbl2big (u);
 }
@@ -360,7 +370,7 @@ scm_i_normbig (SCM b)
   /* presume b is a bignum */
   if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b)))
     {
-      long val = mpz_get_si (SCM_I_BIG_MPZ (b));
+      scm_t_inum val = mpz_get_si (SCM_I_BIG_MPZ (b));
       if (SCM_FIXABLE (val))
         b = SCM_I_MAKINUM (val);
     }
@@ -373,13 +383,13 @@ scm_i_mpz2num (mpz_t b)
   /* convert a mpz number to a SCM number. */
   if (mpz_fits_slong_p (b))
     {
-      long val = mpz_get_si (b);
+      scm_t_inum val = mpz_get_si (b);
       if (SCM_FIXABLE (val))
         return SCM_I_MAKINUM (val);
     }
 
   {
-    SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
+    SCM z = make_bignum ();
     mpz_init_set (SCM_I_BIG_MPZ (z), b);
     return z;
   }
@@ -398,7 +408,7 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
     {
       if (scm_is_eq (denominator, SCM_INUM0))
        scm_num_overflow ("make-ratio");
-      if (scm_is_eq (denominator, SCM_I_MAKINUM(1)))
+      if (scm_is_eq (denominator, SCM_INUM1))
        return numerator;
     }
   else 
@@ -422,15 +432,15 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
   */
   if (SCM_I_INUMP (numerator))
     {
-      long  x = SCM_I_INUM (numerator);
+      scm_t_inum x = SCM_I_INUM (numerator);
       if (scm_is_eq (numerator, SCM_INUM0))
        return SCM_INUM0;
       if (SCM_I_INUMP (denominator))
        {
-         long y;
+         scm_t_inum y;
          y = SCM_I_INUM (denominator);
          if (x == y)
-           return SCM_I_MAKINUM(1);
+           return SCM_INUM1;
          if ((x % y) == 0)
            return SCM_I_MAKINUM (x / y);
        }
@@ -450,14 +460,14 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
     {
       if (SCM_I_INUMP (denominator))
        {
-         long yy = SCM_I_INUM (denominator);
+         scm_t_inum yy = SCM_I_INUM (denominator);
          if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
            return scm_divide (numerator, denominator);
        }
       else
        {
          if (scm_is_eq (numerator, denominator))
-           return SCM_I_MAKINUM(1);
+           return SCM_INUM1;
          if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator),
                               SCM_I_BIG_MPZ (denominator)))
            return scm_divide(numerator, denominator);
@@ -468,7 +478,7 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
    */
   {
     SCM divisor = scm_gcd (numerator, denominator);
-    if (!(scm_is_eq (divisor, SCM_I_MAKINUM(1))))
+    if (!(scm_is_eq (divisor, SCM_INUM1)))
       {
        numerator = scm_divide (numerator, divisor);
        denominator = scm_divide (denominator, divisor);
@@ -488,26 +498,39 @@ scm_i_fraction2double (SCM z)
                                         SCM_FRACTION_DENOMINATOR (z)));
 }
 
-SCM_DEFINE (scm_exact_p, "exact?", 1, 0, 0, 
-            (SCM x),
+SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0, 
+                      (SCM x),
            "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
            "otherwise.")
 #define FUNC_NAME s_scm_exact_p
 {
-  if (SCM_I_INUMP (x))
-    return SCM_BOOL_T;
-  if (SCM_BIGP (x))
+  if (SCM_INEXACTP (x))
+    return SCM_BOOL_F;
+  else if (SCM_NUMBERP (x))
     return SCM_BOOL_T;
-  if (SCM_FRACTIONP (x))
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_exact_p, x, 1, s_scm_exact_p);
+}
+#undef FUNC_NAME
+
+
+SCM_PRIMITIVE_GENERIC (scm_inexact_p, "inexact?", 1, 0, 0,
+            (SCM x),
+           "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
+           "else.")
+#define FUNC_NAME s_scm_inexact_p
+{
+  if (SCM_INEXACTP (x))
     return SCM_BOOL_T;
-  if (SCM_NUMBERP (x))
+  else if (SCM_NUMBERP (x))
     return SCM_BOOL_F;
-  SCM_WRONG_TYPE_ARG (1, x);
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p);
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0, 
+SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0, 
             (SCM n),
            "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
            "otherwise.")
@@ -515,7 +538,7 @@ SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
 {
   if (SCM_I_INUMP (n))
     {
-      long val = SCM_I_INUM (n);
+      scm_t_inum val = SCM_I_INUM (n);
       return scm_from_bool ((val & 1L) != 0);
     }
   else if (SCM_BIGP (n))
@@ -524,25 +547,24 @@ SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
       scm_remember_upto_here_1 (n);
       return scm_from_bool (odd_p);
     }
-  else if (scm_is_true (scm_inf_p (n)))
-    return SCM_BOOL_T;
   else if (SCM_REALP (n))
     {
-      double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
-      if (rem == 1.0)
-       return SCM_BOOL_T;
-      else if (rem == 0.0)
-       return SCM_BOOL_F;
-      else
-       SCM_WRONG_TYPE_ARG (1, n);
+      double val = SCM_REAL_VALUE (n);
+      if (DOUBLE_IS_FINITE (val))
+       {
+         double rem = fabs (fmod (val, 2.0));
+         if (rem == 1.0)
+           return SCM_BOOL_T;
+         else if (rem == 0.0)
+           return SCM_BOOL_F;
+       }
     }
-  else
-    SCM_WRONG_TYPE_ARG (1, n);
+  SCM_WTA_DISPATCH_1 (g_scm_odd_p, n, 1, s_scm_odd_p);
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_even_p, "even?", 1, 0, 0, 
+SCM_PRIMITIVE_GENERIC (scm_even_p, "even?", 1, 0, 0, 
             (SCM n),
            "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
            "otherwise.")
@@ -550,7 +572,7 @@ SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
 {
   if (SCM_I_INUMP (n))
     {
-      long val = SCM_I_INUM (n);
+      scm_t_inum val = SCM_I_INUM (n);
       return scm_from_bool ((val & 1L) == 0);
     }
   else if (SCM_BIGP (n))
@@ -559,52 +581,64 @@ SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
       scm_remember_upto_here_1 (n);
       return scm_from_bool (even_p);
     }
-  else if (scm_is_true (scm_inf_p (n)))
-    return SCM_BOOL_T;
   else if (SCM_REALP (n))
     {
-      double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
-      if (rem == 1.0)
-       return SCM_BOOL_F;
-      else if (rem == 0.0)
-       return SCM_BOOL_T;
-      else
-       SCM_WRONG_TYPE_ARG (1, n);
+      double val = SCM_REAL_VALUE (n);
+      if (DOUBLE_IS_FINITE (val))
+       {
+         double rem = fabs (fmod (val, 2.0));
+         if (rem == 1.0)
+           return SCM_BOOL_F;
+         else if (rem == 0.0)
+           return SCM_BOOL_T;
+       }
     }
+  SCM_WTA_DISPATCH_1 (g_scm_even_p, n, 1, s_scm_even_p);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_finite_p, "finite?", 1, 0, 0,
+                      (SCM x),
+           "Return @code{#t} if the real number @var{x} is neither\n"
+           "infinite nor a NaN, @code{#f} otherwise.")
+#define FUNC_NAME s_scm_finite_p
+{
+  if (SCM_REALP (x))
+    return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x)));
+  else if (scm_is_real (x))
+    return SCM_BOOL_T;
   else
-    SCM_WRONG_TYPE_ARG (1, n);
+    SCM_WTA_DISPATCH_1 (g_scm_finite_p, x, 1, s_scm_finite_p);
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0, 
-            (SCM x),
-           "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
-           "or @samp{-inf.0}, @code{#f} otherwise.")
+SCM_PRIMITIVE_GENERIC (scm_inf_p, "inf?", 1, 0, 0, 
+                      (SCM x),
+       "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
+        "@samp{-inf.0}.  Otherwise return @code{#f}.")
 #define FUNC_NAME s_scm_inf_p
 {
   if (SCM_REALP (x))
-    return scm_from_bool (xisinf (SCM_REAL_VALUE (x)));
-  else if (SCM_COMPLEXP (x))
-    return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x))
-                         || xisinf (SCM_COMPLEX_IMAG (x)));
-  else
+    return scm_from_bool (isinf (SCM_REAL_VALUE (x)));
+  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);
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0, 
-            (SCM n),
-           "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
-           "otherwise.")
+SCM_PRIMITIVE_GENERIC (scm_nan_p, "nan?", 1, 0, 0, 
+                      (SCM x),
+           "Return @code{#t} if the real number @var{x} is a NaN,\n"
+            "or @code{#f} otherwise.")
 #define FUNC_NAME s_scm_nan_p
 {
-  if (SCM_REALP (n))
-    return scm_from_bool (xisnan (SCM_REAL_VALUE (n)));
-  else if (SCM_COMPLEXP (n))
-    return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n))
-                    || xisnan (SCM_COMPLEX_IMAG (n)));
-  else
+  if (SCM_REALP (x))
+    return scm_from_bool (isnan (SCM_REAL_VALUE (x)));
+  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);
 }
 #undef FUNC_NAME
 
@@ -617,8 +651,6 @@ static double guile_NaN;
 static void
 guile_ieee_init (void)
 {
-#if defined (HAVE_ISINF) || defined (HAVE_FINITE)
-
 /* Some version of gcc on some old version of Linux used to crash when
    trying to make Inf and NaN.  */
 
@@ -645,10 +677,6 @@ guile_ieee_init (void)
     }
 #endif
 
-#endif
-
-#if defined (HAVE_ISNAN)
-
 #ifdef NAN
   /* C99 NAN, when available */
   guile_NaN = NAN;
@@ -661,8 +689,6 @@ guile_ieee_init (void)
 #else
   guile_NaN = guile_Inf / guile_Inf;
 #endif
-
-#endif
 }
 
 SCM_DEFINE (scm_inf, "inf", 0, 0, 0, 
@@ -699,17 +725,17 @@ SCM_DEFINE (scm_nan, "nan", 0, 0, 0,
 SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
                       (SCM x),
                       "Return the absolute value of @var{x}.")
-#define FUNC_NAME
+#define FUNC_NAME s_scm_abs
 {
   if (SCM_I_INUMP (x))
     {
-      long int xx = SCM_I_INUM (x);
+      scm_t_inum xx = SCM_I_INUM (x);
       if (xx >= 0)
        return x;
       else if (SCM_POSFIXABLE (-xx))
        return SCM_I_MAKINUM (-xx);
       else
-       return scm_i_long2big (-xx);
+       return scm_i_inum2big (-xx);
     }
   else if (SCM_BIGP (x))
     {
@@ -741,27 +767,26 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
 #undef FUNC_NAME
 
 
-SCM_GPROC (s_quotient, "quotient", 2, 0, 0, scm_quotient, g_quotient);
-/* "Return the quotient of the numbers @var{x} and @var{y}."
- */
-SCM
-scm_quotient (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_quotient, "quotient", 2, 0, 0,
+                      (SCM x, SCM y),
+       "Return the quotient of the numbers @var{x} and @var{y}.")
+#define FUNC_NAME s_scm_quotient
 {
-  if (SCM_I_INUMP (x))
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      long xx = SCM_I_INUM (x);
-      if (SCM_I_INUMP (y))
+      scm_t_inum xx = SCM_I_INUM (x);
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_quotient);
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_quotient);
          else
            {
-             long z = xx / yy;
-             if (SCM_FIXABLE (z))
+             scm_t_inum z = xx / yy;
+             if (SCM_LIKELY (SCM_FIXABLE (z)))
                return SCM_I_MAKINUM (z);
              else
-               return scm_i_long2big (z);
+               return scm_i_inum2big (z);
            }
        }
       else if (SCM_BIGP (y))
@@ -775,19 +800,19 @@ scm_quotient (SCM x, SCM y)
               return SCM_I_MAKINUM (-1);
             }
          else
-           return SCM_I_MAKINUM (0);
+           return SCM_INUM0;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
+       SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
     }
   else if (SCM_BIGP (x))
     {
-      if (SCM_I_INUMP (y))
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_quotient);
-         else if (yy == 1)
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_quotient);
+         else if (SCM_UNLIKELY (yy == 1))
            return x;
          else
            {
@@ -815,32 +840,35 @@ scm_quotient (SCM x, SCM y)
          return scm_i_normbig (result);
        }
       else
-       SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
+       SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG1, s_quotient);
+    SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient);
 }
+#undef FUNC_NAME
 
-SCM_GPROC (s_remainder, "remainder", 2, 0, 0, scm_remainder, g_remainder);
-/* "Return the remainder of the numbers @var{x} and @var{y}.\n"
- * "@lisp\n"
- * "(remainder 13 4) @result{} 1\n"
- * "(remainder -13 4) @result{} -1\n"
- * "@end lisp"
- */
-SCM
-scm_remainder (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_remainder, "remainder", 2, 0, 0,
+                      (SCM x, SCM y),
+       "Return the remainder of the numbers @var{x} and @var{y}.\n"
+       "@lisp\n"
+       "(remainder 13 4) @result{} 1\n"
+       "(remainder -13 4) @result{} -1\n"
+       "@end lisp")
+#define FUNC_NAME s_scm_remainder
 {
-  if (SCM_I_INUMP (x))
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      if (SCM_I_INUMP (y))
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_remainder);
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_remainder);
          else
            {
-             long z = SCM_I_INUM (x) % yy;
+             /* C99 specifies that "%" is the remainder corresponding to a
+                 quotient rounded towards zero, and that's also traditional
+                 for machine division, so z here should be well defined.  */
+             scm_t_inum z = SCM_I_INUM (x) % yy;
              return SCM_I_MAKINUM (z);
            }
        }
@@ -852,21 +880,21 @@ scm_remainder (SCM x, SCM y)
             {
               /* Special case:  x == fixnum-min && y == abs (fixnum-min) */
              scm_remember_upto_here_1 (y);
-              return SCM_I_MAKINUM (0);
+              return SCM_INUM0;
             }
          else
            return x;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
+       SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
     }
   else if (SCM_BIGP (x))
     {
-      if (SCM_I_INUMP (y))
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_remainder);
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_remainder);
          else
            {
              SCM result = scm_i_mkbig ();
@@ -887,38 +915,38 @@ scm_remainder (SCM x, SCM y)
          return scm_i_normbig (result);
        }
       else
-       SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
+       SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder);
     }
   else
-    SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG1, s_remainder);
+    SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_modulo, "modulo", 2, 0, 0, scm_modulo, g_modulo);
-/* "Return the modulo of the numbers @var{x} and @var{y}.\n"
- * "@lisp\n"
- * "(modulo 13 4) @result{} 1\n"
- * "(modulo -13 4) @result{} 3\n"
- * "@end lisp"
- */
-SCM
-scm_modulo (SCM x, SCM y)
+SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
+                      (SCM x, SCM y),
+       "Return the modulo of the numbers @var{x} and @var{y}.\n"
+       "@lisp\n"
+       "(modulo 13 4) @result{} 1\n"
+       "(modulo -13 4) @result{} 3\n"
+       "@end lisp")
+#define FUNC_NAME s_scm_modulo
 {
-  if (SCM_I_INUMP (x))
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      long xx = SCM_I_INUM (x);
-      if (SCM_I_INUMP (y))
+      scm_t_inum xx = SCM_I_INUM (x);
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_modulo);
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_modulo);
          else
            {
              /* C99 specifies that "%" is the remainder corresponding to a
                  quotient rounded towards zero, and that's also traditional
                  for machine division, so z here should be well defined.  */
-             long z = xx % yy;
-             long result;
+             scm_t_inum z = xx % yy;
+             scm_t_inum result;
 
              if (yy < 0)
                {
@@ -977,15 +1005,15 @@ scm_modulo (SCM x, SCM y)
            }
        }
       else
-       SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
+       SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
     }
   else if (SCM_BIGP (x))
     {
-      if (SCM_I_INUMP (y))
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         if (yy == 0)
-           scm_num_overflow (s_modulo);
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_modulo);
          else
            {
              SCM result = scm_i_mkbig ();
@@ -1000,32 +1028,1255 @@ scm_modulo (SCM x, SCM y)
              return scm_i_normbig (result);
            }
        }
-      else if (SCM_BIGP (y))
+      else if (SCM_BIGP (y))
+       {
+         SCM result = scm_i_mkbig ();
+         int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
+         SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
+         mpz_mod (SCM_I_BIG_MPZ (result),
+                  SCM_I_BIG_MPZ (x),
+                  SCM_I_BIG_MPZ (pos_y));
+        
+         scm_remember_upto_here_1 (x);
+         if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
+           mpz_add (SCM_I_BIG_MPZ (result),
+                    SCM_I_BIG_MPZ (y),
+                    SCM_I_BIG_MPZ (result));
+         scm_remember_upto_here_2 (y, pos_y);
+         return scm_i_normbig (result);
+       }
+      else
+       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);
+}
+#undef FUNC_NAME
+
+static SCM scm_i_inexact_euclidean_quotient (double x, double y);
+static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} such that\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "where @math{0 <= @var{r} < abs(@var{y})}.\n"
+                      "@lisp\n"
+                      "(euclidean-quotient 123 10) @result{} 12\n"
+                      "(euclidean-quotient 123 -10) @result{} -12\n"
+                      "(euclidean-quotient -123 10) @result{} -13\n"
+                      "(euclidean-quotient -123 -10) @result{} 13\n"
+                      "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
+                      "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_quotient
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_quotient);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             if (xx < qq * yy)
+               {
+                 if (yy > 0)
+                   qq--;
+                 else
+                   qq++;
+               }
+             return SCM_I_MAKINUM (qq);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         if (SCM_I_INUM (x) >= 0)
+           return SCM_INUM0;
+         else
+           return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y)));
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_quotient);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             if (yy > 0)
+               mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy);
+             else
+               {
+                 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+               }
+             scm_remember_upto_here_1 (x);
+             return scm_i_normbig (q);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM q = scm_i_mkbig ();
+         if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           mpz_fdiv_q (SCM_I_BIG_MPZ (q),
+                       SCM_I_BIG_MPZ (x),
+                       SCM_I_BIG_MPZ (y));
+         else
+           mpz_cdiv_q (SCM_I_BIG_MPZ (q),
+                       SCM_I_BIG_MPZ (x),
+                       SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_i_normbig (q);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_quotient
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                           s_scm_euclidean_quotient);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_quotient
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_quotient (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
+                       s_scm_euclidean_quotient);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_quotient (double x, double y)
+{
+  if (SCM_LIKELY (y > 0))
+    return scm_from_double (floor (x / y));
+  else if (SCM_LIKELY (y < 0))
+    return scm_from_double (ceil (x / y));
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_quotient);  /* or return a NaN? */
+  else
+    return scm_nan ();
+}
+
+/* Compute exact euclidean_quotient the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_quotient (SCM x, SCM y)
+{
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1,
+                       s_scm_euclidean_quotient);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2,
+                       s_scm_euclidean_quotient);
+  else if (scm_is_true (scm_positive_p (y)))
+    return scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    return scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_quotient);
+}
+
+static SCM scm_i_inexact_euclidean_remainder (double x, double y);
+static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the real number @var{r} such that\n"
+                      "@math{0 <= @var{r} < abs(@var{y})} and\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "for some integer @var{q}.\n"
+                      "@lisp\n"
+                      "(euclidean-remainder 123 10) @result{} 3\n"
+                      "(euclidean-remainder 123 -10) @result{} 3\n"
+                      "(euclidean-remainder -123 10) @result{} 7\n"
+                      "(euclidean-remainder -123 -10) @result{} 7\n"
+                      "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
+                      "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_remainder
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_remainder);
+         else
+           {
+             scm_t_inum rr = SCM_I_INUM (x) % yy;
+             if (rr >= 0)
+               return SCM_I_MAKINUM (rr);
+             else if (yy > 0)
+               return SCM_I_MAKINUM (rr + yy);
+             else
+               return SCM_I_MAKINUM (rr - yy);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         scm_t_inum xx = SCM_I_INUM (x);
+         if (xx >= 0)
+           return x;
+         else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             return scm_i_normbig (r);
+           }
+         else
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r));
+             return scm_i_normbig (r);
+           }
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_remainder);
+         else
+           {
+             scm_t_inum rr;
+             if (yy < 0)
+               yy = -yy;
+             rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy);
+             scm_remember_upto_here_1 (x);
+             return SCM_I_MAKINUM (rr);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM r = scm_i_mkbig ();
+         mpz_mod (SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x),
+                  SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_i_normbig (r);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_remainder
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                           s_scm_euclidean_remainder);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_remainder
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_remainder (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
+                       s_scm_euclidean_remainder);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_remainder (double x, double y)
+{
+  double q;
+
+  /* Although it would be more efficient to use fmod here, we can't
+     because it would in some cases produce results inconsistent with
+     scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
+     even close).  In particular, when x is very close to a multiple of
+     y, then r might be either 0.0 or abs(y)-epsilon, but those two
+     cases must correspond to different choices of q.  If r = 0.0 then q
+     must be x/y, and if r = abs(y) then q must be (x-r)/y.  If quotient
+     chooses one and remainder chooses the other, it would be bad.  This
+     problem was observed with x = 130.0 and y = 10/7. */
+  if (SCM_LIKELY (y > 0))
+    q = floor (x / y);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x / y);
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_remainder);  /* or return a NaN? */
+  else
+    return scm_nan ();
+  return scm_from_double (x - q * y);
+}
+
+/* Compute exact euclidean_remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_remainder (SCM x, SCM y)
+{
+  SCM q;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1,
+                       s_scm_euclidean_remainder);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2,
+                       s_scm_euclidean_remainder);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_remainder);
+  return scm_difference (x, scm_product (y, q));
+}
+
+
+static SCM scm_i_inexact_euclidean_divide (double x, double y);
+static SCM scm_i_slow_exact_euclidean_divide (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} and the real number @var{r}\n"
+                      "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "and @math{0 <= @var{r} < abs(@var{y})}.\n"
+                      "@lisp\n"
+                      "(euclidean/ 123 10) @result{} 12 and 3\n"
+                      "(euclidean/ 123 -10) @result{} -12 and 3\n"
+                      "(euclidean/ -123 10) @result{} -13 and 7\n"
+                      "(euclidean/ -123 -10) @result{} 13 and 7\n"
+                      "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+                      "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_euclidean_divide
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_divide);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * yy;
+             if (rr < 0)
+               {
+                 if (yy > 0)
+                   { rr += yy; qq--; }
+                 else
+                   { rr -= yy; qq++; }
+               }
+             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         scm_t_inum xx = SCM_I_INUM (x);
+         if (xx >= 0)
+           return scm_values (scm_list_2 (SCM_INUM0, x));
+         else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             return scm_values
+               (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r)));
+           }
+         else
+           {
+             SCM r = scm_i_mkbig ();
+             mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx);
+             scm_remember_upto_here_1 (y);
+             mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r));
+             return scm_values (scm_list_2 (SCM_INUM1, scm_i_normbig (r)));
+           }
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_divide
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_divide (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+                           s_scm_euclidean_divide);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_euclidean_divide);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             SCM r = scm_i_mkbig ();
+             if (yy > 0)
+               mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                               SCM_I_BIG_MPZ (x), yy);
+             else
+               {
+                 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                                 SCM_I_BIG_MPZ (x), -yy);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+               }
+             scm_remember_upto_here_1 (x);
+             return scm_values (scm_list_2 (scm_i_normbig (q),
+                                            scm_i_normbig (r)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         SCM q = scm_i_mkbig ();
+         SCM r = scm_i_mkbig ();
+         if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+           mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                        SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+         else
+           mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                        SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+         scm_remember_upto_here_2 (x, y);
+         return scm_values (scm_list_2 (scm_i_normbig (q),
+                                        scm_i_normbig (r)));
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_divide
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_euclidean_divide (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+                           s_scm_euclidean_divide);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_euclidean_divide
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+                           s_scm_euclidean_divide);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_euclidean_divide
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_euclidean_divide (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
+                       s_scm_euclidean_divide);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_euclidean_divide (double x, double y)
+{
+  double q, r;
+
+  if (SCM_LIKELY (y > 0))
+    q = floor (x / y);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x / y);
+  else if (y == 0)
+    scm_num_overflow (s_scm_euclidean_divide);  /* or return a NaN? */
+  else
+    q = guile_NaN;
+  r = x - q * y;
+  return scm_values (scm_list_2 (scm_from_double (q),
+                                scm_from_double (r)));
+}
+
+/* Compute exact euclidean quotient and remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_euclidean_divide (SCM x, SCM y)
+{
+  SCM q, r;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1,
+                       s_scm_euclidean_divide);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2,
+                       s_scm_euclidean_divide);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_divide (x, y));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_divide (x, y));
+  else
+    scm_num_overflow (s_scm_euclidean_divide);
+  r = scm_difference (x, scm_product (q, y));
+  return scm_values (scm_list_2 (q, r));
+}
+
+static SCM scm_i_inexact_centered_quotient (double x, double y);
+static SCM scm_i_bigint_centered_quotient (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_quotient (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} such that\n"
+                      "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
+                      "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+                      "@lisp\n"
+                      "(centered-quotient 123 10) @result{} 12\n"
+                      "(centered-quotient 123 -10) @result{} -12\n"
+                      "(centered-quotient -123 10) @result{} -12\n"
+                      "(centered-quotient -123 -10) @result{} 12\n"
+                      "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
+                      "(centered-quotient 16/3 -10/7) @result{} -4\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_quotient
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_quotient);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * yy;
+             if (SCM_LIKELY (xx > 0))
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr >= (yy + 1) / 2)
+                       qq++;
+                   }
+                 else
+                   {
+                     if (rr >= (1 - yy) / 2)
+                       qq--;
+                   }
+               }
+             else
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr < -yy / 2)
+                       qq--;
+                   }
+                 else
+                   {
+                     if (rr < yy / 2)
+                       qq++;
+                   }
+               }
+             return SCM_I_MAKINUM (qq);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         /* Pass a denormalized bignum version of x (even though it
+            can fit in a fixnum) to scm_i_bigint_centered_quotient */
+         return scm_i_bigint_centered_quotient
+           (scm_i_long2big (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                           s_scm_centered_quotient);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_quotient);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             scm_t_inum rr;
+             /* Arrange for rr to initially be non-positive,
+                because that simplifies the test to see
+                if it is within the needed bounds. */
+             if (yy > 0)
+               {
+                 rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                       SCM_I_BIG_MPZ (x), yy);
+                 scm_remember_upto_here_1 (x);
+                 if (rr < -yy / 2)
+                   mpz_sub_ui (SCM_I_BIG_MPZ (q),
+                               SCM_I_BIG_MPZ (q), 1);
+               }
+             else
+               {
+                 rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                       SCM_I_BIG_MPZ (x), -yy);
+                 scm_remember_upto_here_1 (x);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+                 if (rr < yy / 2)
+                   mpz_add_ui (SCM_I_BIG_MPZ (q),
+                               SCM_I_BIG_MPZ (q), 1);
+               }
+             return scm_i_normbig (q);
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_quotient (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_quotient (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                           s_scm_centered_quotient);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       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);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_quotient
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_quotient (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
+                       s_scm_centered_quotient);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_quotient (double x, double y)
+{
+  if (SCM_LIKELY (y > 0))
+    return scm_from_double (floor (x/y + 0.5));
+  else if (SCM_LIKELY (y < 0))
+    return scm_from_double (ceil (x/y - 0.5));
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_quotient);  /* or return a NaN? */
+  else
+    return scm_nan ();
+}
+
+/* Assumes that both x and y are bigints, though
+   x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_quotient (SCM x, SCM y)
+{
+  SCM q, r, min_r;
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+  q = scm_i_mkbig ();
+  r = scm_i_mkbig ();
+
+  /* min_r will eventually become -abs(y)/2 */
+  min_r = scm_i_mkbig ();
+  mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+                  SCM_I_BIG_MPZ (y), 1);
+
+  /* Arrange for rr to initially be non-positive,
+     because that simplifies the test to see
+     if it is within the needed bounds. */
+  if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+    {
+      mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      scm_remember_upto_here_2 (x, y);
+      mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+       mpz_sub_ui (SCM_I_BIG_MPZ (q),
+                   SCM_I_BIG_MPZ (q), 1);
+    }
+  else
+    {
+      mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      scm_remember_upto_here_2 (x, y);
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+       mpz_add_ui (SCM_I_BIG_MPZ (q),
+                   SCM_I_BIG_MPZ (q), 1);
+    }
+  scm_remember_upto_here_2 (r, min_r);
+  return scm_i_normbig (q);
+}
+
+/* Compute exact centered quotient the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_quotient (SCM x, SCM y)
+{
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1,
+                       s_scm_centered_quotient);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2,
+                       s_scm_centered_quotient);
+  else if (scm_is_true (scm_positive_p (y)))
+    return scm_floor (scm_sum (scm_divide (x, y),
+                              exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    return scm_ceiling (scm_difference (scm_divide (x, y),
+                                       exactly_one_half));
+  else
+    scm_num_overflow (s_scm_centered_quotient);
+}
+
+static SCM scm_i_inexact_centered_remainder (double x, double y);
+static SCM scm_i_bigint_centered_remainder (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_remainder (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the real number @var{r} such that\n"
+                      "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
+                      "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "for some integer @var{q}.\n"
+                      "@lisp\n"
+                      "(centered-remainder 123 10) @result{} 3\n"
+                      "(centered-remainder 123 -10) @result{} 3\n"
+                      "(centered-remainder -123 10) @result{} -3\n"
+                      "(centered-remainder -123 -10) @result{} -3\n"
+                      "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
+                      "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_remainder
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_remainder);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum rr = xx % yy;
+             if (SCM_LIKELY (xx > 0))
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr >= (yy + 1) / 2)
+                       rr -= yy;
+                   }
+                 else
+                   {
+                     if (rr >= (1 - yy) / 2)
+                       rr += yy;
+                   }
+               }
+             else
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr < -yy / 2)
+                       rr += yy;
+                   }
+                 else
+                   {
+                     if (rr < yy / 2)
+                       rr -= yy;
+                   }
+               }
+             return SCM_I_MAKINUM (rr);
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         /* Pass a denormalized bignum version of x (even though it
+            can fit in a fixnum) to scm_i_bigint_centered_remainder */
+         return scm_i_bigint_centered_remainder
+           (scm_i_long2big (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                           s_scm_centered_remainder);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_remainder);
+         else
+           {
+             scm_t_inum rr;
+             /* Arrange for rr to initially be non-positive,
+                because that simplifies the test to see
+                if it is within the needed bounds. */
+             if (yy > 0)
+               {
+                 rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy);
+                 scm_remember_upto_here_1 (x);
+                 if (rr < -yy / 2)
+                   rr += yy;
+               }
+             else
+               {
+                 rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy);
+                 scm_remember_upto_here_1 (x);
+                 if (rr < yy / 2)
+                   rr -= yy;
+               }
+             return SCM_I_MAKINUM (rr);
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_remainder (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_remainder (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                           s_scm_centered_remainder);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       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);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_remainder
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_remainder (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
+                       s_scm_centered_remainder);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_remainder (double x, double y)
+{
+  double q;
+
+  /* Although it would be more efficient to use fmod here, we can't
+     because it would in some cases produce results inconsistent with
+     scm_i_inexact_centered_quotient, such that x != r + q * y (not even
+     close).  In particular, when x-y/2 is very close to a multiple of
+     y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
+     two cases must correspond to different choices of q.  If quotient
+     chooses one and remainder chooses the other, it would be bad. */
+  if (SCM_LIKELY (y > 0))
+    q = floor (x/y + 0.5);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x/y - 0.5);
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_remainder);  /* or return a NaN? */
+  else
+    return scm_nan ();
+  return scm_from_double (x - q * y);
+}
+
+/* Assumes that both x and y are bigints, though
+   x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_remainder (SCM x, SCM y)
+{
+  SCM r, min_r;
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+  r = scm_i_mkbig ();
+
+  /* min_r will eventually become -abs(y)/2 */
+  min_r = scm_i_mkbig ();
+  mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+                  SCM_I_BIG_MPZ (y), 1);
+
+  /* Arrange for rr to initially be non-positive,
+     because that simplifies the test to see
+     if it is within the needed bounds. */
+  if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+    {
+      mpz_cdiv_r (SCM_I_BIG_MPZ (r),
+                 SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+       mpz_add (SCM_I_BIG_MPZ (r),
+                SCM_I_BIG_MPZ (r),
+                SCM_I_BIG_MPZ (y));
+    }
+  else
+    {
+      mpz_fdiv_r (SCM_I_BIG_MPZ (r),
+                 SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+       mpz_sub (SCM_I_BIG_MPZ (r),
+                SCM_I_BIG_MPZ (r),
+                SCM_I_BIG_MPZ (y));
+    }
+  scm_remember_upto_here_2 (x, y);
+  return scm_i_normbig (r);
+}
+
+/* Compute exact centered_remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_remainder (SCM x, SCM y)
+{
+  SCM q;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1,
+                       s_scm_centered_remainder);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2,
+                       s_scm_centered_remainder);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_sum (scm_divide (x, y), exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_difference (scm_divide (x, y), exactly_one_half));
+  else
+    scm_num_overflow (s_scm_centered_remainder);
+  return scm_difference (x, scm_product (y, q));
+}
+
+
+static SCM scm_i_inexact_centered_divide (double x, double y);
+static SCM scm_i_bigint_centered_divide (SCM x, SCM y);
+static SCM scm_i_slow_exact_centered_divide (SCM x, SCM y);
+
+SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return the integer @var{q} and the real number @var{r}\n"
+                      "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
+                      "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
+                      "@lisp\n"
+                      "(centered/ 123 10) @result{} 12 and 3\n"
+                      "(centered/ 123 -10) @result{} -12 and 3\n"
+                      "(centered/ -123 10) @result{} -12 and -3\n"
+                      "(centered/ -123 -10) @result{} 12 and -3\n"
+                      "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
+                      "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
+                      "@end lisp")
+#define FUNC_NAME s_scm_centered_divide
+{
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_divide);
+         else
+           {
+             scm_t_inum xx = SCM_I_INUM (x);
+             scm_t_inum qq = xx / yy;
+             scm_t_inum rr = xx - qq * yy;
+             if (SCM_LIKELY (xx > 0))
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr >= (yy + 1) / 2)
+                       { qq++; rr -= yy; }
+                   }
+                 else
+                   {
+                     if (rr >= (1 - yy) / 2)
+                       { qq--; rr += yy; }
+                   }
+               }
+             else
+               {
+                 if (SCM_LIKELY (yy > 0))
+                   {
+                     if (rr < -yy / 2)
+                       { qq--; rr += yy; }
+                   }
+                 else
+                   {
+                     if (rr < yy / 2)
+                       { qq++; rr -= yy; }
+                   }
+               }
+             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       {
+         /* Pass a denormalized bignum version of x (even though it
+            can fit in a fixnum) to scm_i_bigint_centered_divide */
+         return scm_i_bigint_centered_divide
+           (scm_i_long2big (SCM_I_INUM (x)), y);
+       }
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_divide
+         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_divide (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+                           s_scm_centered_divide);
+    }
+  else if (SCM_BIGP (x))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
+       {
+         scm_t_inum yy = SCM_I_INUM (y);
+         if (SCM_UNLIKELY (yy == 0))
+           scm_num_overflow (s_scm_centered_divide);
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             scm_t_inum rr;
+             /* Arrange for rr to initially be non-positive,
+                because that simplifies the test to see
+                if it is within the needed bounds. */
+             if (yy > 0)
+               {
+                 rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                       SCM_I_BIG_MPZ (x), yy);
+                 scm_remember_upto_here_1 (x);
+                 if (rr < -yy / 2)
+                   {
+                     mpz_sub_ui (SCM_I_BIG_MPZ (q),
+                                 SCM_I_BIG_MPZ (q), 1);
+                     rr += yy;
+                   }
+               }
+             else
+               {
+                 rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                       SCM_I_BIG_MPZ (x), -yy);
+                 scm_remember_upto_here_1 (x);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+                 if (rr < yy / 2)
+                   {
+                     mpz_add_ui (SCM_I_BIG_MPZ (q),
+                                 SCM_I_BIG_MPZ (q), 1);
+                     rr -= yy;
+                   }
+               }
+             return scm_values (scm_list_2 (scm_i_normbig (q),
+                                            SCM_I_MAKINUM (rr)));
+           }
+       }
+      else if (SCM_BIGP (y))
+       return scm_i_bigint_centered_divide (x, y);
+      else if (SCM_REALP (y))
+       return scm_i_inexact_centered_divide
+         (scm_i_big2dbl (x), SCM_REAL_VALUE (y));
+      else if (SCM_FRACTIONP (y))
+       return scm_i_slow_exact_centered_divide (x, y);
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+                           s_scm_centered_divide);
+    }
+  else if (SCM_REALP (x))
+    {
+      if (SCM_REALP (y) || SCM_I_INUMP (y) ||
+         SCM_BIGP (y) || SCM_FRACTIONP (y))
+       return scm_i_inexact_centered_divide
+         (SCM_REAL_VALUE (x), scm_to_double (y));
+     else
+       SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+                           s_scm_centered_divide);
+    }
+  else if (SCM_FRACTIONP (x))
+    {
+      if (SCM_REALP (y))
+       return scm_i_inexact_centered_divide
+         (scm_i_fraction2double (x), SCM_REAL_VALUE (y));
+      else
+       return scm_i_slow_exact_centered_divide (x, y);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
+                       s_scm_centered_divide);
+}
+#undef FUNC_NAME
+
+static SCM
+scm_i_inexact_centered_divide (double x, double y)
+{
+  double q, r;
+
+  if (SCM_LIKELY (y > 0))
+    q = floor (x/y + 0.5);
+  else if (SCM_LIKELY (y < 0))
+    q = ceil (x/y - 0.5);
+  else if (y == 0)
+    scm_num_overflow (s_scm_centered_divide);  /* or return a NaN? */
+  else
+    q = guile_NaN;
+  r = x - q * y;
+  return scm_values (scm_list_2 (scm_from_double (q),
+                                scm_from_double (r)));
+}
+
+/* Assumes that both x and y are bigints, though
+   x might be able to fit into a fixnum. */
+static SCM
+scm_i_bigint_centered_divide (SCM x, SCM y)
+{
+  SCM q, r, min_r;
+
+  /* Note that x might be small enough to fit into a
+     fixnum, so we must not let it escape into the wild */
+  q = scm_i_mkbig ();
+  r = scm_i_mkbig ();
+
+  /* min_r will eventually become -abs(y/2) */
+  min_r = scm_i_mkbig ();
+  mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r),
+                  SCM_I_BIG_MPZ (y), 1);
+
+  /* Arrange for rr to initially be non-positive,
+     because that simplifies the test to see
+     if it is within the needed bounds. */
+  if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
+    {
+      mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r));
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
+       {
+         mpz_sub_ui (SCM_I_BIG_MPZ (q),
+                     SCM_I_BIG_MPZ (q), 1);
+         mpz_add (SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (y));
+       }
+    }
+  else
+    {
+      mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
+      if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0)
        {
-           {
-             SCM result = scm_i_mkbig ();
-             int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
-             SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
-             mpz_mod (SCM_I_BIG_MPZ (result),
-                      SCM_I_BIG_MPZ (x),
-                      SCM_I_BIG_MPZ (pos_y));
-        
-             scm_remember_upto_here_1 (x);
-             if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
-               mpz_add (SCM_I_BIG_MPZ (result),
-                        SCM_I_BIG_MPZ (y),
-                        SCM_I_BIG_MPZ (result));
-             scm_remember_upto_here_2 (y, pos_y);
-             return scm_i_normbig (result);
-           }
+         mpz_add_ui (SCM_I_BIG_MPZ (q),
+                     SCM_I_BIG_MPZ (q), 1);
+         mpz_sub (SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (r),
+                  SCM_I_BIG_MPZ (y));
        }
-      else
-       SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
     }
+  scm_remember_upto_here_2 (x, y);
+  return scm_values (scm_list_2 (scm_i_normbig (q),
+                                scm_i_normbig (r)));
+}
+
+/* Compute exact centered quotient and remainder the slow way.
+   We use this only if both arguments are exact,
+   and at least one of them is a fraction */
+static SCM
+scm_i_slow_exact_centered_divide (SCM x, SCM y)
+{
+  SCM q, r;
+
+  if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1,
+                       s_scm_centered_divide);
+  else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)))
+    SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+                       s_scm_centered_divide);
+  else if (scm_is_true (scm_positive_p (y)))
+    q = scm_floor (scm_sum (scm_divide (x, y),
+                           exactly_one_half));
+  else if (scm_is_true (scm_negative_p (y)))
+    q = scm_ceiling (scm_difference (scm_divide (x, y),
+                                    exactly_one_half));
   else
-    SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo);
+    scm_num_overflow (s_scm_centered_divide);
+  r = scm_difference (x, scm_product (q, y));
+  return scm_values (scm_list_2 (q, r));
 }
 
+
 SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1,
                        (SCM x, SCM y, SCM rest),
                        "Return the greatest common divisor of all parameter values.\n"
@@ -1054,19 +2305,19 @@ scm_gcd (SCM x, SCM y)
     {
       if (SCM_I_INUMP (y))
         {
-          long xx = SCM_I_INUM (x);
-          long yy = SCM_I_INUM (y);
-          long u = xx < 0 ? -xx : xx;
-          long v = yy < 0 ? -yy : yy;
-          long result;
+          scm_t_inum xx = SCM_I_INUM (x);
+          scm_t_inum yy = SCM_I_INUM (y);
+          scm_t_inum u = xx < 0 ? -xx : xx;
+          scm_t_inum v = yy < 0 ? -yy : yy;
+          scm_t_inum result;
           if (xx == 0)
            result = v;
          else if (yy == 0)
            result = u;
          else
            {
-             long k = 1;
-             long t;
+             scm_t_inum k = 1;
+             scm_t_inum t;
              /* Determine a common factor 2^k */
              while (!(1 & (u | v)))
                {
@@ -1096,7 +2347,7 @@ scm_gcd (SCM x, SCM y)
            }
           return (SCM_POSFIXABLE (result)
                  ? SCM_I_MAKINUM (result)
-                 : scm_i_long2big (result));
+                 : scm_i_inum2big (result));
         }
       else if (SCM_BIGP (y))
         {
@@ -1110,8 +2361,8 @@ scm_gcd (SCM x, SCM y)
     {
       if (SCM_I_INUMP (y))
         {
-          unsigned long result;
-          long yy;
+          scm_t_bits result;
+          scm_t_inum yy;
         big_inum:
           yy = SCM_I_INUM (y);
           if (yy == 0)
@@ -1122,7 +2373,7 @@ scm_gcd (SCM x, SCM y)
           scm_remember_upto_here_1 (x);
           return (SCM_POSFIXABLE (result) 
                  ? SCM_I_MAKINUM (result)
-                 : scm_from_ulong (result));
+                 : scm_from_unsigned_integer (result));
         }
       else if (SCM_BIGP (y))
         {
@@ -1189,7 +2440,7 @@ scm_lcm (SCM n1, SCM n2)
         inumbig:
           {
             SCM result = scm_i_mkbig ();
-            long nn1 = SCM_I_INUM (n1);
+            scm_t_inum nn1 = SCM_I_INUM (n1);
             if (nn1 == 0) return SCM_INUM0;
             if (nn1 < 0) nn1 = - nn1;
             mpz_lcm_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n2), nn1);
@@ -1279,7 +2530,7 @@ SCM_DEFINE (scm_i_logand, "logand", 0, 2, 1,
 SCM scm_logand (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logand
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1298,7 +2549,7 @@ SCM scm_logand (SCM n1, SCM n2)
       nn1 = SCM_I_INUM (n1);
       if (SCM_I_INUMP (n2))
        {
-         long nn2 = SCM_I_INUM (n2);
+         scm_t_inum nn2 = SCM_I_INUM (n2);
          return SCM_I_MAKINUM (nn1 & nn2);
        }
       else if SCM_BIGP (n2)
@@ -1369,7 +2620,7 @@ SCM_DEFINE (scm_i_logior, "logior", 0, 2, 1,
 SCM scm_logior (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logior
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1459,7 +2710,7 @@ SCM_DEFINE (scm_i_logxor, "logxor", 0, 2, 1,
 SCM scm_logxor (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logxor
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1476,7 +2727,7 @@ SCM scm_logxor (SCM n1, SCM n2)
       nn1 = SCM_I_INUM (n1);
       if (SCM_I_INUMP (n2))
        {
-         long nn2 = SCM_I_INUM (n2);
+         scm_t_inum nn2 = SCM_I_INUM (n2);
          return SCM_I_MAKINUM (nn1 ^ nn2);
        }
       else if (SCM_BIGP (n2))
@@ -1534,14 +2785,14 @@ SCM_DEFINE (scm_logtest, "logtest", 2, 0, 0,
            "@end lisp")
 #define FUNC_NAME s_scm_logtest
 {
-  long int nj;
+  scm_t_inum nj;
 
   if (SCM_I_INUMP (j))
     {
       nj = SCM_I_INUM (j);
       if (SCM_I_INUMP (k))
        {
-         long nk = SCM_I_INUM (k);
+         scm_t_inum nk = SCM_I_INUM (k);
          return scm_from_bool (nj & nk);
        }
       else if (SCM_BIGP (k))
@@ -1783,8 +3034,9 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
            "Return @var{n} raised to the power @var{k}.  @var{k} must be an\n"
            "exact integer, @var{n} can be any number.\n"
            "\n"
-           "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
-           "in the usual way.  @math{@var{n}^0} is 1, as usual, and that\n"
+           "Negative @var{k} is supported, and results in\n"
+           "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
+           "@math{@var{n}^0} is 1, as usual, and that\n"
            "includes @math{0^0} is 1.\n"
            "\n"
            "@lisp\n"
@@ -1795,14 +3047,26 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
            "@end lisp")
 #define FUNC_NAME s_scm_integer_expt
 {
-  long i2 = 0;
+  scm_t_inum i2 = 0;
   SCM z_i2 = SCM_BOOL_F;
   int i2_is_big = 0;
   SCM acc = SCM_I_MAKINUM (1L);
 
-  /* 0^0 == 1 according to R5RS */
-  if (scm_is_eq (n, SCM_INUM0) || scm_is_eq (n, acc))
-    return scm_is_false (scm_zero_p(k)) ? n : acc;
+  SCM_VALIDATE_NUMBER (SCM_ARG1, n);
+  if (!SCM_I_INUMP (k) && !SCM_BIGP (k))
+    SCM_WRONG_TYPE_ARG (2, k);
+
+  if (scm_is_true (scm_zero_p (n)))
+    {
+      if (scm_is_true (scm_zero_p (k)))  /* 0^0 == 1 per R5RS */
+       return acc;        /* return exact 1, regardless of n */
+      else if (scm_is_true (scm_positive_p (k)))
+       return n;
+      else  /* return NaN for (0 ^ k) for negative k per R6RS */
+       return scm_nan ();
+    }
+  else if (scm_is_eq (n, acc))
+    return acc;
   else if (scm_is_eq (n, SCM_I_MAKINUM (-1L)))
     return scm_is_false (scm_even_p (k)) ? n : acc;
 
@@ -1890,7 +3154,7 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
 
   if (SCM_I_INUMP (n))
     {
-      long nn = SCM_I_INUM (n);
+      scm_t_inum nn = SCM_I_INUM (n);
 
       if (bits_to_shift > 0)
         {
@@ -1905,7 +3169,7 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
             return n;
 
           if (bits_to_shift < SCM_I_FIXNUM_BIT-1
-              && ((unsigned long)
+              && ((scm_t_bits)
                   (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
                   <= 1))
             {
@@ -1913,7 +3177,7 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
             }
           else
             {
-              SCM result = scm_i_long2big (nn);
+              SCM result = scm_i_inum2big (nn);
               mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
                             bits_to_shift);
               return result;
@@ -1923,7 +3187,7 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
         {
           bits_to_shift = -bits_to_shift;
           if (bits_to_shift >= SCM_LONG_BIT)
-            return (nn >= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
+            return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM(-1));
           else
             return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
         }
@@ -1986,7 +3250,7 @@ SCM_DEFINE (scm_bit_extract, "bit-extract", 3, 0, 0,
 
   if (SCM_I_INUMP (n))
     {
-      long int in = SCM_I_INUM (n);
+      scm_t_inum in = SCM_I_INUM (n);
 
       /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
          SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
@@ -1998,7 +3262,7 @@ SCM_DEFINE (scm_bit_extract, "bit-extract", 3, 0, 0,
           * special case requires us to produce a result that has
           * more bits than can be stored in a fixnum.
           */
-          SCM result = scm_i_long2big (in);
+          SCM result = scm_i_inum2big (in);
           mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
                            bits);
           return result;
@@ -2057,8 +3321,8 @@ SCM_DEFINE (scm_logcount, "logcount", 1, 0, 0,
 {
   if (SCM_I_INUMP (n))
     {
-      unsigned long int c = 0;
-      long int nn = SCM_I_INUM (n);
+      unsigned long c = 0;
+      scm_t_inum nn = SCM_I_INUM (n);
       if (nn < 0)
         nn = -1 - nn;
       while (nn)
@@ -2105,9 +3369,9 @@ SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0,
 {
   if (SCM_I_INUMP (n))
     {
-      unsigned long int c = 0;
+      unsigned long c = 0;
       unsigned int l = 4;
-      long int nn = SCM_I_INUM (n);
+      scm_t_inum nn = SCM_I_INUM (n);
       if (nn < 0)
        nn = -1 - nn;
       while (nn)
@@ -2179,7 +3443,7 @@ void init_fx_radix(double *fx_list, int radix)
 }
 
 /* use this array as a way to generate a single digit */
-static const char*number_chars="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+static const char number_chars[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 
 static size_t
 idbl2str (double f, char *a, int radix)
@@ -2214,7 +3478,7 @@ idbl2str (double f, char *a, int radix)
       goto zero;       /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
     }
 
-  if (xisinf (f))
+  if (isinf (f))
     {
       if (f < 0)
        strcpy (a, "-inf.0");
@@ -2222,7 +3486,7 @@ idbl2str (double f, char *a, int radix)
        strcpy (a, "+inf.0");
       return ch+6;
     }
-  else if (xisnan (f))
+  else if (isnan (f))
     {
       strcpy (a, "+nan.0");
       return ch+6;
@@ -2374,7 +3638,7 @@ icmplx2str (double real, double imag, char *str, int radix)
     {
       /* Don't output a '+' for negative numbers or for Inf and
         NaN.  They will provide their own sign. */
-      if (0 <= imag && !xisinf (imag) && !xisnan (imag))
+      if (0 <= imag && !isinf (imag) && !isnan (imag))
        str[i++] = '+';
       i += idbl2str (imag, &str[i], radix);
       str[i++] = 'i';
@@ -2421,6 +3685,9 @@ scm_iuint2str (scm_t_uintmax num, int rad, char *p)
   size_t i;
   scm_t_uintmax n = num;
 
+  if (rad < 2 || rad > 36)
+    scm_out_of_range ("scm_iuint2str", scm_from_int (rad));
+
   for (n /= rad; n > 0; n /= rad)
     j++;
 
@@ -2431,7 +3698,7 @@ scm_iuint2str (scm_t_uintmax num, int rad, char *p)
       int d = n % rad;
 
       n /= rad;
-      p[i] = d + ((d < 10) ? '0' : 'a' - 10);
+      p[i] = number_chars[d];
     }
   return j;
 }
@@ -2518,7 +3785,7 @@ scm_i_print_fraction (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
 {
   SCM str;
   str = scm_number_to_string (sexp, SCM_UNDEFINED);
-  scm_lfwrite_str (str, port);
+  scm_display (str, port);
   scm_remember_upto_here_1 (str);
   return !0;
 }
@@ -2564,11 +3831,26 @@ enum t_exactness {NO_EXACTNESS, INEXACT, EXACT};
 
 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
 
-/* In non ASCII-style encodings the following macro might not work. */
-#define XDIGIT2UINT(d)                                                  \
-  (uc_is_property_decimal_digit ((int) (unsigned char) d)               \
-   ? (d) - '0'                                                          \
-   : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
+/* Caller is responsible for checking that the return value is in range
+   for the given radix, which should be <= 36. */
+static unsigned int
+char_decimal_value (scm_t_uint32 c)
+{
+  /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
+     that's certainly above any valid decimal, so we take advantage of
+     that to elide some tests. */
+  unsigned int d = (unsigned int) uc_decimal_value (c);
+
+  /* If that failed, try extended hexadecimals, then. Only accept ascii
+     hexadecimals. */
+  if (d >= 10U)
+    {
+      c = uc_tolower (c);
+      if (c >= (scm_t_uint32) 'a')
+        d = c - (scm_t_uint32)'a' + 10U;
+    }
+  return d;
+}
 
 static SCM
 mem2uinteger (SCM mem, unsigned int *p_idx,
@@ -2587,9 +3869,7 @@ mem2uinteger (SCM mem, unsigned int *p_idx,
     return SCM_BOOL_F;
 
   c = scm_i_string_ref (mem, idx);
-  if (!uc_is_property_ascii_hex_digit ((scm_t_uint32) c))
-    return SCM_BOOL_F;
-  digit_value = XDIGIT2UINT (c);
+  digit_value = char_decimal_value (c);
   if (digit_value >= radix)
     return SCM_BOOL_F;
 
@@ -2598,21 +3878,21 @@ mem2uinteger (SCM mem, unsigned int *p_idx,
   while (idx != len)
     {
       scm_t_wchar c = scm_i_string_ref (mem, idx);
-      if (uc_is_property_ascii_hex_digit ((scm_t_uint32) c))
-       {
-         if (hash_seen)
-           break;
-         digit_value = XDIGIT2UINT (c);
-         if (digit_value >= radix)
-           break;
-       }
-      else if (c == '#')
+      if (c == '#')
        {
          hash_seen = 1;
          digit_value = 0;
        }
+      else if (hash_seen)
+        break;
       else
-       break;
+        {
+          digit_value = char_decimal_value (c);
+          /* This check catches non-decimals in addition to out-of-range
+             decimals.  */
+          if (digit_value >= radix)
+           break;
+       }
 
       idx++;
       if (SCM_MOST_POSITIVE_FIXNUM / radix < shift)
@@ -2669,7 +3949,7 @@ mem2decimal_from_point (SCM result, SCM mem,
       scm_t_bits shift = 1;
       scm_t_bits add = 0;
       unsigned int digit_value;
-      SCM big_shift = SCM_I_MAKINUM (1);
+      SCM big_shift = SCM_INUM1;
 
       idx++;
       while (idx != len)
@@ -2857,7 +4137,7 @@ mem2ureal (SCM mem, unsigned int *p_idx,
       else if (!uc_is_property_decimal_digit ((scm_t_uint32) scm_i_string_ref (mem, idx+1)))
        return SCM_BOOL_F;
       else
-       result = mem2decimal_from_point (SCM_I_MAKINUM (0), mem,
+       result = mem2decimal_from_point (SCM_INUM0, mem,
                                         p_idx, &x);
     }
   else
@@ -2908,7 +4188,7 @@ mem2ureal (SCM mem, unsigned int *p_idx,
   /* When returning an inexact zero, make sure it is represented as a
      floating point value so that we can change its sign. 
   */
-  if (scm_is_eq (result, SCM_I_MAKINUM(0)) && *p_exactness == INEXACT)
+  if (scm_is_eq (result, SCM_INUM0) && *p_exactness == INEXACT)
     result = scm_from_double (0.0);
 
   return result;
@@ -2959,7 +4239,7 @@ mem2complex (SCM mem, unsigned int idx,
          if (idx != len)
            return SCM_BOOL_F;
          
-         return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign));
+         return scm_make_rectangular (SCM_INUM0, SCM_I_MAKINUM (sign));
        }
       else
        return SCM_BOOL_F;
@@ -2983,7 +4263,7 @@ mem2complex (SCM mem, unsigned int idx,
            return SCM_BOOL_F;
          if (idx != len)
            return SCM_BOOL_F;
-         return scm_make_rectangular (SCM_I_MAKINUM (0), ureal);
+         return scm_make_rectangular (SCM_INUM0, ureal);
 
        case '@':
          /* polar input: <real>@<real>. */
@@ -3041,7 +4321,7 @@ mem2complex (SCM mem, unsigned int idx,
 
              if (scm_is_false (imag))
                imag = SCM_I_MAKINUM (sign);
-             else if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
+             else if (sign == -1 && scm_is_false (scm_nan_p (imag)))
                imag = scm_difference (imag, SCM_UNDEFINED);
 
              if (idx == len)
@@ -3194,40 +4474,6 @@ SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
 /*** END strs->nums ***/
 
 
-SCM
-scm_bigequal (SCM x, SCM y)
-{
-  int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
-  scm_remember_upto_here_2 (x, y);
-  return scm_from_bool (0 == result);
-}
-
-SCM
-scm_real_equalp (SCM x, SCM y)
-{
-  return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
-}
-
-SCM
-scm_complex_equalp (SCM x, SCM y)
-{
-  return scm_from_bool (SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y)
-                  && SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y));
-}
-
-SCM
-scm_i_fraction_equalp (SCM x, SCM y)
-{
-  if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x),
-                              SCM_FRACTION_NUMERATOR (y)))
-      || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x),
-                                 SCM_FRACTION_DENOMINATOR (y))))
-    return SCM_BOOL_F;
-  else
-    return SCM_BOOL_T;
-}
-
-
 SCM_DEFINE (scm_number_p, "number?", 1, 0, 0, 
             (SCM x),
            "Return @code{#t} if @var{x} is a number, @code{#f}\n"
@@ -3260,8 +4506,8 @@ SCM_DEFINE (scm_real_p, "real?", 1, 0, 0,
            "fulfilled if @var{x} is an integer number.")
 #define FUNC_NAME s_scm_real_p
 {
-  /* we can't represent irrational numbers. */
-  return scm_rational_p (x);
+  return scm_from_bool
+    (SCM_I_INUMP (x) || SCM_REALP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x));
 }
 #undef FUNC_NAME
 
@@ -3273,18 +4519,12 @@ SCM_DEFINE (scm_rational_p, "rational?", 1, 0, 0,
            "fulfilled if @var{x} is an integer number.")
 #define FUNC_NAME s_scm_rational_p
 {
-  if (SCM_I_INUMP (x))
-    return SCM_BOOL_T;
-  else if (SCM_IMP (x))
-    return SCM_BOOL_F;
-  else if (SCM_BIGP (x))
-    return SCM_BOOL_T;
-  else if (SCM_FRACTIONP (x))
+  if (SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))
     return SCM_BOOL_T;
   else if (SCM_REALP (x))
-    /* due to their limited precision, all floating point numbers are
-       rational as well. */
-    return SCM_BOOL_T;
+    /* due to their limited precision, finite floating point numbers are
+       rational as well. (finite means neither infinity nor a NaN) */
+    return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x)));
   else
     return SCM_BOOL_F;
 }
@@ -3296,53 +4536,48 @@ SCM_DEFINE (scm_integer_p, "integer?", 1, 0, 0,
            "else.")
 #define FUNC_NAME s_scm_integer_p
 {
-  double r;
-  if (SCM_I_INUMP (x))
-    return SCM_BOOL_T;
-  if (SCM_IMP (x))
-    return SCM_BOOL_F;
-  if (SCM_BIGP (x))
+  if (SCM_I_INUMP (x) || SCM_BIGP (x))
     return SCM_BOOL_T;
-  if (!SCM_INEXACTP (x))
-    return SCM_BOOL_F;
-  if (SCM_COMPLEXP (x))
+  else if (SCM_REALP (x))
+    {
+      double val = SCM_REAL_VALUE (x);
+      return scm_from_bool (!isinf (val) && (val == floor (val)));
+    }
+  else
     return SCM_BOOL_F;
-  r = SCM_REAL_VALUE (x);
-  /* +/-inf passes r==floor(r), making those #t */
-  if (r == floor (r))
-    return SCM_BOOL_T;
-  return SCM_BOOL_F;
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_inexact_p, "inexact?", 1, 0, 0, 
-            (SCM x),
-           "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
-           "else.")
-#define FUNC_NAME s_scm_inexact_p
+SCM scm_i_num_eq_p (SCM, SCM, SCM);
+SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p, "=", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return @code{#t} if all parameters are numerically equal.")
+#define FUNC_NAME s_scm_i_num_eq_p
 {
-  if (SCM_INEXACTP (x))
+  if (SCM_UNBNDP (x) || SCM_UNBNDP (y))
     return SCM_BOOL_T;
-  if (SCM_NUMBERP (x))
-    return SCM_BOOL_F;
-  SCM_WRONG_TYPE_ARG (1, x);
+  while (!scm_is_null (rest))
+    {
+      if (scm_is_false (scm_num_eq_p (x, y)))
+        return SCM_BOOL_F;
+      x = y;
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_num_eq_p (x, y);
 }
 #undef FUNC_NAME
-
-
-SCM_GPROC1 (s_eq_p, "=", scm_tc7_rpsubr, scm_num_eq_p, g_eq_p);
-/* "Return @code{#t} if all parameters are numerically equal."  */
 SCM
 scm_num_eq_p (SCM x, SCM y)
 {
  again:
   if (SCM_I_INUMP (x))
     {
-      long xx = SCM_I_INUM (x);
+      scm_t_signed_bits xx = SCM_I_INUM (x);
       if (SCM_I_INUMP (y))
        {
-         long yy = SCM_I_INUM (y);
+         scm_t_signed_bits yy = SCM_I_INUM (y);
          return scm_from_bool (xx == yy);
        }
       else if (SCM_BIGP (y))
@@ -3361,13 +4596,13 @@ scm_num_eq_p (SCM x, SCM y)
              An alternative (for any size system actually) would be to check
              yy is an integer (with floor) and is in range of an inum
              (compare against appropriate powers of 2) then test
-             xx==(long)yy.  It's just a matter of which casts/comparisons
-             might be fastest or easiest for the cpu.  */
+             xx==(scm_t_signed_bits)yy.  It's just a matter of which
+             casts/comparisons might be fastest or easiest for the cpu.  */
 
           double yy = SCM_REAL_VALUE (y);
           return scm_from_bool ((double) xx == yy
                                && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
-                                   || xx == (long) yy));
+                                   || xx == (scm_t_signed_bits) yy));
         }
       else if (SCM_COMPLEXP (y))
        return scm_from_bool (((double) xx == SCM_COMPLEX_REAL (y))
@@ -3375,7 +4610,7 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return SCM_BOOL_F;
       else
-       SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
+       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))
     {
@@ -3390,7 +4625,7 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_REALP (y))
        {
          int cmp;
-         if (xisnan (SCM_REAL_VALUE (y)))
+         if (isnan (SCM_REAL_VALUE (y)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
          scm_remember_upto_here_1 (x);
@@ -3401,7 +4636,7 @@ scm_num_eq_p (SCM x, SCM y)
          int cmp;
          if (0.0 != SCM_COMPLEX_IMAG (y))
            return SCM_BOOL_F;
-         if (xisnan (SCM_COMPLEX_REAL (y)))
+         if (isnan (SCM_COMPLEX_REAL (y)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_COMPLEX_REAL (y));
          scm_remember_upto_here_1 (x);
@@ -3410,7 +4645,7 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return SCM_BOOL_F;
       else
-       SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
+       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))
     {
@@ -3418,15 +4653,15 @@ scm_num_eq_p (SCM x, SCM y)
       if (SCM_I_INUMP (y))
         {
           /* see comments with inum/real above */
-          long yy = SCM_I_INUM (y);
+          scm_t_signed_bits yy = SCM_I_INUM (y);
           return scm_from_bool (xx == (double) yy
                                && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
-                                   || (long) xx == yy));
+                                   || (scm_t_signed_bits) xx == yy));
         }
       else if (SCM_BIGP (y))
        {
          int cmp;
-         if (xisnan (SCM_REAL_VALUE (x)))
+         if (isnan (SCM_REAL_VALUE (x)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
          scm_remember_upto_here_1 (y);
@@ -3440,15 +4675,15 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
         {
           double  xx = SCM_REAL_VALUE (x);
-          if (xisnan (xx))
+          if (isnan (xx))
             return SCM_BOOL_F;
-          if (xisinf (xx))
+          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_eq_p, x, y, SCM_ARGn, s_eq_p);
+       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))
     {
@@ -3460,7 +4695,7 @@ scm_num_eq_p (SCM x, SCM y)
          int cmp;
          if (0.0 != SCM_COMPLEX_IMAG (x))
            return SCM_BOOL_F;
-         if (xisnan (SCM_COMPLEX_REAL (x)))
+         if (isnan (SCM_COMPLEX_REAL (x)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_COMPLEX_REAL (x));
          scm_remember_upto_here_1 (y);
@@ -3478,15 +4713,15 @@ scm_num_eq_p (SCM x, SCM y)
           if (SCM_COMPLEX_IMAG (x) != 0.0)
             return SCM_BOOL_F;
           xx = SCM_COMPLEX_REAL (x);
-          if (xisnan (xx))
+          if (isnan (xx))
             return SCM_BOOL_F;
-          if (xisinf (xx))
+          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_eq_p, x, y, SCM_ARGn, s_eq_p);
+       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))
     {
@@ -3497,9 +4732,9 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_REALP (y))
         {
           double yy = SCM_REAL_VALUE (y);
-          if (xisnan (yy))
+          if (isnan (yy))
             return SCM_BOOL_F;
-          if (xisinf (yy))
+          if (isinf (yy))
             return scm_from_bool (0.0 < yy);
           y = scm_inexact_to_exact (y);  /* with y as frac or int */
           goto again;
@@ -3510,9 +4745,9 @@ 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 (xisnan (yy))
+          if (isnan (yy))
             return SCM_BOOL_F;
-          if (xisinf (yy))
+          if (isinf (yy))
             return scm_from_bool (0.0 < yy);
           y = scm_inexact_to_exact (y);  /* with y as frac or int */
           goto again;
@@ -3520,10 +4755,10 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
        return scm_i_fraction_equalp (x, y);
       else
-       SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
+       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_eq_p, x, y, SCM_ARG1, s_eq_p);
+    SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARG1, s_scm_i_num_eq_p);
 }
 
 
@@ -3533,20 +4768,36 @@ scm_num_eq_p (SCM x, SCM y)
    mpq_cmp.  flonum/frac compares likewise, but with the slight complication
    of the float exponent to take into account.  */
 
-SCM_GPROC1 (s_less_p, "<", scm_tc7_rpsubr, scm_less_p, g_less_p);
-/* "Return @code{#t} if the list of parameters is monotonically\n"
- * "increasing."
- */
+SCM_INTERNAL SCM scm_i_num_less_p (SCM, SCM, SCM);
+SCM_PRIMITIVE_GENERIC (scm_i_num_less_p, "<", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return @code{#t} if the list of parameters is monotonically\n"
+                       "increasing.")
+#define FUNC_NAME s_scm_i_num_less_p
+{
+  if (SCM_UNBNDP (x) || SCM_UNBNDP (y))
+    return SCM_BOOL_T;
+  while (!scm_is_null (rest))
+    {
+      if (scm_is_false (scm_less_p (x, y)))
+        return SCM_BOOL_F;
+      x = y;
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_less_p (x, y);
+}
+#undef FUNC_NAME
 SCM
 scm_less_p (SCM x, SCM y)
 {
  again:
   if (SCM_I_INUMP (x))
     {
-      long xx = SCM_I_INUM (x);
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_I_INUMP (y))
        {
-         long yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          return scm_from_bool (xx < yy);
        }
       else if (SCM_BIGP (y))
@@ -3566,7 +4817,7 @@ scm_less_p (SCM x, SCM y)
           goto again;
         }
       else
-       SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
+       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))
     {
@@ -3585,7 +4836,7 @@ scm_less_p (SCM x, SCM y)
       else if (SCM_REALP (y))
        {
          int cmp;
-         if (xisnan (SCM_REAL_VALUE (y)))
+         if (isnan (SCM_REAL_VALUE (y)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
          scm_remember_upto_here_1 (x);
@@ -3594,7 +4845,7 @@ scm_less_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
         goto int_frac;
       else
-       SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
+       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))
     {
@@ -3603,7 +4854,7 @@ scm_less_p (SCM x, SCM y)
       else if (SCM_BIGP (y))
        {
          int cmp;
-         if (xisnan (SCM_REAL_VALUE (x)))
+         if (isnan (SCM_REAL_VALUE (x)))
            return SCM_BOOL_F;
          cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
          scm_remember_upto_here_1 (y);
@@ -3614,15 +4865,15 @@ scm_less_p (SCM x, SCM y)
       else if (SCM_FRACTIONP (y))
         {
           double  xx = SCM_REAL_VALUE (x);
-         if (xisnan (xx))
+         if (isnan (xx))
            return SCM_BOOL_F;
-          if (xisinf (xx))
+          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_less_p, x, y, SCM_ARGn, s_less_p);
+       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))
     {
@@ -3636,9 +4887,9 @@ scm_less_p (SCM x, SCM y)
       else if (SCM_REALP (y))
         {
           double yy = SCM_REAL_VALUE (y);
-          if (xisnan (yy))
+          if (isnan (yy))
             return SCM_BOOL_F;
-          if (xisinf (yy))
+          if (isinf (yy))
             return scm_from_bool (0.0 < yy);
           y = scm_inexact_to_exact (y);  /* with y as frac or int */
           goto again;
@@ -3655,43 +4906,75 @@ scm_less_p (SCM x, SCM y)
           goto again;
         }
       else
-       SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
+       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_less_p, x, y, SCM_ARG1, s_less_p);
+    SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARG1, s_scm_i_num_less_p);
 }
 
 
-SCM_GPROC1 (s_scm_gr_p, ">", scm_tc7_rpsubr, scm_gr_p, g_gr_p);
-/* "Return @code{#t} if the list of parameters is monotonically\n"
- * "decreasing."
- */
-#define FUNC_NAME s_scm_gr_p
+SCM scm_i_num_gr_p (SCM, SCM, SCM);
+SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p, ">", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return @code{#t} if the list of parameters is monotonically\n"
+                       "decreasing.")
+#define FUNC_NAME s_scm_i_num_gr_p
+{
+  if (SCM_UNBNDP (x) || SCM_UNBNDP (y))
+    return SCM_BOOL_T;
+  while (!scm_is_null (rest))
+    {
+      if (scm_is_false (scm_gr_p (x, y)))
+        return SCM_BOOL_F;
+      x = y;
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_gr_p (x, y);
+}
+#undef FUNC_NAME
+#define FUNC_NAME s_scm_i_num_gr_p
 SCM
 scm_gr_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_gr_p, x, y, SCM_ARG1, FUNC_NAME);
+    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_gr_p, x, y, SCM_ARG2, FUNC_NAME);
+    SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p, x, y, SCM_ARG2, FUNC_NAME);
   else
     return scm_less_p (y, x);
 }
 #undef FUNC_NAME
 
 
-SCM_GPROC1 (s_scm_leq_p, "<=", scm_tc7_rpsubr, scm_leq_p, g_leq_p);
-/* "Return @code{#t} if the list of parameters is monotonically\n"
- * "non-decreasing."
- */
-#define FUNC_NAME s_scm_leq_p
+SCM scm_i_num_leq_p (SCM, SCM, SCM);
+SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p, "<=", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return @code{#t} if the list of parameters is monotonically\n"
+                       "non-decreasing.")
+#define FUNC_NAME s_scm_i_num_leq_p
+{
+  if (SCM_UNBNDP (x) || SCM_UNBNDP (y))
+    return SCM_BOOL_T;
+  while (!scm_is_null (rest))
+    {
+      if (scm_is_false (scm_leq_p (x, y)))
+        return SCM_BOOL_F;
+      x = y;
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_leq_p (x, y);
+}
+#undef FUNC_NAME
+#define FUNC_NAME s_scm_i_num_leq_p
 SCM
 scm_leq_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_leq_p, x, y, SCM_ARG1, FUNC_NAME);
+    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_leq_p, x, y, SCM_ARG2, FUNC_NAME);
+    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
@@ -3700,18 +4983,34 @@ scm_leq_p (SCM x, SCM y)
 #undef FUNC_NAME
 
 
-SCM_GPROC1 (s_scm_geq_p, ">=", scm_tc7_rpsubr, scm_geq_p, g_geq_p);
-/* "Return @code{#t} if the list of parameters is monotonically\n"
- * "non-increasing."
- */
-#define FUNC_NAME s_scm_geq_p
+SCM scm_i_num_geq_p (SCM, SCM, SCM);
+SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p, ">=", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return @code{#t} if the list of parameters is monotonically\n"
+                       "non-increasing.")
+#define FUNC_NAME s_scm_i_num_geq_p
+{
+  if (SCM_UNBNDP (x) || SCM_UNBNDP (y))
+    return SCM_BOOL_T;
+  while (!scm_is_null (rest))
+    {
+      if (scm_is_false (scm_geq_p (x, y)))
+        return SCM_BOOL_F;
+      x = y;
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_geq_p (x, y);
+}
+#undef FUNC_NAME
+#define FUNC_NAME s_scm_i_num_geq_p
 SCM
 scm_geq_p (SCM x, SCM y)
 {
   if (!SCM_NUMBERP (x))
-    SCM_WTA_DISPATCH_2 (g_geq_p, x, y, SCM_ARG1, FUNC_NAME);
+    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_geq_p, x, y, SCM_ARG2, FUNC_NAME);
+    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
@@ -3720,12 +5019,11 @@ scm_geq_p (SCM x, SCM y)
 #undef FUNC_NAME
 
 
-SCM_GPROC (s_zero_p, "zero?", 1, 0, 0, scm_zero_p, g_zero_p);
-/* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
- * "zero."
- */
-SCM
-scm_zero_p (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_zero_p, "zero?", 1, 0, 0,
+                      (SCM z),
+       "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
+       "zero.")
+#define FUNC_NAME s_scm_zero_p
 {
   if (SCM_I_INUMP (z))
     return scm_from_bool (scm_is_eq (z, SCM_INUM0));
@@ -3739,16 +5037,16 @@ scm_zero_p (SCM z)
   else if (SCM_FRACTIONP (z))
     return SCM_BOOL_F;
   else
-    SCM_WTA_DISPATCH_1 (g_zero_p, z, SCM_ARG1, s_zero_p);
+    SCM_WTA_DISPATCH_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_positive_p, "positive?", 1, 0, 0, scm_positive_p, g_positive_p);
-/* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
- * "zero."
- */
-SCM
-scm_positive_p (SCM x)
+SCM_PRIMITIVE_GENERIC (scm_positive_p, "positive?", 1, 0, 0,
+                      (SCM x),
+       "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
+       "zero.")
+#define FUNC_NAME s_scm_positive_p
 {
   if (SCM_I_INUMP (x))
     return scm_from_bool (SCM_I_INUM (x) > 0);
@@ -3763,16 +5061,16 @@ scm_positive_p (SCM x)
   else if (SCM_FRACTIONP (x))
     return scm_positive_p (SCM_FRACTION_NUMERATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_positive_p, x, SCM_ARG1, s_positive_p);
+    SCM_WTA_DISPATCH_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_negative_p, "negative?", 1, 0, 0, scm_negative_p, g_negative_p);
-/* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
- * "zero."
- */
-SCM
-scm_negative_p (SCM x)
+SCM_PRIMITIVE_GENERIC (scm_negative_p, "negative?", 1, 0, 0,
+                      (SCM x),
+       "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
+       "zero.")
+#define FUNC_NAME s_scm_negative_p
 {
   if (SCM_I_INUMP (x))
     return scm_from_bool (SCM_I_INUM (x) < 0);
@@ -3787,8 +5085,9 @@ scm_negative_p (SCM x)
   else if (SCM_FRACTIONP (x))
     return scm_negative_p (SCM_FRACTION_NUMERATOR (x));
   else
-    SCM_WTA_DISPATCH_1 (g_negative_p, x, SCM_ARG1, s_negative_p);
+    SCM_WTA_DISPATCH_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p);
 }
+#undef FUNC_NAME
 
 
 /* scm_min and scm_max return an inexact when either argument is inexact, as
@@ -3829,10 +5128,10 @@ scm_max (SCM x, SCM y)
   
   if (SCM_I_INUMP (x))
     {
-      long xx = SCM_I_INUM (x);
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_I_INUMP (y))
        {
-         long yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          return (xx < yy) ? y : x;
        }
       else if (SCM_BIGP (y))
@@ -3905,7 +5204,7 @@ scm_max (SCM x, SCM y)
             calling isnan is unavoidable, since it's the only way to know
             which of x or y causes any compares to be false */
          double xx = SCM_REAL_VALUE (x);
-         return (xisnan (xx) || xx > SCM_REAL_VALUE (y)) ? x : y;
+         return (isnan (xx) || xx > SCM_REAL_VALUE (y)) ? x : y;
        }
       else if (SCM_FRACTIONP (y))
        {
@@ -3975,10 +5274,10 @@ scm_min (SCM x, SCM y)
   
   if (SCM_I_INUMP (x))
     {
-      long xx = SCM_I_INUM (x);
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_I_INUMP (y))
        {
-         long yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          return (xx < yy) ? x : y;
        }
       else if (SCM_BIGP (y))
@@ -4051,7 +5350,7 @@ scm_min (SCM x, SCM y)
             calling isnan is unavoidable, since it's the only way to know
             which of x or y causes any compares to be false */
          double xx = SCM_REAL_VALUE (x);
-         return (xisnan (xx) || xx < SCM_REAL_VALUE (y)) ? x : y;
+         return (isnan (xx) || xx < SCM_REAL_VALUE (y)) ? x : y;
        }
       else if (SCM_FRACTIONP (y))
        {
@@ -4121,10 +5420,10 @@ scm_sum (SCM x, SCM y)
     {
       if (SCM_LIKELY (SCM_I_INUMP (y)))
         {
-          long xx = SCM_I_INUM (x);
-          long yy = SCM_I_INUM (y);
-          long int z = xx + yy;
-          return SCM_FIXABLE (z) ? SCM_I_MAKINUM (z) : scm_i_long2big (z);
+          scm_t_inum xx = SCM_I_INUM (x);
+          scm_t_inum yy = SCM_I_INUM (y);
+          scm_t_inum z = xx + yy;
+          return SCM_FIXABLE (z) ? SCM_I_MAKINUM (z) : scm_i_inum2big (z);
         }
       else if (SCM_BIGP (y))
         {
@@ -4133,12 +5432,12 @@ scm_sum (SCM x, SCM y)
         }
       else if (SCM_REALP (y))
         {
-          long int xx = SCM_I_INUM (x);
+          scm_t_inum xx = SCM_I_INUM (x);
           return scm_from_double (xx + SCM_REAL_VALUE (y));
         }
       else if (SCM_COMPLEXP (y))
         {
-          long int xx = SCM_I_INUM (x);
+          scm_t_inum xx = SCM_I_INUM (x);
           return scm_c_make_rectangular (xx + SCM_COMPLEX_REAL (y),
                                    SCM_COMPLEX_IMAG (y));
         }
@@ -4152,7 +5451,7 @@ scm_sum (SCM x, SCM y)
       {
        if (SCM_I_INUMP (y))
          {
-           long int inum;
+           scm_t_inum inum;
            int bigsgn;
          add_big_inum:
            inum = SCM_I_INUM (y);      
@@ -4291,7 +5590,7 @@ SCM_DEFINE (scm_oneplus, "1+", 1, 0, 0,
            "Return @math{@var{x}+1}.")
 #define FUNC_NAME s_scm_oneplus
 {
-  return scm_sum (x, SCM_I_MAKINUM (1));
+  return scm_sum (x, SCM_INUM1);
 }
 #undef FUNC_NAME
 
@@ -4326,11 +5625,11 @@ scm_difference (SCM x, SCM y)
       else 
         if (SCM_I_INUMP (x))
           {
-            long xx = -SCM_I_INUM (x);
+            scm_t_inum xx = -SCM_I_INUM (x);
             if (SCM_FIXABLE (xx))
               return SCM_I_MAKINUM (xx);
             else
-              return scm_i_long2big (xx);
+              return scm_i_inum2big (xx);
           }
         else if (SCM_BIGP (x))
           /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
@@ -4352,21 +5651,25 @@ scm_difference (SCM x, SCM y)
     {
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long int xx = SCM_I_INUM (x);
-         long int yy = SCM_I_INUM (y);
-         long int z = xx - yy;
+         scm_t_inum xx = SCM_I_INUM (x);
+         scm_t_inum yy = SCM_I_INUM (y);
+         scm_t_inum z = xx - yy;
          if (SCM_FIXABLE (z))
            return SCM_I_MAKINUM (z);
          else
-           return scm_i_long2big (z);
+           return scm_i_inum2big (z);
        }
       else if (SCM_BIGP (y))
        {
          /* inum-x - big-y */
-         long xx = SCM_I_INUM (x);
+         scm_t_inum xx = SCM_I_INUM (x);
 
          if (xx == 0)
-           return scm_i_clonebig (y, 0);
+           {
+             /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
+                bignum, but negating that gives a fixnum.  */
+             return scm_i_normbig (scm_i_clonebig (y, 0));
+           }
          else
            {
              int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
@@ -4391,12 +5694,12 @@ scm_difference (SCM x, SCM y)
        }
       else if (SCM_REALP (y))
        {
-         long int xx = SCM_I_INUM (x);
+         scm_t_inum xx = SCM_I_INUM (x);
          return scm_from_double (xx - SCM_REAL_VALUE (y));
        }
       else if (SCM_COMPLEXP (y))
        {
-         long int xx = SCM_I_INUM (x);
+         scm_t_inum xx = SCM_I_INUM (x);
          return scm_c_make_rectangular (xx - SCM_COMPLEX_REAL (y),
                                   - SCM_COMPLEX_IMAG (y));
        }
@@ -4413,13 +5716,13 @@ scm_difference (SCM x, SCM y)
       if (SCM_I_INUMP (y))
        {
          /* big-x - inum-y */
-         long yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
 
          scm_remember_upto_here_1 (x);
          if (sgn_x == 0)
            return (SCM_FIXABLE (-yy) ?
-                   SCM_I_MAKINUM (-yy) : scm_from_long (-yy));
+                   SCM_I_MAKINUM (-yy) : scm_from_inum (-yy));
          else
            {
              SCM result = scm_i_mkbig ();
@@ -4551,7 +5854,7 @@ SCM_DEFINE (scm_oneminus, "1-", 1, 0, 0,
            "Return @math{@var{x}-1}.")
 #define FUNC_NAME s_scm_oneminus
 {
-  return scm_difference (x, SCM_I_MAKINUM (1));
+  return scm_difference (x, SCM_INUM1);
 }
 #undef FUNC_NAME
 
@@ -4589,7 +5892,7 @@ scm_product (SCM x, SCM y)
   
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      long xx;
+      scm_t_inum xx;
 
     intbig:
       xx = SCM_I_INUM (x);
@@ -4598,18 +5901,29 @@ scm_product (SCM x, SCM y)
        {
         case 0: return x; break;
         case 1: return y; break;
+         /*
+          * The following case (x = -1) is important for more than
+          * just optimization.  It handles the case of negating
+          * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
+          * which is a bignum that must be changed back into a fixnum.
+          * Failure to do so will cause the following to return #f:
+          * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
+          */
+        case -1:
+         return scm_difference(y, SCM_UNDEFINED);
+         break;
        }
 
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
-         long kk = xx * yy;
+         scm_t_inum yy = SCM_I_INUM (y);
+         scm_t_inum kk = xx * yy;
          SCM k = SCM_I_MAKINUM (kk);
          if ((kk == SCM_I_INUM (k)) && (kk / xx == yy))
            return k;
          else
            {
-             SCM result = scm_i_long2big (xx);
+             SCM result = scm_i_inum2big (xx);
              mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), yy);
              return scm_i_normbig (result);
            }
@@ -4821,7 +6135,7 @@ do_divide (SCM x, SCM y, int inexact)
        SCM_WTA_DISPATCH_0 (g_divide, s_divide);
       else if (SCM_I_INUMP (x))
        {
-         long xx = SCM_I_INUM (x);
+         scm_t_inum xx = SCM_I_INUM (x);
          if (xx == 1 || xx == -1)
            return x;
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
@@ -4832,14 +6146,14 @@ do_divide (SCM x, SCM y, int inexact)
            {
              if (inexact)
                return scm_from_double (1.0 / (double) xx);
-             else return scm_i_make_ratio (SCM_I_MAKINUM(1), x);
+             else return scm_i_make_ratio (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 (SCM_I_MAKINUM(1), x);
+         else return scm_i_make_ratio (SCM_INUM1, x);
        }
       else if (SCM_REALP (x))
        {
@@ -4877,10 +6191,10 @@ do_divide (SCM x, SCM y, int inexact)
 
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      long xx = SCM_I_INUM (x);
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
-         long yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          if (yy == 0)
            {
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
@@ -4897,11 +6211,11 @@ do_divide (SCM x, SCM y, int inexact)
            }
          else
            {
-             long z = xx / yy;
+             scm_t_inum z = xx / yy;
              if (SCM_FIXABLE (z))
                return SCM_I_MAKINUM (z);
              else
-               return scm_i_long2big (z);
+               return scm_i_inum2big (z);
            }
        }
       else if (SCM_BIGP (y))
@@ -4952,7 +6266,7 @@ do_divide (SCM x, SCM y, int inexact)
     {
       if (SCM_I_INUMP (y))
        {
-         long int yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
          if (yy == 0)
            {
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
@@ -4975,7 +6289,7 @@ do_divide (SCM x, SCM y, int inexact)
                 middle ground: test, then if divisible, use the faster div
                 func. */
 
-             long abs_yy = yy < 0 ? -yy : yy;
+             scm_t_inum abs_yy = yy < 0 ? -yy : yy;
              int divisible_p = mpz_divisible_ui_p (SCM_I_BIG_MPZ (x), abs_yy);
 
              if (divisible_p)
@@ -4997,47 +6311,33 @@ do_divide (SCM x, SCM y, int inexact)
        }
       else if (SCM_BIGP (y))
        {
-         int y_is_zero = (mpz_sgn (SCM_I_BIG_MPZ (y)) == 0);
-         if (y_is_zero)
+         /* big_x / big_y */
+         if (inexact)
            {
-#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
-             scm_num_overflow (s_divide);
-#else
-             int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
-             scm_remember_upto_here_1 (x);
-             return (sgn == 0) ? scm_nan () : scm_inf ();
-#endif
+             /* 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
            {
-             /* 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))
@@ -5066,7 +6366,7 @@ do_divide (SCM x, SCM y, int inexact)
       double rx = SCM_REAL_VALUE (x);
       if (SCM_I_INUMP (y))
        {
-         long int yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
          if (yy == 0)
            scm_num_overflow (s_divide);
@@ -5106,7 +6406,7 @@ do_divide (SCM x, SCM y, int inexact)
       double ix = SCM_COMPLEX_IMAG (x);
       if (SCM_I_INUMP (y))
        {
-         long int yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
          if (yy == 0)
            scm_num_overflow (s_divide);
@@ -5162,7 +6462,7 @@ do_divide (SCM x, SCM y, int inexact)
     {
       if (SCM_I_INUMP (y)) 
        {
-         long int yy = SCM_I_INUM (y);
+         scm_t_inum yy = SCM_I_INUM (y);
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
          if (yy == 0)
            scm_num_overflow (s_divide);
@@ -5280,8 +6580,6 @@ SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-static SCM exactly_one_half;
-
 SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
            (SCM x),
            "Round the number @var{x} towards the nearest integer. "
@@ -5303,7 +6601,7 @@ SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
       /* Adjust so that the rounding is towards even.  */
       if (scm_is_true (scm_num_eq_p (plus_half, result))
           && scm_is_true (scm_odd_p (result)))
-        return scm_difference (result, SCM_I_MAKINUM (1));
+        return scm_difference (result, SCM_INUM1);
       else
         return result;
     }
@@ -5333,7 +6631,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0,
          /* For negative x, we need to return q-1 unless x is an
             integer.  But fractions are never integer, per our
             assumptions. */
-         return scm_difference (q, SCM_I_MAKINUM (1));
+         return scm_difference (q, SCM_INUM1);
        }
     }
   else
@@ -5364,7 +6662,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
          /* For positive x, we need to return q+1 unless x is an
             integer.  But fractions are never integer, per our
             assumptions. */
-         return scm_sum (q, SCM_I_MAKINUM (1));
+         return scm_sum (q, SCM_INUM1);
        }
     }
   else
@@ -5378,19 +6676,44 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
    Written by Jerry D. Hedden, (C) FSF.
    See the file `COPYING' for terms applying to this program. */
 
-SCM_DEFINE (scm_expt, "expt", 2, 0, 0,
-            (SCM x, SCM y),
-           "Return @var{x} raised to the power of @var{y}.") 
+SCM_PRIMITIVE_GENERIC (scm_expt, "expt", 2, 0, 0,
+                      (SCM x, SCM y),
+                      "Return @var{x} raised to the power of @var{y}.")
 #define FUNC_NAME s_scm_expt
 {
-  if (!SCM_INEXACTP (y) && scm_is_integer (y))
-    return scm_integer_expt (x, y);
+  if (scm_is_integer (y))
+    {
+      if (scm_is_true (scm_exact_p (y)))
+       return scm_integer_expt (x, y);
+      else
+       {
+         /* Here we handle the case where the exponent is an inexact
+            integer.  We make the exponent exact in order to use
+            scm_integer_expt, and thus avoid the spurious imaginary
+            parts that may result from round-off errors in the general
+            e^(y log x) method below (for example when squaring a large
+            negative number).  In this case, we must return an inexact
+            result for correctness.  We also make the base inexact so
+            that scm_integer_expt will use fast inexact arithmetic
+            internally.  Note that making the base inexact is not
+            sufficient to guarantee an inexact result, because
+            scm_integer_expt will return an exact 1 when the exponent
+            is 0, even if the base is inexact. */
+         return scm_exact_to_inexact
+           (scm_integer_expt (scm_exact_to_inexact (x),
+                              scm_inexact_to_exact (y)));
+       }
+    }
   else if (scm_is_real (x) && scm_is_real (y) && scm_to_double (x) >= 0.0)
     {
       return scm_from_double (pow (scm_to_double (x), scm_to_double (y)));
     }
-  else
+  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);
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt);
 }
 #undef FUNC_NAME
 
@@ -5615,7 +6938,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_asinh, "asinh", 1, 0, 0,
   else if (scm_is_number (z))
     return scm_log (scm_sum (z,
                              scm_sqrt (scm_sum (scm_product (z, z),
-                                                SCM_I_MAKINUM (1)))));
+                                                SCM_INUM1))));
   else
     SCM_WTA_DISPATCH_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh);
 }
@@ -5631,7 +6954,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_acosh, "acosh", 1, 0, 0,
   else if (scm_is_number (z))
     return scm_log (scm_sum (z,
                              scm_sqrt (scm_difference (scm_product (z, z),
-                                                       SCM_I_MAKINUM (1)))));
+                                                       SCM_INUM1))));
   else
     SCM_WTA_DISPATCH_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh);
 }
@@ -5645,8 +6968,8 @@ SCM_PRIMITIVE_GENERIC (scm_sys_atanh, "atanh", 1, 0, 0,
   if (scm_is_real (z) && scm_to_double (z) >= -1.0 && scm_to_double (z) <= 1.0)
     return scm_from_double (atanh (scm_to_double (z)));
   else if (scm_is_number (z))
-    return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z),
-                                            scm_difference (SCM_I_MAKINUM (1), z))),
+    return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1, z),
+                                            scm_difference (SCM_INUM1, z))),
                        SCM_I_MAKINUM (2));
   else
     SCM_WTA_DISPATCH_1 (g_scm_sys_atanh, z, 1, s_scm_sys_atanh);
@@ -5661,9 +6984,10 @@ scm_c_make_rectangular (double re, double im)
   else
     {
       SCM z;
-      SCM_NEWSMOB (z, scm_tc16_complex,
-                  scm_gc_malloc_pointerless (sizeof (scm_t_complex),
+
+      z = PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex),
                                              "complex"));
+      SCM_SET_CELL_TYPE (z, scm_tc16_complex);
       SCM_COMPLEX_REAL (z) = re;
       SCM_COMPLEX_IMAG (z) = im;
       return z;
@@ -5715,100 +7039,86 @@ SCM_DEFINE (scm_make_polar, "make-polar", 2, 0, 0,
 #undef FUNC_NAME
 
 
-SCM_GPROC (s_real_part, "real-part", 1, 0, 0, scm_real_part, g_real_part);
-/* "Return the real part of the number @var{z}."
- */
-SCM
-scm_real_part (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_real_part, "real-part", 1, 0, 0,
+                      (SCM z),
+                      "Return the real part of the number @var{z}.")
+#define FUNC_NAME s_scm_real_part
 {
-  if (SCM_I_INUMP (z))
-    return z;
-  else if (SCM_BIGP (z))
-    return z;
-  else if (SCM_REALP (z))
-    return z;
-  else if (SCM_COMPLEXP (z))
+  if (SCM_COMPLEXP (z))
     return scm_from_double (SCM_COMPLEX_REAL (z));
-  else if (SCM_FRACTIONP (z))
+  else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z))
     return z;
   else
-    SCM_WTA_DISPATCH_1 (g_real_part, z, SCM_ARG1, s_real_part);
+    SCM_WTA_DISPATCH_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_imag_part, "imag-part", 1, 0, 0, scm_imag_part, g_imag_part);
-/* "Return the imaginary part of the number @var{z}."
- */
-SCM
-scm_imag_part (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_imag_part, "imag-part", 1, 0, 0,
+                      (SCM z),
+                      "Return the imaginary part of the number @var{z}.")
+#define FUNC_NAME s_scm_imag_part
 {
-  if (SCM_I_INUMP (z))
-    return SCM_INUM0;
-  else if (SCM_BIGP (z))
-    return SCM_INUM0;
-  else if (SCM_REALP (z))
-    return scm_flo0;
-  else if (SCM_COMPLEXP (z))
+  if (SCM_COMPLEXP (z))
     return scm_from_double (SCM_COMPLEX_IMAG (z));
-  else if (SCM_FRACTIONP (z))
+  else if (SCM_REALP (z))
+    return flo0;
+  else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
     return SCM_INUM0;
   else
-    SCM_WTA_DISPATCH_1 (g_imag_part, z, SCM_ARG1, s_imag_part);
+    SCM_WTA_DISPATCH_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part);
 }
+#undef FUNC_NAME
 
-SCM_GPROC (s_numerator, "numerator", 1, 0, 0, scm_numerator, g_numerator);
-/* "Return the numerator of the number @var{z}."
- */
-SCM
-scm_numerator (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_numerator, "numerator", 1, 0, 0,
+                      (SCM z),
+                      "Return the numerator of the number @var{z}.")
+#define FUNC_NAME s_scm_numerator
 {
-  if (SCM_I_INUMP (z))
-    return z;
-  else if (SCM_BIGP (z))
+  if (SCM_I_INUMP (z) || SCM_BIGP (z))
     return z;
   else if (SCM_FRACTIONP (z))
     return SCM_FRACTION_NUMERATOR (z);
   else if (SCM_REALP (z))
     return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
   else
-    SCM_WTA_DISPATCH_1 (g_numerator, z, SCM_ARG1, s_numerator);
+    SCM_WTA_DISPATCH_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_denominator, "denominator", 1, 0, 0, scm_denominator, g_denominator);
-/* "Return the denominator of the number @var{z}."
- */
-SCM
-scm_denominator (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_denominator, "denominator", 1, 0, 0,
+                      (SCM z),
+                      "Return the denominator of the number @var{z}.")
+#define FUNC_NAME s_scm_denominator
 {
-  if (SCM_I_INUMP (z))
-    return SCM_I_MAKINUM (1);
-  else if (SCM_BIGP (z)) 
-    return SCM_I_MAKINUM (1);
+  if (SCM_I_INUMP (z) || SCM_BIGP (z)) 
+    return SCM_INUM1;
   else if (SCM_FRACTIONP (z))
     return SCM_FRACTION_DENOMINATOR (z);
   else if (SCM_REALP (z))
     return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
   else
-    SCM_WTA_DISPATCH_1 (g_denominator, z, SCM_ARG1, s_denominator);
+    SCM_WTA_DISPATCH_1 (g_scm_denominator, z, SCM_ARG1, s_scm_denominator);
 }
+#undef FUNC_NAME
 
-SCM_GPROC (s_magnitude, "magnitude", 1, 0, 0, scm_magnitude, g_magnitude);
-/* "Return the magnitude of the number @var{z}. This is the same as\n"
- * "@code{abs} for real arguments, but also allows complex numbers."
- */
-SCM
-scm_magnitude (SCM z)
+
+SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0,
+                      (SCM z),
+       "Return the magnitude of the number @var{z}. This is the same as\n"
+       "@code{abs} for real arguments, but also allows complex numbers.")
+#define FUNC_NAME s_scm_magnitude
 {
   if (SCM_I_INUMP (z))
     {
-      long int zz = SCM_I_INUM (z);
+      scm_t_inum zz = SCM_I_INUM (z);
       if (zz >= 0)
        return z;
       else if (SCM_POSFIXABLE (-zz))
        return SCM_I_MAKINUM (-zz);
       else
-       return scm_i_long2big (-zz);
+       return scm_i_inum2big (-zz);
     }
   else if (SCM_BIGP (z))
     {
@@ -5831,24 +7141,24 @@ scm_magnitude (SCM z)
                             SCM_FRACTION_DENOMINATOR (z));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_magnitude, z, SCM_ARG1, s_magnitude);
+    SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_angle, "angle", 1, 0, 0, scm_angle, g_angle);
-/* "Return the angle of the complex number @var{z}."
- */
-SCM
-scm_angle (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0,
+                      (SCM z),
+                      "Return the angle of the complex number @var{z}.")
+#define FUNC_NAME s_scm_angle
 {
   /* atan(0,-1) is pi and it'd be possible to have that as a constant like
-     scm_flo0 to save allocating a new flonum with scm_from_double each time.
+     flo0 to save allocating a new flonum with scm_from_double each time.
      But if atan2 follows the floating point rounding mode, then the value
      is not a constant.  Maybe it'd be close enough though.  */
   if (SCM_I_INUMP (z))
     {
       if (SCM_I_INUM (z) >= 0)
-        return scm_flo0;
+        return flo0;
       else
        return scm_from_double (atan2 (0.0, -1.0));
     }
@@ -5859,12 +7169,12 @@ scm_angle (SCM z)
       if (sgn < 0)
        return scm_from_double (atan2 (0.0, -1.0));
       else
-        return scm_flo0;
+        return flo0;
     }
   else if (SCM_REALP (z))
     {
       if (SCM_REAL_VALUE (z) >= 0)
-        return scm_flo0;
+        return flo0;
       else
         return scm_from_double (atan2 (0.0, -1.0));
     }
@@ -5873,19 +7183,19 @@ scm_angle (SCM z)
   else if (SCM_FRACTIONP (z))
     {
       if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
-       return scm_flo0;
+       return flo0;
       else return scm_from_double (atan2 (0.0, -1.0));
     }
   else
-    SCM_WTA_DISPATCH_1 (g_angle, z, SCM_ARG1, s_angle);
+    SCM_WTA_DISPATCH_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
 }
+#undef FUNC_NAME
 
 
-SCM_GPROC (s_exact_to_inexact, "exact->inexact", 1, 0, 0, scm_exact_to_inexact, g_exact_to_inexact);
-/* Convert the number @var{x} to its inexact representation.\n" 
- */
-SCM
-scm_exact_to_inexact (SCM z)
+SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact, "exact->inexact", 1, 0, 0,
+                      (SCM z),
+       "Convert the number @var{z} to its inexact representation.\n")
+#define FUNC_NAME s_scm_exact_to_inexact
 {
   if (SCM_I_INUMP (z))
     return scm_from_double ((double) SCM_I_INUM (z));
@@ -5896,22 +7206,21 @@ scm_exact_to_inexact (SCM z)
   else if (SCM_INEXACTP (z))
     return z;
   else
-    SCM_WTA_DISPATCH_1 (g_exact_to_inexact, z, 1, s_exact_to_inexact);
+    SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact, z, 1, s_scm_exact_to_inexact);
 }
+#undef FUNC_NAME
 
 
-SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, 
-            (SCM z),
-           "Return an exact number that is numerically closest to @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, 
+                      (SCM z),
+       "Return an exact number that is numerically closest to @var{z}.")
 #define FUNC_NAME s_scm_inexact_to_exact
 {
-  if (SCM_I_INUMP (z))
-    return z;
-  else if (SCM_BIGP (z))
+  if (SCM_I_INUMP (z) || SCM_BIGP (z))
     return z;
   else if (SCM_REALP (z))
     {
-      if (xisinf (SCM_REAL_VALUE (z)) || xisnan (SCM_REAL_VALUE (z)))
+      if (!DOUBLE_IS_FINITE (SCM_REAL_VALUE (z)))
        SCM_OUT_OF_RANGE (1, z);
       else
        {
@@ -5933,7 +7242,7 @@ SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
   else if (SCM_FRACTIONP (z))
     return z;
   else
-    SCM_WRONG_TYPE_ARG (1, z);
+    SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact, z, 1, s_scm_inexact_to_exact);
 }
 #undef FUNC_NAME
 
@@ -5964,9 +7273,9 @@ SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
 
       SCM ex = scm_inexact_to_exact (x);
       SCM int_part = scm_floor (ex);
-      SCM tt = SCM_I_MAKINUM (1);
-      SCM a1 = SCM_I_MAKINUM (0), a2 = SCM_I_MAKINUM (1), a = SCM_I_MAKINUM (0);
-      SCM b1 = SCM_I_MAKINUM (1), b2 = SCM_I_MAKINUM (0), b = SCM_I_MAKINUM (0);
+      SCM tt = SCM_INUM1;
+      SCM a1 = SCM_INUM0, a2 = SCM_INUM1, a = SCM_INUM0;
+      SCM b1 = SCM_INUM1, b2 = SCM_INUM0, b = SCM_INUM0;
       SCM rx;
       int i = 0;
 
@@ -6200,8 +7509,6 @@ scm_i_range_error (SCM bad_val, SCM min, SCM max)
 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
 #include "libguile/conv-integer.i.c"
 
-#if SCM_HAVE_T_INT64
-
 #define TYPE                     scm_t_int64
 #define TYPE_MIN                 SCM_T_INT64_MIN
 #define TYPE_MAX                 SCM_T_INT64_MAX
@@ -6218,8 +7525,6 @@ scm_i_range_error (SCM bad_val, SCM min, SCM max)
 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
 #include "libguile/conv-uinteger.i.c"
 
-#endif
-
 void
 scm_to_mpz (SCM val, mpz_t rop)
 {
@@ -6267,20 +7572,28 @@ scm_to_double (SCM val)
 SCM
 scm_from_double (double val)
 {
-  SCM z = scm_double_cell (scm_tc16_real, 0, 0, 0);
+  SCM z;
+
+  z = PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real"));
+
+  SCM_SET_CELL_TYPE (z, scm_tc16_real);
   SCM_REAL_VALUE (z) = val;
+
   return z;
 }
 
-#if SCM_ENABLE_DISCOURAGED == 1
+#if SCM_ENABLE_DEPRECATED == 1
 
 float
-scm_num2float (SCM num, unsigned long int pos, const char *s_caller)
+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 (!xisinf (res))
+      if (!isinf (res))
        return res;
       else
        scm_out_of_range (NULL, num);
@@ -6290,12 +7603,15 @@ scm_num2float (SCM num, unsigned long int pos, const char *s_caller)
 }
 
 double
-scm_num2double (SCM num, unsigned long int pos, const char *s_caller)
+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 (!xisinf (res))
+      if (!isinf (res))
        return res;
       else
        scm_out_of_range (NULL, num);
@@ -6366,9 +7682,9 @@ scm_is_number (SCM z)
    real-only case, and because we have to test SCM_COMPLEXP anyway so may as
    well use it to go straight to the applicable C func.  */
 
-SCM_DEFINE (scm_log, "log", 1, 0, 0,
-            (SCM z),
-           "Return the natural logarithm of @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0,
+                      (SCM z),
+                      "Return the natural logarithm of @var{z}.")
 #define FUNC_NAME s_scm_log
 {
   if (SCM_COMPLEXP (z))
@@ -6382,7 +7698,7 @@ SCM_DEFINE (scm_log, "log", 1, 0, 0,
                                      atan2 (im, re));
 #endif
     }
-  else
+  else if (SCM_NUMBERP (z))
     {
       /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
          although the value itself overflows.  */
@@ -6393,13 +7709,15 @@ SCM_DEFINE (scm_log, "log", 1, 0, 0,
       else
         return scm_c_make_rectangular (l, M_PI);
     }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log);
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
-            (SCM z),
-           "Return the base 10 logarithm of @var{z}.")
+SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0,
+                      (SCM z),
+                      "Return the base 10 logarithm of @var{z}.")
 #define FUNC_NAME s_scm_log10
 {
   if (SCM_COMPLEXP (z))
@@ -6407,7 +7725,8 @@ SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
       /* Mingw has clog() but not clog10().  (Maybe it'd be worth using
          clog() and a multiply by M_LOG10E, rather than the fallback
          log10+hypot+atan2.)  */
-#if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
+#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
+      && defined SCM_COMPLEX_VALUE
       return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z)));
 #else
       double re = SCM_COMPLEX_REAL (z);
@@ -6416,7 +7735,7 @@ SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
                                      M_LOG10E * atan2 (im, re));
 #endif
     }
-  else
+  else if (SCM_NUMBERP (z))
     {
       /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
          although the value itself overflows.  */
@@ -6427,14 +7746,16 @@ SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
       else
         return scm_c_make_rectangular (l, M_LOG10E * M_PI);
     }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10);
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_exp, "exp", 1, 0, 0,
-            (SCM z),
-           "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
-           "base of natural logarithms (2.71828@dots{}).")
+SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0,
+                      (SCM z),
+       "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
+       "base of natural logarithms (2.71828@dots{}).")
 #define FUNC_NAME s_scm_exp
 {
   if (SCM_COMPLEXP (z))
@@ -6446,50 +7767,55 @@ SCM_DEFINE (scm_exp, "exp", 1, 0, 0,
                                SCM_COMPLEX_IMAG (z));
 #endif
     }
-  else
+  else if (SCM_NUMBERP (z))
     {
       /* When z is a negative bignum the conversion to double overflows,
          giving -infinity, but that's ok, the exp is still 0.0.  */
       return scm_from_double (exp (scm_to_double (z)));
     }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_exp, z, 1, s_scm_exp);
 }
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
-            (SCM x),
-           "Return the square root of @var{z}.  Of the two possible roots\n"
-           "(positive and negative), the one with the a positive real part\n"
-           "is returned, or if that's zero then a positive imaginary part.\n"
-           "Thus,\n"
-           "\n"
-           "@example\n"
-           "(sqrt 9.0)       @result{} 3.0\n"
-           "(sqrt -9.0)      @result{} 0.0+3.0i\n"
-           "(sqrt 1.0+1.0i)  @result{} 1.09868411346781+0.455089860562227i\n"
-           "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
-           "@end example")
+SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
+                      (SCM z),
+       "Return the square root of @var{z}.  Of the two possible roots\n"
+       "(positive and negative), the one with the a positive real part\n"
+       "is returned, or if that's zero then a positive imaginary part.\n"
+       "Thus,\n"
+       "\n"
+       "@example\n"
+       "(sqrt 9.0)       @result{} 3.0\n"
+       "(sqrt -9.0)      @result{} 0.0+3.0i\n"
+       "(sqrt 1.0+1.0i)  @result{} 1.09868411346781+0.455089860562227i\n"
+       "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
+       "@end example")
 #define FUNC_NAME s_scm_sqrt
 {
-  if (SCM_COMPLEXP (x))
+  if (SCM_COMPLEXP (z))
     {
-#if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
-      return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x)));
+#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT   \
+      && defined SCM_COMPLEX_VALUE
+      return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z)));
 #else
-      double re = SCM_COMPLEX_REAL (x);
-      double im = SCM_COMPLEX_IMAG (x);
+      double re = SCM_COMPLEX_REAL (z);
+      double im = SCM_COMPLEX_IMAG (z);
       return scm_c_make_polar (sqrt (hypot (re, im)),
                                0.5 * atan2 (im, re));
 #endif
     }
-  else
+  else if (SCM_NUMBERP (z))
     {
-      double xx = scm_to_double (x);
+      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);
 }
 #undef FUNC_NAME
 
@@ -6513,7 +7839,7 @@ scm_init_numbers ()
 
   scm_add_feature ("complex");
   scm_add_feature ("inexact");
-  scm_flo0 = scm_from_double (0.0);
+  flo0 = scm_from_double (0.0);
 
   /* determine floating point precision */
   for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
@@ -6523,11 +7849,10 @@ scm_init_numbers ()
     }
 #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;
+  scm_dblprec[10-2] = (DBL_DIG > 20) ? 20 : DBL_DIG;
 #endif
 
-  exactly_one_half = scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
-                                                      SCM_I_MAKINUM (2)));
+  exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
 #include "libguile/numbers.x"
 }