Rename {euclidean,centered}_quo_rem to {euclidean,centered}_divide
[bpt/guile.git] / libguile / numbers.c
index 4a458c4..41d178b 100644 (file)
@@ -1,22 +1,23 @@
-/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008 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.
  *
  *
  * This library is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public
- * License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 3 of
+ * the License, or (at your option) any later version.
  *
- * This library is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Lesser General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public
  * License along with this library; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ * 02110-1301 USA
  */
 
 \f
 
  */
 
-/* tell glibc (2.3) to give prototype for C99 trunc(), csqrt(), etc */
-#define _GNU_SOURCE
-
-#if HAVE_CONFIG_H
+#ifdef HAVE_CONFIG_H
 #  include <config.h>
 #endif
 
 #include <math.h>
-#include <ctype.h>
 #include <string.h>
+#include <unicase.h>
+#include <unictype.h>
 
 #if HAVE_COMPLEX_H
 #include <complex.h>
@@ -61,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"
@@ -68,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;
-}
 
+#if !defined (HAVE_ASINH)
+static double asinh (double x) { return log (x + sqrt (x * x + 1)); }
 #endif
+#if !defined (HAVE_ACOSH)
+static double acosh (double x) { return log (x + sqrt (x * x - 1)); }
+#endif
+#if !defined (HAVE_ATANH)
+static double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); }
 #endif
-
 
 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
    an explicit check.  In some future gmp (don't know what version number),
    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
@@ -187,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;
 }
@@ -210,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;
 }
@@ -219,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));
@@ -240,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;
 }
@@ -266,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);
 }
@@ -351,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);
     }
@@ -364,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;
   }
@@ -389,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 
@@ -413,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);
        }
@@ -441,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);
@@ -459,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);
@@ -479,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.")
@@ -506,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))
@@ -515,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.")
@@ -541,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))
@@ -550,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
 
@@ -608,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.  */
 
@@ -620,7 +661,7 @@ guile_ieee_init (void)
      before trying to use it.  (But in practice we believe this is not a
      problem on any system guile is likely to target.)  */
   guile_Inf = INFINITY;
-#elif HAVE_DINFINITY
+#elif defined HAVE_DINFINITY
   /* OSF */
   extern unsigned int DINFINITY[2];
   guile_Inf = (*((double *) (DINFINITY)));
@@ -636,14 +677,10 @@ guile_ieee_init (void)
     }
 #endif
 
-#endif
-
-#if defined (HAVE_ISNAN)
-
 #ifdef NAN
   /* C99 NAN, when available */
   guile_NaN = NAN;
-#elif HAVE_DQNAN
+#elif defined HAVE_DQNAN
   {
     /* OSF */
     extern unsigned int DQNAN[2];
@@ -652,8 +689,6 @@ guile_ieee_init (void)
 #else
   guile_NaN = guile_Inf / guile_Inf;
 #endif
-
-#endif
 }
 
 SCM_DEFINE (scm_inf, "inf", 0, 0, 0, 
@@ -690,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))
     {
@@ -732,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))
@@ -766,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
            {
@@ -806,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);
            }
        }
@@ -843,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 ();
@@ -878,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)
                {
@@ -948,102 +985,1339 @@ scm_modulo (SCM x, SCM y)
                }
              else
                {
-                 result = scm_i_mkbig ();
-                 /* do this after the last scm_op */
-                 mpz_init_set_si (z_x, xx);
-                 mpz_mod (SCM_I_BIG_MPZ (result),
-                          z_x,
-                          SCM_I_BIG_MPZ (y));        
-                 scm_remember_upto_here_1 (y);
+                 result = scm_i_mkbig ();
+                 /* do this after the last scm_op */
+                 mpz_init_set_si (z_x, xx);
+                 mpz_mod (SCM_I_BIG_MPZ (result),
+                          z_x,
+                          SCM_I_BIG_MPZ (y));        
+                 scm_remember_upto_here_1 (y);
+               }
+        
+             if ((sgn_y < 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_1 (y);
+             /* and do this before the next one */
+             mpz_clear (z_x);
+             return scm_i_normbig (result);
+           }
+       }
+      else
+       SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo);
+    }
+  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_modulo);
+         else
+           {
+             SCM result = scm_i_mkbig ();
+             mpz_mod_ui (SCM_I_BIG_MPZ (result),
+                         SCM_I_BIG_MPZ (x),
+                         (yy < 0) ? - yy : yy);
+             scm_remember_upto_here_1 (x);
+             if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
+               mpz_sub_ui (SCM_I_BIG_MPZ (result),
+                           SCM_I_BIG_MPZ (result),
+                           - yy);
+             return scm_i_normbig (result);
+           }
+       }
+      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; }
+                   }
                }
-        
-             if ((sgn_y < 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_1 (y);
-             /* and do this before the next one */
-             mpz_clear (z_x);
-             return scm_i_normbig (result);
+             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_modulo, x, y, SCM_ARG2, s_modulo);
+       SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2,
+                           s_scm_centered_divide);
     }
   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_centered_divide);
          else
            {
-             SCM result = scm_i_mkbig ();
-             mpz_mod_ui (SCM_I_BIG_MPZ (result),
-                         SCM_I_BIG_MPZ (x),
-                         (yy < 0) ? - yy : yy);
-             scm_remember_upto_here_1 (x);
-             if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
-               mpz_sub_ui (SCM_I_BIG_MPZ (result),
-                           SCM_I_BIG_MPZ (result),
-                           - yy);
-             return scm_i_normbig (result);
+             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_GPROC1 (s_gcd, "gcd", scm_tc7_asubr, scm_gcd, g_gcd);
-/* "Return the greatest common divisor of all arguments.\n"
- * "If called without arguments, 0 is returned."
- */
+
+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"
+                       "If called without arguments, 0 is returned.")
+#define FUNC_NAME s_scm_i_gcd
+{
+  while (!scm_is_null (rest))
+    { x = scm_gcd (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_gcd (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_gcd s_scm_i_gcd
+#define g_gcd g_scm_i_gcd
+
 SCM
 scm_gcd (SCM x, SCM y)
 {
   if (SCM_UNBNDP (y))
-    return SCM_UNBNDP (x) ? SCM_INUM0 : x;
+    return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x);
   
   if (SCM_I_INUMP (x))
     {
       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)))
                {
@@ -1073,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))
         {
@@ -1087,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)
@@ -1099,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))
         {
@@ -1117,10 +2391,24 @@ scm_gcd (SCM x, SCM y)
     SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG1, s_gcd);
 }
 
-SCM_GPROC1 (s_lcm, "lcm", scm_tc7_asubr, scm_lcm, g_lcm);
-/* "Return the least common multiple of the arguments.\n"
- * "If called without arguments, 1 is returned."
- */
+SCM_PRIMITIVE_GENERIC (scm_i_lcm, "lcm", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the least common multiple of the arguments.\n"
+                       "If called without arguments, 1 is returned.")
+#define FUNC_NAME s_scm_i_lcm
+{
+  while (!scm_is_null (rest))
+    { x = scm_lcm (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_lcm (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_lcm s_scm_i_lcm
+#define g_lcm g_scm_i_lcm
+
 SCM
 scm_lcm (SCM n1, SCM n2)
 {
@@ -1152,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);
@@ -1218,17 +2506,31 @@ scm_lcm (SCM n1, SCM n2)
 
 */
 
-SCM_DEFINE1 (scm_logand, "logand", scm_tc7_asubr,
-             (SCM n1, SCM n2),
-            "Return the bitwise AND of the integer arguments.\n\n"
-            "@lisp\n"
-            "(logand) @result{} -1\n"
-            "(logand 7) @result{} 7\n"
-            "(logand #b111 #b011 #b001) @result{} 1\n"
-            "@end lisp")
+SCM_DEFINE (scm_i_logand, "logand", 0, 2, 1,
+            (SCM x, SCM y, SCM rest),
+            "Return the bitwise AND of the integer arguments.\n\n"
+            "@lisp\n"
+            "(logand) @result{} -1\n"
+            "(logand 7) @result{} 7\n"
+            "(logand #b111 #b011 #b001) @result{} 1\n"
+            "@end lisp")
+#define FUNC_NAME s_scm_i_logand
+{
+  while (!scm_is_null (rest))
+    { x = scm_logand (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_logand (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_scm_logand s_scm_i_logand
+
+SCM scm_logand (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logand
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1247,7 +2549,7 @@ SCM_DEFINE1 (scm_logand, "logand", scm_tc7_asubr,
       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)
@@ -1294,17 +2596,31 @@ SCM_DEFINE1 (scm_logand, "logand", scm_tc7_asubr,
 #undef FUNC_NAME
 
 
-SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
-             (SCM n1, SCM n2),
-            "Return the bitwise OR of the integer arguments.\n\n"
-            "@lisp\n"
-            "(logior) @result{} 0\n"
-            "(logior 7) @result{} 7\n"
-            "(logior #b000 #b001 #b011) @result{} 3\n"
-           "@end lisp")
+SCM_DEFINE (scm_i_logior, "logior", 0, 2, 1,
+            (SCM x, SCM y, SCM rest),
+            "Return the bitwise OR of the integer arguments.\n\n"
+            "@lisp\n"
+            "(logior) @result{} 0\n"
+            "(logior 7) @result{} 7\n"
+            "(logior #b000 #b001 #b011) @result{} 3\n"
+            "@end lisp")
+#define FUNC_NAME s_scm_i_logior
+{
+  while (!scm_is_null (rest))
+    { x = scm_logior (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_logior (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_scm_logior s_scm_i_logior
+
+SCM scm_logior (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logior
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1368,8 +2684,8 @@ SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
 #undef FUNC_NAME
 
 
-SCM_DEFINE1 (scm_logxor, "logxor", scm_tc7_asubr,
-             (SCM n1, SCM n2),
+SCM_DEFINE (scm_i_logxor, "logxor", 0, 2, 1,
+            (SCM x, SCM y, SCM rest),
             "Return the bitwise XOR of the integer arguments.  A bit is\n"
             "set in the result if it is set in an odd number of arguments.\n"
             "@lisp\n"
@@ -1378,9 +2694,23 @@ SCM_DEFINE1 (scm_logxor, "logxor", scm_tc7_asubr,
             "(logxor #b000 #b001 #b011) @result{} 2\n"
             "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
            "@end lisp")
+#define FUNC_NAME s_scm_i_logxor
+{
+  while (!scm_is_null (rest))
+    { x = scm_logxor (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_logxor (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_scm_logxor s_scm_i_logxor
+
+SCM scm_logxor (SCM n1, SCM n2)
 #define FUNC_NAME s_scm_logxor
 {
-  long int nn1;
+  scm_t_inum nn1;
 
   if (SCM_UNBNDP (n2))
     {
@@ -1397,7 +2727,7 @@ SCM_DEFINE1 (scm_logxor, "logxor", scm_tc7_asubr,
       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))
@@ -1455,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))
@@ -1704,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"
@@ -1716,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;
 
@@ -1811,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)
         {
@@ -1826,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))
             {
@@ -1834,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;
@@ -1844,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));
         }
@@ -1907,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". */
@@ -1919,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;
@@ -1978,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)
@@ -2026,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)
@@ -2100,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)
@@ -2135,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");
@@ -2143,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;
@@ -2295,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';
@@ -2342,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++;
 
@@ -2352,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;
 }
@@ -2439,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 (scm_i_string_chars (str), scm_i_string_length (str), port);
+  scm_display (str, port);
   scm_remember_upto_here_1 (str);
   return !0;
 }
@@ -2485,14 +3831,29 @@ 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)                                  \
-  (isdigit ((int) (unsigned char) d)                    \
-   ? (d) - '0'                                          \
-   : 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 (const char* mem, size_t len, unsigned int *p_idx,
+mem2uinteger (SCM mem, unsigned int *p_idx,
              unsigned int radix, enum t_exactness *p_exactness)
 {
   unsigned int idx = *p_idx;
@@ -2502,14 +3863,13 @@ mem2uinteger (const char* mem, size_t len, unsigned int *p_idx,
   unsigned int digit_value;
   SCM result;
   char c;
+  size_t len = scm_i_string_length (mem);
 
   if (idx == len)
     return SCM_BOOL_F;
 
-  c = mem[idx];
-  if (!isxdigit ((int) (unsigned char) c))
-    return SCM_BOOL_F;
-  digit_value = XDIGIT2UINT (c);
+  c = scm_i_string_ref (mem, idx);
+  digit_value = char_decimal_value (c);
   if (digit_value >= radix)
     return SCM_BOOL_F;
 
@@ -2517,22 +3877,22 @@ mem2uinteger (const char* mem, size_t len, unsigned int *p_idx,
   result = SCM_I_MAKINUM (digit_value);
   while (idx != len)
     {
-      char c = mem[idx];
-      if (isxdigit ((int) (unsigned char) c))
-       {
-         if (hash_seen)
-           break;
-         digit_value = XDIGIT2UINT (c);
-         if (digit_value >= radix)
-           break;
-       }
-      else if (c == '#')
+      scm_t_wchar c = scm_i_string_ref (mem, idx);
+      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)
@@ -2571,31 +3931,31 @@ mem2uinteger (const char* mem, size_t len, unsigned int *p_idx,
  * has already been seen in the digits before the point.
  */
 
-/* In non ASCII-style encodings the following macro might not work. */
-#define DIGIT2UINT(d) ((d) - '0')
+#define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
 
 static SCM
-mem2decimal_from_point (SCM result, const char* mem, size_t len
+mem2decimal_from_point (SCM result, SCM mem
                        unsigned int *p_idx, enum t_exactness *p_exactness)
 {
   unsigned int idx = *p_idx;
   enum t_exactness x = *p_exactness;
+  size_t len = scm_i_string_length (mem);
 
   if (idx == len)
     return result;
 
-  if (mem[idx] == '.')
+  if (scm_i_string_ref (mem, idx) == '.')
     {
       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)
        {
-         char c = mem[idx];
-         if (isdigit ((int) (unsigned char) c))
+         scm_t_wchar c = scm_i_string_ref (mem, idx);
+         if (uc_is_property_decimal_digit ((scm_t_uint32) c))
            {
              if (x == INEXACT)
                return SCM_BOOL_F;
@@ -2645,13 +4005,13 @@ mem2decimal_from_point (SCM result, const char* mem, size_t len,
     {
       int sign = 1;
       unsigned int start;
-      char c;
+      scm_t_wchar c;
       int exponent;
       SCM e;
 
       /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
 
-      switch (mem[idx])
+      switch (scm_i_string_ref (mem, idx))
        {
        case 'd': case 'D':
        case 'e': case 'E':
@@ -2659,32 +4019,41 @@ mem2decimal_from_point (SCM result, const char* mem, size_t len,
        case 'l': case 'L':
        case 's': case 'S':
          idx++;
+          if (idx == len)
+            return SCM_BOOL_F;
+
          start = idx;
-         c = mem[idx];
+         c = scm_i_string_ref (mem, idx);
          if (c == '-')
            {
              idx++;
+              if (idx == len)
+                return SCM_BOOL_F;
+
              sign = -1;
-             c = mem[idx];
+             c = scm_i_string_ref (mem, idx);
            }
          else if (c == '+')
            {
              idx++;
+              if (idx == len)
+                return SCM_BOOL_F;
+
              sign = 1;
-             c = mem[idx];
+             c = scm_i_string_ref (mem, idx);
            }
          else
            sign = 1;
 
-         if (!isdigit ((int) (unsigned char) c))
+         if (!uc_is_property_decimal_digit ((scm_t_uint32) c))
            return SCM_BOOL_F;
 
          idx++;
          exponent = DIGIT2UINT (c);
          while (idx != len)
            {
-             char c = mem[idx];
-             if (isdigit ((int) (unsigned char) c))
+             scm_t_wchar c = scm_i_string_ref (mem, idx);
+             if (uc_is_property_decimal_digit ((scm_t_uint32) c))
                {
                  idx++;
                  if (exponent <= SCM_MAXEXP)
@@ -2697,7 +4066,7 @@ mem2decimal_from_point (SCM result, const char* mem, size_t len,
          if (exponent > SCM_MAXEXP)
            {
              size_t exp_len = idx - start;
-             SCM exp_string = scm_from_locale_stringn (&mem[start], exp_len);
+             SCM exp_string = scm_i_substring_copy (mem, start, start + exp_len);
              SCM exp_num = scm_string_to_number (exp_string, SCM_UNDEFINED);
              scm_out_of_range ("string->number", exp_num);
            }
@@ -2729,63 +4098,67 @@ mem2decimal_from_point (SCM result, const char* mem, size_t len,
 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
 
 static SCM
-mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
+mem2ureal (SCM mem, unsigned int *p_idx,
           unsigned int radix, enum t_exactness *p_exactness)
 {
   unsigned int idx = *p_idx;
   SCM result;
+  size_t len = scm_i_string_length (mem);
+
+  /* Start off believing that the number will be exact.  This changes
+     to INEXACT if we see a decimal point or a hash. */
+  enum t_exactness x = EXACT;
 
   if (idx == len)
     return SCM_BOOL_F;
 
-  if (idx+5 <= len && !strncmp (mem+idx, "inf.0", 5))
+  if (idx+5 <= len && !scm_i_string_strcmp (mem, idx, "inf.0"))
     {
       *p_idx = idx+5;
       return scm_inf ();
     }
 
-  if (idx+4 < len && !strncmp (mem+idx, "nan.", 4))
+  if (idx+4 < len && !scm_i_string_strcmp (mem, idx, "nan."))
     {
-      enum t_exactness x = EXACT;
-
       /* Cobble up the fractional part.  We might want to set the
         NaN's mantissa from it. */
       idx += 4;
-      mem2uinteger (mem, len, &idx, 10, &x);
+      mem2uinteger (mem, &idx, 10, &x);
       *p_idx = idx;
       return scm_nan ();
     }
 
-  if (mem[idx] == '.')
+  if (scm_i_string_ref (mem, idx) == '.')
     {
       if (radix != 10)
        return SCM_BOOL_F;
       else if (idx + 1 == len)
        return SCM_BOOL_F;
-      else if (!isdigit ((int) (unsigned char) mem[idx + 1]))
+      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, len,
-                                        p_idx, p_exactness);
+       result = mem2decimal_from_point (SCM_INUM0, mem,
+                                        p_idx, &x);
     }
   else
     {
-      enum t_exactness x = EXACT;
       SCM uinteger;
 
-      uinteger = mem2uinteger (mem, len, &idx, radix, &x);
+      uinteger = mem2uinteger (mem, &idx, radix, &x);
       if (scm_is_false (uinteger))
        return SCM_BOOL_F;
 
       if (idx == len)
        result = uinteger;
-      else if (mem[idx] == '/')
+      else if (scm_i_string_ref (mem, idx) == '/')
        {
          SCM divisor;
 
          idx++;
+          if (idx == len)
+            return SCM_BOOL_F;
 
-         divisor = mem2uinteger (mem, len, &idx, radix, &x);
+         divisor = mem2uinteger (mem, &idx, radix, &x);
          if (scm_is_false (divisor))
            return SCM_BOOL_F;
 
@@ -2794,7 +4167,7 @@ mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
        }
       else if (radix == 10)
        {
-         result = mem2decimal_from_point (uinteger, mem, len, &idx, &x);
+         result = mem2decimal_from_point (uinteger, mem, &idx, &x);
          if (scm_is_false (result))
            return SCM_BOOL_F;
        }
@@ -2802,14 +4175,20 @@ mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
        result = uinteger;
 
       *p_idx = idx;
-      if (x == INEXACT)
-       *p_exactness = x;
     }
 
+  /* Update *p_exactness if the number just read was inexact.  This is
+     important for complex numbers, so that a complex number is
+     treated as inexact overall if either its real or imaginary part
+     is inexact.
+  */
+  if (x == INEXACT)
+    *p_exactness = x;
+
   /* 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;
@@ -2819,17 +4198,18 @@ mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
 
 static SCM
-mem2complex (const char* mem, size_t len, unsigned int idx,
+mem2complex (SCM mem, unsigned int idx,
             unsigned int radix, enum t_exactness *p_exactness)
 {
-  char c;
+  scm_t_wchar c;
   int sign = 0;
   SCM ureal;
+  size_t len = scm_i_string_length (mem);
 
   if (idx == len)
     return SCM_BOOL_F;
 
-  c = mem[idx];
+  c = scm_i_string_ref (mem, idx);
   if (c == '+')
     {
       idx++;
@@ -2844,7 +4224,7 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
   if (idx == len)
     return SCM_BOOL_F;
 
-  ureal = mem2ureal (mem, len, &idx, radix, p_exactness);
+  ureal = mem2ureal (mem, &idx, radix, p_exactness);
   if (scm_is_false (ureal))
     {
       /* input must be either +i or -i */
@@ -2852,13 +4232,14 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
       if (sign == 0)
        return SCM_BOOL_F;
 
-      if (mem[idx] == 'i' || mem[idx] == 'I')
+      if (scm_i_string_ref (mem, idx) == 'i'
+         || scm_i_string_ref (mem, idx) == 'I')
        {
          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;
@@ -2871,7 +4252,7 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
       if (idx == len)
        return ureal;
 
-      c = mem[idx];
+      c = scm_i_string_ref (mem, idx);
       switch (c)
        {
        case 'i': case 'I':
@@ -2882,7 +4263,7 @@ mem2complex (const char* mem, size_t len, 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>. */
@@ -2896,21 +4277,25 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
              SCM angle;
              SCM result;
 
-             c = mem[idx];
+             c = scm_i_string_ref (mem, idx);
              if (c == '+')
                {
                  idx++;
+                  if (idx == len)
+                    return SCM_BOOL_F;
                  sign = 1;
                }
              else if (c == '-')
                {
                  idx++;
+                  if (idx == len)
+                    return SCM_BOOL_F;
                  sign = -1;
                }
              else
                sign = 1;
 
-             angle = mem2ureal (mem, len, &idx, radix, p_exactness);
+             angle = mem2ureal (mem, &idx, radix, p_exactness);
              if (scm_is_false (angle))
                return SCM_BOOL_F;
              if (idx != len)
@@ -2932,16 +4317,17 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
          else
            {
              int sign = (c == '+') ? 1 : -1;
-             SCM imag = mem2ureal (mem, len, &idx, radix, p_exactness);
+             SCM imag = mem2ureal (mem, &idx, radix, p_exactness);
 
              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)
                return SCM_BOOL_F;
-             if (mem[idx] != 'i' && mem[idx] != 'I')
+             if (scm_i_string_ref (mem, idx) != 'i'
+                 && scm_i_string_ref (mem, idx) != 'I')
                return SCM_BOOL_F;
 
              idx++;
@@ -2962,19 +4348,19 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
 enum t_radix {NO_RADIX=0, DUAL=2, OCT=8, DEC=10, HEX=16};
 
 SCM
-scm_c_locale_stringn_to_number (const char* mem, size_t len,
-                               unsigned int default_radix)
+scm_i_string_to_number (SCM mem, unsigned int default_radix)
 {
   unsigned int idx = 0;
   unsigned int radix = NO_RADIX;
   enum t_exactness forced_x = NO_EXACTNESS;
   enum t_exactness implicit_x = EXACT;
   SCM result;
+  size_t len = scm_i_string_length (mem);
 
   /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
-  while (idx + 2 < len && mem[idx] == '#')
+  while (idx + 2 < len && scm_i_string_ref (mem, idx) == '#')
     {
-      switch (mem[idx + 1])
+      switch (scm_i_string_ref (mem, idx + 1))
        {
        case 'b': case 'B':
          if (radix != NO_RADIX)
@@ -3014,9 +4400,9 @@ scm_c_locale_stringn_to_number (const char* mem, size_t len,
 
   /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
   if (radix == NO_RADIX)
-    result = mem2complex (mem, len, idx, default_radix, &implicit_x);
+    result = mem2complex (mem, idx, default_radix, &implicit_x);
   else
-    result = mem2complex (mem, len, idx, (unsigned int) radix, &implicit_x);
+    result = mem2complex (mem, idx, (unsigned int) radix, &implicit_x);
 
   if (scm_is_false (result))
     return SCM_BOOL_F;
@@ -3047,6 +4433,15 @@ scm_c_locale_stringn_to_number (const char* mem, size_t len,
     }
 }
 
+SCM
+scm_c_locale_stringn_to_number (const char* mem, size_t len,
+                               unsigned int default_radix)
+{
+  SCM str = scm_from_locale_stringn (mem, len);
+
+  return scm_i_string_to_number (str, default_radix);
+}
+
 
 SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
             (SCM string, SCM radix),
@@ -3069,9 +4464,7 @@ SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
   else
     base = scm_to_unsigned_integer (radix, 2, INT_MAX);
 
-  answer = scm_c_locale_stringn_to_number (scm_i_string_chars (string),
-                                          scm_i_string_length (string),
-                                          base);
+  answer = scm_i_string_to_number (string, base);
   scm_remember_upto_here_1 (string);
   return answer;
 }
@@ -3081,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"
@@ -3147,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
 
@@ -3160,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;
 }
@@ -3183,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))
@@ -3248,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))
@@ -3262,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))
     {
@@ -3277,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);
@@ -3288,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);
@@ -3297,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))
     {
@@ -3305,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);
@@ -3327,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))
     {
@@ -3347,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);
@@ -3365,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))
     {
@@ -3384,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;
@@ -3397,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;
@@ -3407,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);
 }
 
 
@@ -3420,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))
@@ -3453,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))
     {
@@ -3472,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);
@@ -3481,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))
     {
@@ -3490,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);
@@ -3501,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))
     {
@@ -3523,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;
@@ -3542,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
@@ -3587,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
@@ -3607,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));
@@ -3626,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);
@@ -3650,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);
@@ -3674,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
@@ -3684,9 +5096,23 @@ scm_negative_p (SCM x)
    unlike scm_less_p above which takes some trouble to preserve all bits in
    its test, such trouble is not required for min and max.  */
 
-SCM_GPROC1 (s_max, "max", scm_tc7_asubr, scm_max, g_max);
-/* "Return the maximum of all parameter values."
- */
+SCM_PRIMITIVE_GENERIC (scm_i_max, "max", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the maximum of all parameter values.")
+#define FUNC_NAME s_scm_i_max
+{
+  while (!scm_is_null (rest))
+    { x = scm_max (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_max (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_max s_scm_i_max
+#define g_max g_scm_i_max
+
 SCM
 scm_max (SCM x, SCM y)
 {
@@ -3702,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))
@@ -3778,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))
        {
@@ -3816,9 +5242,23 @@ scm_max (SCM x, SCM y)
 }
 
 
-SCM_GPROC1 (s_min, "min", scm_tc7_asubr, scm_min, g_min);
-/* "Return the minium of all parameter values."
- */
+SCM_PRIMITIVE_GENERIC (scm_i_min, "min", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the minimum of all parameter values.")
+#define FUNC_NAME s_scm_i_min
+{
+  while (!scm_is_null (rest))
+    { x = scm_min (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_min (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_min s_scm_i_min
+#define g_min g_scm_i_min
+
 SCM
 scm_min (SCM x, SCM y)
 {
@@ -3834,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))
@@ -3910,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))
        {
@@ -3941,17 +5381,31 @@ scm_min (SCM x, SCM y)
           goto use_less;
        }
       else
-       SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
+       SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
     }
   else
     SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARG1, s_min);
 }
 
 
-SCM_GPROC1 (s_sum, "+", scm_tc7_asubr, scm_sum, g_sum);
-/* "Return the sum of all parameter values.  Return 0 if called without\n"
- * "any parameters." 
- */
+SCM_PRIMITIVE_GENERIC (scm_i_sum, "+", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the sum of all parameter values.  Return 0 if called without\n"
+                       "any parameters." )
+#define FUNC_NAME s_scm_i_sum
+{
+  while (!scm_is_null (rest))
+    { x = scm_sum (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_sum (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_sum s_scm_i_sum
+#define g_sum g_scm_i_sum
+
 SCM
 scm_sum (SCM x, SCM y)
 {
@@ -3966,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))
         {
@@ -3978,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));
         }
@@ -3997,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);      
@@ -4136,18 +5590,33 @@ 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
 
 
-SCM_GPROC1 (s_difference, "-", scm_tc7_asubr, scm_difference, g_difference);
-/* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
- * the sum of all but the first argument are subtracted from the first
- * argument.  */
-#define FUNC_NAME s_difference
+SCM_PRIMITIVE_GENERIC (scm_i_difference, "-", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
+                       "the sum of all but the first argument are subtracted from the first\n"
+                       "argument.")
+#define FUNC_NAME s_scm_i_difference
+{
+  while (!scm_is_null (rest))
+    { x = scm_difference (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_difference (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_difference s_scm_i_difference
+#define g_difference g_scm_i_difference
+
 SCM
 scm_difference (SCM x, SCM y)
+#define FUNC_NAME s_difference
 {
   if (SCM_UNLIKELY (SCM_UNBNDP (y)))
     {
@@ -4156,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
@@ -4182,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));
@@ -4221,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));
        }
@@ -4243,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 ();
@@ -4381,15 +5854,29 @@ 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
 
 
-SCM_GPROC1 (s_product, "*", scm_tc7_asubr, scm_product, g_product);
-/* "Return the product of all arguments.  If called without arguments,\n"
- * "1 is returned."
- */
+SCM_PRIMITIVE_GENERIC (scm_i_product, "*", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the product of all arguments.  If called without arguments,\n"
+                       "1 is returned.")
+#define FUNC_NAME s_scm_i_product
+{
+  while (!scm_is_null (rest))
+    { x = scm_product (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_product (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_product s_scm_i_product
+#define g_product g_scm_i_product
+
 SCM
 scm_product (SCM x, SCM y)
 {
@@ -4405,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);
@@ -4414,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);
            }
@@ -4606,13 +6104,28 @@ arising out of or in connection with the use or performance of
 this software.
 ****************************************************************/
 
-SCM_GPROC1 (s_divide, "/", scm_tc7_asubr, scm_divide, g_divide);
-/* Divide the first argument by the product of the remaining
-   arguments.  If called with one argument @var{z1}, 1/@var{z1} is
-   returned.  */
-#define FUNC_NAME s_divide
+SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Divide the first argument by the product of the remaining\n"
+                       "arguments.  If called with one argument @var{z1}, 1/@var{z1} is\n"
+                       "returned.")
+#define FUNC_NAME s_scm_i_divide
+{
+  while (!scm_is_null (rest))
+    { x = scm_divide (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_divide (x, y);
+}
+#undef FUNC_NAME
+                       
+#define s_divide s_scm_i_divide
+#define g_divide g_scm_i_divide
+
 static SCM
-scm_i_divide (SCM x, SCM y, int inexact)
+do_divide (SCM x, SCM y, int inexact)
+#define FUNC_NAME s_divide
 {
   double a;
 
@@ -4622,7 +6135,7 @@ scm_i_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
@@ -4633,14 +6146,14 @@ scm_i_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))
        {
@@ -4678,10 +6191,10 @@ scm_i_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
@@ -4698,11 +6211,11 @@ scm_i_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))
@@ -4753,7 +6266,7 @@ scm_i_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
@@ -4776,7 +6289,7 @@ scm_i_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)
@@ -4798,47 +6311,33 @@ scm_i_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))
@@ -4867,7 +6366,7 @@ scm_i_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);
@@ -4907,7 +6406,7 @@ scm_i_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);
@@ -4963,7 +6462,7 @@ scm_i_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);
@@ -5005,61 +6504,16 @@ scm_i_divide (SCM x, SCM y, int inexact)
 SCM
 scm_divide (SCM x, SCM y)
 {
-  return scm_i_divide (x, y, 0);
+  return do_divide (x, y, 0);
 }
 
 static SCM scm_divide2real (SCM x, SCM y)
 {
-  return scm_i_divide (x, y, 1);
+  return do_divide (x, y, 1);
 }
 #undef FUNC_NAME
 
 
-double
-scm_asinh (double x)
-{
-#if HAVE_ASINH
-  return asinh (x);
-#else
-#define asinh scm_asinh
-  return log (x + sqrt (x * x + 1));
-#endif
-}
-SCM_GPROC1 (s_asinh, "$asinh", scm_tc7_dsubr, (SCM (*)()) asinh, g_asinh);
-/* "Return the inverse hyperbolic sine of @var{x}."
- */
-
-
-double
-scm_acosh (double x)
-{
-#if HAVE_ACOSH
-  return acosh (x);
-#else
-#define acosh scm_acosh
-  return log (x + sqrt (x * x - 1));
-#endif
-}
-SCM_GPROC1 (s_acosh, "$acosh", scm_tc7_dsubr, (SCM (*)()) acosh, g_acosh);
-/* "Return the inverse hyperbolic cosine of @var{x}."
- */
-
-
-double
-scm_atanh (double x)
-{
-#if HAVE_ATANH
-  return atanh (x);
-#else
-#define atanh scm_atanh
-  return 0.5 * log ((1 + x) / (1 - x));
-#endif
-}
-SCM_GPROC1 (s_atanh, "$atanh", scm_tc7_dsubr, (SCM (*)()) atanh, g_atanh);
-/* "Return the inverse hyperbolic tangent of @var{x}."
- */
-
-
 double
 scm_c_truncate (double x)
 {
@@ -5126,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. "
@@ -5149,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;
     }
@@ -5179,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
@@ -5210,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
@@ -5218,108 +6670,309 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-SCM_GPROC1 (s_i_sqrt, "$sqrt", scm_tc7_dsubr, (SCM (*)()) sqrt, g_i_sqrt);
-/* "Return the square root of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_abs, "$abs", scm_tc7_dsubr, (SCM (*)()) fabs, g_i_abs);
-/* "Return the absolute value of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_exp, "$exp", scm_tc7_dsubr, (SCM (*)()) exp, g_i_exp);
-/* "Return the @var{x}th power of e."
- */
-SCM_GPROC1 (s_i_log, "$log", scm_tc7_dsubr, (SCM (*)()) log, g_i_log);
-/* "Return the natural logarithm of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_sin, "$sin", scm_tc7_dsubr, (SCM (*)()) sin, g_i_sin);
-/* "Return the sine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_cos, "$cos", scm_tc7_dsubr, (SCM (*)()) cos, g_i_cos);
-/* "Return the cosine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_tan, "$tan", scm_tc7_dsubr, (SCM (*)()) tan, g_i_tan);
-/* "Return the tangent of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_asin, "$asin", scm_tc7_dsubr, (SCM (*)()) asin, g_i_asin);
-/* "Return the arc sine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_acos, "$acos", scm_tc7_dsubr, (SCM (*)()) acos, g_i_acos);
-/* "Return the arc cosine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_atan, "$atan", scm_tc7_dsubr, (SCM (*)()) atan, g_i_atan);
-/* "Return the arc tangent of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_sinh, "$sinh", scm_tc7_dsubr, (SCM (*)()) sinh, g_i_sinh);
-/* "Return the hyperbolic sine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_cosh, "$cosh", scm_tc7_dsubr, (SCM (*)()) cosh, g_i_cosh);
-/* "Return the hyperbolic cosine of the real number @var{x}."
- */
-SCM_GPROC1 (s_i_tanh, "$tanh", scm_tc7_dsubr, (SCM (*)()) tanh, g_i_tanh);
-/* "Return the hyperbolic tangent of the real number @var{x}."
- */
+/* sin/cos/tan/asin/acos/atan
+   sinh/cosh/tanh/asinh/acosh/atanh
+   Derived from "Transcen.scm", Complex trancendental functions for SCM.
+   Written by Jerry D. Hedden, (C) FSF.
+   See the file `COPYING' for terms applying to this program. */
 
-struct dpair
+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
 {
-  double 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 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
 
-static void scm_two_doubles (SCM x,
-                            SCM y,
-                            const char *sstring,
-                            struct dpair * xy);
+SCM_PRIMITIVE_GENERIC (scm_sin, "sin", 1, 0, 0,
+                       (SCM z),
+                       "Compute the sine of @var{z}.")
+#define FUNC_NAME s_scm_sin
+{
+  if (scm_is_real (z))
+    return scm_from_double (sin (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_c_make_rectangular (sin (x) * cosh (y),
+                                     cos (x) * sinh (y));
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_sin, z, 1, s_scm_sin);
+}
+#undef FUNC_NAME
 
-static void
-scm_two_doubles (SCM x, SCM y, const char *sstring, struct dpair *xy)
+SCM_PRIMITIVE_GENERIC (scm_cos, "cos", 1, 0, 0,
+                       (SCM z),
+                       "Compute the cosine of @var{z}.")
+#define FUNC_NAME s_scm_cos
 {
-  if (SCM_I_INUMP (x))
-    xy->x = SCM_I_INUM (x);
-  else if (SCM_BIGP (x))
-    xy->x = scm_i_big2dbl (x);
-  else if (SCM_REALP (x))
-    xy->x = SCM_REAL_VALUE (x);
-  else if (SCM_FRACTIONP (x))
-    xy->x = scm_i_fraction2double (x);
+  if (scm_is_real (z))
+    return scm_from_double (cos (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_c_make_rectangular (cos (x) * cosh (y),
+                                     -sin (x) * sinh (y));
+    }
   else
-    scm_wrong_type_arg (sstring, SCM_ARG1, x);
-
-  if (SCM_I_INUMP (y))
-    xy->y = SCM_I_INUM (y);
-  else if (SCM_BIGP (y))
-    xy->y = scm_i_big2dbl (y);
-  else if (SCM_REALP (y))
-    xy->y = SCM_REAL_VALUE (y);
-  else if (SCM_FRACTIONP (y))
-    xy->y = scm_i_fraction2double (y);
+    SCM_WTA_DISPATCH_1 (g_scm_cos, z, 1, s_scm_cos);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_tan, "tan", 1, 0, 0,
+                       (SCM z),
+                       "Compute the tangent of @var{z}.")
+#define FUNC_NAME s_scm_tan
+{
+  if (scm_is_real (z))
+    return scm_from_double (tan (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y, w;
+      x = 2.0 * SCM_COMPLEX_REAL (z);
+      y = 2.0 * SCM_COMPLEX_IMAG (z);
+      w = cos (x) + cosh (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
+      if (w == 0.0)
+        scm_num_overflow (s_scm_tan);
+#endif
+      return scm_c_make_rectangular (sin (x) / w, sinh (y) / w);
+    }
   else
-    scm_wrong_type_arg (sstring, SCM_ARG2, y);
+    SCM_WTA_DISPATCH_1 (g_scm_tan, z, 1, s_scm_tan);
 }
+#undef FUNC_NAME
 
+SCM_PRIMITIVE_GENERIC (scm_sinh, "sinh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the hyperbolic sine of @var{z}.")
+#define FUNC_NAME s_scm_sinh
+{
+  if (scm_is_real (z))
+    return scm_from_double (sinh (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_c_make_rectangular (sinh (x) * cos (y),
+                                     cosh (x) * sin (y));
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_sinh, z, 1, s_scm_sinh);
+}
+#undef FUNC_NAME
 
-SCM_DEFINE (scm_sys_expt, "$expt", 2, 0, 0,
-            (SCM x, SCM y),
-           "Return @var{x} raised to the power of @var{y}. This\n"
-           "procedure does not accept complex arguments.") 
-#define FUNC_NAME s_scm_sys_expt
+SCM_PRIMITIVE_GENERIC (scm_cosh, "cosh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the hyperbolic cosine of @var{z}.")
+#define FUNC_NAME s_scm_cosh
+{
+  if (scm_is_real (z))
+    return scm_from_double (cosh (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_c_make_rectangular (cosh (x) * cos (y),
+                                     sinh (x) * sin (y));
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_cosh, z, 1, s_scm_cosh);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_tanh, "tanh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the hyperbolic tangent of @var{z}.")
+#define FUNC_NAME s_scm_tanh
+{
+  if (scm_is_real (z))
+    return scm_from_double (tanh (scm_to_double (z)));
+  else if (SCM_COMPLEXP (z))
+    { double x, y, w;
+      x = 2.0 * SCM_COMPLEX_REAL (z);
+      y = 2.0 * SCM_COMPLEX_IMAG (z);
+      w = cosh (x) + cos (y);
+#ifndef ALLOW_DIVIDE_BY_ZERO
+      if (w == 0.0)
+        scm_num_overflow (s_scm_tanh);
+#endif
+      return scm_c_make_rectangular (sinh (x) / w, sin (y) / w);
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_tanh, z, 1, s_scm_tanh);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_asin, "asin", 1, 0, 0,
+                       (SCM z),
+                       "Compute the arc sine of @var{z}.")
+#define FUNC_NAME s_scm_asin
 {
-  struct dpair xy;
-  scm_two_doubles (x, y, FUNC_NAME, &xy);
-  return scm_from_double (pow (xy.x, xy.y));
+  if (scm_is_real (z))
+    {
+      double w = scm_to_double (z);
+      if (w >= -1.0 && w <= 1.0)
+        return scm_from_double (asin (w));
+      else
+        return scm_product (scm_c_make_rectangular (0, -1),
+                            scm_sys_asinh (scm_c_make_rectangular (0, w)));
+    }
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_product (scm_c_make_rectangular (0, -1),
+                          scm_sys_asinh (scm_c_make_rectangular (-y, x)));
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_asin, z, 1, s_scm_asin);
 }
 #undef FUNC_NAME
 
+SCM_PRIMITIVE_GENERIC (scm_acos, "acos", 1, 0, 0,
+                       (SCM z),
+                       "Compute the arc cosine of @var{z}.")
+#define FUNC_NAME s_scm_acos
+{
+  if (scm_is_real (z))
+    {
+      double w = scm_to_double (z);
+      if (w >= -1.0 && w <= 1.0)
+        return scm_from_double (acos (w));
+      else
+        return scm_sum (scm_from_double (acos (0.0)),
+                        scm_product (scm_c_make_rectangular (0, 1),
+                                     scm_sys_asinh (scm_c_make_rectangular (0, w))));
+    }
+  else if (SCM_COMPLEXP (z))
+    { double x, y;
+      x = SCM_COMPLEX_REAL (z);
+      y = SCM_COMPLEX_IMAG (z);
+      return scm_sum (scm_from_double (acos (0.0)),
+                      scm_product (scm_c_make_rectangular (0, 1),
+                                   scm_sys_asinh (scm_c_make_rectangular (-y, x))));
+    }
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_acos, z, 1, s_scm_acos);
+}
+#undef FUNC_NAME
 
-SCM_DEFINE (scm_sys_atan2, "$atan2", 2, 0, 0,
-            (SCM x, SCM y),
-           "Return the arc tangent of the two arguments @var{x} and\n"
-           "@var{y}. This is similar to calculating the arc tangent of\n"
-           "@var{x} / @var{y}, except that the signs of both arguments\n"
-           "are used to determine the quadrant of the result. This\n"
-           "procedure does not accept complex arguments.")
-#define FUNC_NAME s_scm_sys_atan2
+SCM_PRIMITIVE_GENERIC (scm_atan, "atan", 1, 1, 0,
+                       (SCM z, SCM y),
+                       "With one argument, compute the arc tangent of @var{z}.\n"
+                       "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
+                       "using the sign of @var{z} and @var{y} to determine the quadrant.")
+#define FUNC_NAME s_scm_atan
 {
-  struct dpair xy;
-  scm_two_doubles (x, y, FUNC_NAME, &xy);
-  return scm_from_double (atan2 (xy.x, xy.y));
+  if (SCM_UNBNDP (y))
+    {
+      if (scm_is_real (z))
+        return scm_from_double (atan (scm_to_double (z)));
+      else if (SCM_COMPLEXP (z))
+        {
+          double v, w;
+          v = SCM_COMPLEX_REAL (z);
+          w = SCM_COMPLEX_IMAG (z);
+          return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v, w - 1.0),
+                                                  scm_c_make_rectangular (v, w + 1.0))),
+                             scm_c_make_rectangular (0, 2));
+        }
+      else
+        SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan);
+    }
+  else if (scm_is_real (z))
+    {
+      if (scm_is_real (y))
+        return scm_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
+      else
+        SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan);
+    }
+  else
+    SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_sys_asinh, "asinh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the inverse hyperbolic sine of @var{z}.")
+#define FUNC_NAME s_scm_sys_asinh
+{
+  if (scm_is_real (z))
+    return scm_from_double (asinh (scm_to_double (z)));
+  else if (scm_is_number (z))
+    return scm_log (scm_sum (z,
+                             scm_sqrt (scm_sum (scm_product (z, z),
+                                                SCM_INUM1))));
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_sys_acosh, "acosh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the inverse hyperbolic cosine of @var{z}.")
+#define FUNC_NAME s_scm_sys_acosh
+{
+  if (scm_is_real (z) && scm_to_double (z) >= 1.0)
+    return scm_from_double (acosh (scm_to_double (z)));
+  else if (scm_is_number (z))
+    return scm_log (scm_sum (z,
+                             scm_sqrt (scm_difference (scm_product (z, z),
+                                                       SCM_INUM1))));
+  else
+    SCM_WTA_DISPATCH_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh);
+}
+#undef FUNC_NAME
+
+SCM_PRIMITIVE_GENERIC (scm_sys_atanh, "atanh", 1, 0, 0,
+                       (SCM z),
+                       "Compute the inverse hyperbolic tangent of @var{z}.")
+#define FUNC_NAME s_scm_sys_atanh
+{
+  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_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);
 }
 #undef FUNC_NAME
 
@@ -5331,8 +6984,10 @@ scm_c_make_rectangular (double re, double im)
   else
     {
       SCM z;
-      SCM_NEWSMOB (z, scm_tc16_complex, scm_gc_malloc (sizeof (scm_t_complex),
-                                                      "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;
@@ -5345,9 +7000,12 @@ SCM_DEFINE (scm_make_rectangular, "make-rectangular", 2, 0, 0,
            "and @var{imaginary-part} parts.")
 #define FUNC_NAME s_scm_make_rectangular
 {
-  struct dpair xy;
-  scm_two_doubles (real_part, imaginary_part, FUNC_NAME, &xy);
-  return scm_c_make_rectangular (xy.x, xy.y);
+  SCM_ASSERT_TYPE (scm_is_real (real_part), real_part,
+                   SCM_ARG1, FUNC_NAME, "real");
+  SCM_ASSERT_TYPE (scm_is_real (imaginary_part), imaginary_part,
+                   SCM_ARG2, FUNC_NAME, "real");
+  return scm_c_make_rectangular (scm_to_double (real_part),
+                                 scm_to_double (imaginary_part));
 }
 #undef FUNC_NAME
 
@@ -5355,7 +7013,12 @@ SCM
 scm_c_make_polar (double mag, double ang)
 {
   double s, c;
-#if HAVE_SINCOS
+
+  /* The sincos(3) function is undocumented an broken on Tru64.  Thus we only
+     use it on Glibc-based systems that have it (it's a GNU extension).  See
+     http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
+     details.  */
+#if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
   sincos (ang, &s, &c);
 #else
   s = sin (ang);
@@ -5369,107 +7032,93 @@ SCM_DEFINE (scm_make_polar, "make-polar", 2, 0, 0,
            "Return the complex number @var{x} * e^(i * @var{y}).")
 #define FUNC_NAME s_scm_make_polar
 {
-  struct dpair xy;
-  scm_two_doubles (x, y, FUNC_NAME, &xy);
-  return scm_c_make_polar (xy.x, xy.y);
+  SCM_ASSERT_TYPE (scm_is_real (x), x, SCM_ARG1, FUNC_NAME, "real");
+  SCM_ASSERT_TYPE (scm_is_real (y), y, SCM_ARG2, FUNC_NAME, "real");
+  return scm_c_make_polar (scm_to_double (x), scm_to_double (y));
 }
 #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))
     {
@@ -5492,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));
     }
@@ -5520,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));
     }
@@ -5534,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));
@@ -5557,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
        {
@@ -5594,13 +7242,23 @@ 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
 
 SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0, 
-            (SCM x, SCM err),
-           "Return an exact number that is within @var{err} of @var{x}.")
+            (SCM x, SCM eps),
+           "Returns the @emph{simplest} rational number differing\n"
+           "from @var{x} by no more than @var{eps}.\n"
+           "\n"
+           "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
+           "exact result when both its arguments are exact.  Thus, you might need\n"
+           "to use @code{inexact->exact} on the arguments.\n"
+           "\n"
+           "@lisp\n"
+           "(rationalize (inexact->exact 1.2) 1/100)\n"
+           "@result{} 6/5\n"
+           "@end lisp")
 #define FUNC_NAME s_scm_rationalize
 {
   if (SCM_I_INUMP (x))
@@ -5615,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;
 
@@ -5632,7 +7290,7 @@ SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
         converges after less than a dozen iterations.
       */
 
-      err = scm_abs (err);
+      eps = scm_abs (eps);
       while (++i < 1000000)
        {
          a = scm_sum (scm_product (a1, tt), a2);    /* a = a1*tt + a2 */
@@ -5640,11 +7298,11 @@ SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
          if (scm_is_false (scm_zero_p (b)) &&         /* b != 0 */
              scm_is_false 
              (scm_gr_p (scm_abs (scm_difference (ex, scm_divide (a, b))),
-                        err)))                      /* abs(x-a/b) <= err */
+                        eps)))                      /* abs(x-a/b) <= eps */
            {
              SCM res = scm_sum (int_part, scm_divide (a, b));
              if (scm_is_false (scm_exact_p (x))
-                 || scm_is_false (scm_exact_p (err)))
+                 || scm_is_false (scm_exact_p (eps)))
                return scm_exact_to_inexact (res);
              else
                return res;
@@ -5843,7 +7501,13 @@ scm_i_range_error (SCM bad_val, SCM min, SCM max)
 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
 #include "libguile/conv-uinteger.i.c"
 
-#if SCM_HAVE_T_INT64
+#define TYPE                     scm_t_wchar
+#define TYPE_MIN                 (scm_t_int32)-1
+#define TYPE_MAX                 (scm_t_int32)0x10ffff
+#define SIZEOF_TYPE              4
+#define SCM_TO_TYPE_PROTO(arg)   scm_to_wchar (arg)
+#define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
+#include "libguile/conv-integer.i.c"
 
 #define TYPE                     scm_t_int64
 #define TYPE_MIN                 SCM_T_INT64_MIN
@@ -5861,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)
 {
@@ -5910,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);
@@ -5933,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);
@@ -6009,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))
@@ -6025,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.  */
@@ -6036,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))
@@ -6050,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);
@@ -6059,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.  */
@@ -6070,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))
@@ -6089,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
 
@@ -6156,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)
@@ -6166,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"
 }