Merge remote-tracking branch 'origin/stable-2.0'

Conflicts:
	libguile/numbers.c
	libguile/vm-i-scheme.c

20 files changed, 1495 insertions, 333 deletions

diff --git a/lib/Makefile.am b/lib/Makefile.am
index 8857a90ce..2ba04b72b 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -21,7 +21,7 @@
 # the same distribution terms as the rest of that program.
 #
 # Generated by gnulib-tool.
-# Reproduce by: gnulib-tool --import --dir=. --local-dir=gnulib-local --lib=libgnu --source-base=lib --m4-base=m4 --doc-base=doc --tests-base=tests --aux-dir=build-aux --lgpl=3 --no-conditional-dependencies --libtool --macro-prefix=gl --no-vc-files accept alignof alloca-opt announce-gen autobuild bind byteswap c-strcase canonicalize-lgpl ceil clock-time close connect dirfd duplocale environ extensions flock floor fpieee frexp fstat full-read full-write func gendocs getaddrinfo getlogin getpeername getsockname getsockopt git-version-gen gitlog-to-changelog gnu-web-doc-update gnupload havelib iconv_open-utf inet_ntop inet_pton isinf isnan ldexp lib-symbol-versions lib-symbol-visibility libunistring listen localcharset locale log1p maintainer-makefile malloc-gnu malloca nl_langinfo nproc open pipe-posix pipe2 poll putenv recv recvfrom regex rename select send sendto setenv setsockopt shutdown socket stat-time stdlib strftime striconveh string sys_stat time times trunc verify vsnprintf warnings wchar
+# Reproduce by: gnulib-tool --import --dir=. --local-dir=gnulib-local --lib=libgnu --source-base=lib --m4-base=m4 --doc-base=doc --tests-base=tests --aux-dir=build-aux --lgpl=3 --no-conditional-dependencies --libtool --macro-prefix=gl --no-vc-files accept alignof alloca-opt announce-gen autobuild bind byteswap c-strcase canonicalize-lgpl ceil clock-time close connect copysign dirfd duplocale environ extensions flock floor fpieee frexp fstat full-read full-write func gendocs getaddrinfo getlogin getpeername getsockname getsockopt git-version-gen gitlog-to-changelog gnu-web-doc-update gnupload havelib iconv_open-utf inet_ntop inet_pton isfinite isinf isnan ldexp lib-symbol-versions lib-symbol-visibility libunistring listen localcharset locale log1p maintainer-makefile malloc-gnu malloca nl_langinfo nproc open pipe-posix pipe2 poll putenv recv recvfrom regex rename select send sendto setenv setsockopt shutdown socket stat-time stdlib strftime striconveh string sys_stat time times trunc verify vsnprintf warnings wchar
 
 AUTOMAKE_OPTIONS = 1.5 gnits subdir-objects
 
@@ -50,6 +50,7 @@ EXTRA_libgnu_la_SOURCES =
 libgnu_la_LDFLAGS = $(AM_LDFLAGS)
 libgnu_la_LDFLAGS += -no-undefined
 libgnu_la_LDFLAGS += $(CEIL_LIBM)
+libgnu_la_LDFLAGS += $(COPYSIGN_LIBM)
 libgnu_la_LDFLAGS += $(FLOOR_LIBM)
 libgnu_la_LDFLAGS += $(FREXP_LIBM)
 libgnu_la_LDFLAGS += $(GETADDRINFO_LIB)
@@ -312,6 +313,15 @@ EXTRA_libgnu_la_SOURCES += connect.c
 
 ## end   gnulib module connect
 
+## begin gnulib module copysign
+
+
+EXTRA_DIST += copysign.c
+
+EXTRA_libgnu_la_SOURCES += copysign.c
+
+## end   gnulib module copysign
+
 ## begin gnulib module dirent
 
 BUILT_SOURCES += dirent.h
@@ -753,6 +763,15 @@ EXTRA_libgnu_la_SOURCES += inet_pton.c
 
 ## end   gnulib module inet_pton
 
+## begin gnulib module isfinite
+
+
+EXTRA_DIST += isfinite.c
+
+EXTRA_libgnu_la_SOURCES += isfinite.c
+
+## end   gnulib module isfinite
+
 ## begin gnulib module isinf
 
 
@@ -789,6 +808,15 @@ EXTRA_libgnu_la_SOURCES += isnan.c isnanf.c
 
 ## end   gnulib module isnanf
 
+## begin gnulib module isnanf-nolibm
+
+
+EXTRA_DIST += float+.h isnan.c isnanf-nolibm.h isnanf.c
+
+EXTRA_libgnu_la_SOURCES += isnan.c isnanf.c
+
+## end   gnulib module isnanf-nolibm
+
 ## begin gnulib module isnanl
 
 
@@ -798,6 +826,15 @@ EXTRA_libgnu_la_SOURCES += isnan.c isnanl.c
 
 ## end   gnulib module isnanl
 
+## begin gnulib module isnanl-nolibm
+
+
+EXTRA_DIST += float+.h isnan.c isnanl-nolibm.h isnanl.c
+
+EXTRA_libgnu_la_SOURCES += isnan.c isnanl.c
+
+## end   gnulib module isnanl-nolibm
+
 ## begin gnulib module langinfo
 
 BUILT_SOURCES += langinfo.h
@@ -1734,6 +1771,15 @@ EXTRA_DIST += signal.in.h
 
 ## end   gnulib module signal-h
 
+## begin gnulib module signbit
+
+
+EXTRA_DIST += float+.h signbitd.c signbitf.c signbitl.c
+
+EXTRA_libgnu_la_SOURCES += signbitd.c signbitf.c signbitl.c
+
+## end   gnulib module signbit
+
 ## begin gnulib module size_max
 
 libgnu_la_SOURCES += size_max.h
diff --git a/lib/copysign.c b/lib/copysign.c
new file mode 100644
index 000000000..61efeb8de
--- /dev/null
+++ b/lib/copysign.c
@@ -0,0 +1,26 @@
+/* Copy sign into another 'double' number.
+   Copyright (C) 2011-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+double
+copysign (double x, double y)
+{
+  return (signbit (x) != signbit (y) ? - x : x);
+}
diff --git a/lib/isfinite.c b/lib/isfinite.c
new file mode 100644
index 000000000..d9eddf54b
--- /dev/null
+++ b/lib/isfinite.c
@@ -0,0 +1,51 @@
+/* Test for finite value (zero, subnormal, or normal, and not infinite or NaN).
+   Copyright (C) 2007-2013 Free Software Foundation, Inc.
+
+   This program 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, or (at your option)
+   any later version.
+
+   This program 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 program; if not, see <http://www.gnu.org/licenses/>.  */
+
+/* Written by Ben Pfaff <blp@gnu.org>, 2007. */
+
+#include <config.h>
+
+#include "isnanf-nolibm.h"
+#include "isnand-nolibm.h"
+#include "isnanl-nolibm.h"
+
+/* The "cc" compiler on HP-UX 11.11, when optimizing, simplifies the test
+   x - y == 0.0  to  x == y, a simplification which is invalid when x and y
+   are Infinity.  Disable this optimization.  */
+#if defined __hpux && !defined __GNUC__
+static float zerof;
+static double zerod;
+static long double zerol;
+#else
+# define zerof 0.f
+# define zerod 0.
+# define zerol 0.L
+#endif
+
+int gl_isfinitef (float x)
+{
+  return !isnanf (x) && x - x == zerof;
+}
+
+int gl_isfinited (double x)
+{
+  return !isnand (x) && x - x == zerod;
+}
+
+int gl_isfinitel (long double x)
+{
+  return !isnanl (x) && x - x == zerol;
+}
diff --git a/lib/isnanf-nolibm.h b/lib/isnanf-nolibm.h
new file mode 100644
index 000000000..56f4fde61
--- /dev/null
+++ b/lib/isnanf-nolibm.h
@@ -0,0 +1,40 @@
+/* Test for NaN that does not need libm.
+   Copyright (C) 2007-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#if HAVE_ISNANF_IN_LIBC
+/* Get declaration of isnan macro or (older) isnanf function.  */
+# include <math.h>
+# if __GNUC__ >= 4
+   /* GCC 4.0 and newer provides three built-ins for isnan.  */
+#  undef isnanf
+#  define isnanf(x) __builtin_isnanf ((float)(x))
+# elif defined isnan
+#  undef isnanf
+#  define isnanf(x) isnan ((float)(x))
+# else
+   /* Get declaration of isnanf(), if not declared in <math.h>.  */
+#  if defined __sgi
+   /* We can't include <ieeefp.h>, because it conflicts with our definition of
+      isnand.  Therefore declare isnanf separately.  */
+extern int isnanf (float x);
+#  endif
+# endif
+#else
+/* Test whether X is a NaN.  */
+# undef isnanf
+# define isnanf rpl_isnanf
+extern int isnanf (float x);
+#endif
diff --git a/lib/isnanl-nolibm.h b/lib/isnanl-nolibm.h
new file mode 100644
index 000000000..c5d0323b4
--- /dev/null
+++ b/lib/isnanl-nolibm.h
@@ -0,0 +1,33 @@
+/* Test for NaN that does not need libm.
+   Copyright (C) 2007-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#if HAVE_ISNANL_IN_LIBC
+/* Get declaration of isnan macro or (older) isnanl function.  */
+# include <math.h>
+# if __GNUC__ >= 4
+   /* GCC 4.0 and newer provides three built-ins for isnan.  */
+#  undef isnanl
+#  define isnanl(x) __builtin_isnanl ((long double)(x))
+# elif defined isnan
+#  undef isnanl
+#  define isnanl(x) isnan ((long double)(x))
+# endif
+#else
+/* Test whether X is a NaN.  */
+# undef isnanl
+# define isnanl rpl_isnanl
+extern int isnanl (long double x);
+#endif
diff --git a/lib/signbitd.c b/lib/signbitd.c
new file mode 100644
index 000000000..1c813dabf
--- /dev/null
+++ b/lib/signbitd.c
@@ -0,0 +1,64 @@
+/* signbit() macro: Determine the sign bit of a floating-point number.
+   Copyright (C) 2007-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+#include <string.h>
+#include "isnand-nolibm.h"
+#include "float+.h"
+
+#ifdef gl_signbitd_OPTIMIZED_MACRO
+# undef gl_signbitd
+#endif
+
+int
+gl_signbitd (double arg)
+{
+#if defined DBL_SIGNBIT_WORD && defined DBL_SIGNBIT_BIT
+  /* The use of a union to extract the bits of the representation of a
+     'long double' is safe in practice, despite of the "aliasing rules" of
+     C99, because the GCC docs say
+       "Even with '-fstrict-aliasing', type-punning is allowed, provided the
+        memory is accessed through the union type."
+     and similarly for other compilers.  */
+# define NWORDS \
+    ((sizeof (double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
+  union { double value; unsigned int word[NWORDS]; } m;
+  m.value = arg;
+  return (m.word[DBL_SIGNBIT_WORD] >> DBL_SIGNBIT_BIT) & 1;
+#elif HAVE_COPYSIGN_IN_LIBC
+  return copysign (1.0, arg) < 0;
+#else
+  /* This does not do the right thing for NaN, but this is irrelevant for
+     most use cases.  */
+  if (isnand (arg))
+    return 0;
+  if (arg < 0.0)
+    return 1;
+  else if (arg == 0.0)
+    {
+      /* Distinguish 0.0 and -0.0.  */
+      static double plus_zero = 0.0;
+      double arg_mem = arg;
+      return (memcmp (&plus_zero, &arg_mem, SIZEOF_DBL) != 0);
+    }
+  else
+    return 0;
+#endif
+}
diff --git a/lib/signbitf.c b/lib/signbitf.c
new file mode 100644
index 000000000..817484b8f
--- /dev/null
+++ b/lib/signbitf.c
@@ -0,0 +1,64 @@
+/* signbit() macro: Determine the sign bit of a floating-point number.
+   Copyright (C) 2007, 2009-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+#include <string.h>
+#include "isnanf-nolibm.h"
+#include "float+.h"
+
+#ifdef gl_signbitf_OPTIMIZED_MACRO
+# undef gl_signbitf
+#endif
+
+int
+gl_signbitf (float arg)
+{
+#if defined FLT_SIGNBIT_WORD && defined FLT_SIGNBIT_BIT
+  /* The use of a union to extract the bits of the representation of a
+     'long double' is safe in practice, despite of the "aliasing rules" of
+     C99, because the GCC docs say
+       "Even with '-fstrict-aliasing', type-punning is allowed, provided the
+        memory is accessed through the union type."
+     and similarly for other compilers.  */
+# define NWORDS \
+    ((sizeof (float) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
+  union { float value; unsigned int word[NWORDS]; } m;
+  m.value = arg;
+  return (m.word[FLT_SIGNBIT_WORD] >> FLT_SIGNBIT_BIT) & 1;
+#elif HAVE_COPYSIGNF_IN_LIBC
+  return copysignf (1.0f, arg) < 0;
+#else
+  /* This does not do the right thing for NaN, but this is irrelevant for
+     most use cases.  */
+  if (isnanf (arg))
+    return 0;
+  if (arg < 0.0f)
+    return 1;
+  else if (arg == 0.0f)
+    {
+      /* Distinguish 0.0f and -0.0f.  */
+      static float plus_zero = 0.0f;
+      float arg_mem = arg;
+      return (memcmp (&plus_zero, &arg_mem, SIZEOF_FLT) != 0);
+    }
+  else
+    return 0;
+#endif
+}
diff --git a/lib/signbitl.c b/lib/signbitl.c
new file mode 100644
index 000000000..159cfceb8
--- /dev/null
+++ b/lib/signbitl.c
@@ -0,0 +1,64 @@
+/* signbit() macro: Determine the sign bit of a floating-point number.
+   Copyright (C) 2007, 2009-2013 Free Software Foundation, Inc.
+
+   This program 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 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+#include <string.h>
+#include "isnanl-nolibm.h"
+#include "float+.h"
+
+#ifdef gl_signbitl_OPTIMIZED_MACRO
+# undef gl_signbitl
+#endif
+
+int
+gl_signbitl (long double arg)
+{
+#if defined LDBL_SIGNBIT_WORD && defined LDBL_SIGNBIT_BIT
+  /* The use of a union to extract the bits of the representation of a
+     'long double' is safe in practice, despite of the "aliasing rules" of
+     C99, because the GCC docs say
+       "Even with '-fstrict-aliasing', type-punning is allowed, provided the
+        memory is accessed through the union type."
+     and similarly for other compilers.  */
+# define NWORDS \
+    ((sizeof (long double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
+  union { long double value; unsigned int word[NWORDS]; } m;
+  m.value = arg;
+  return (m.word[LDBL_SIGNBIT_WORD] >> LDBL_SIGNBIT_BIT) & 1;
+#elif HAVE_COPYSIGNL_IN_LIBC
+  return copysignl (1.0L, arg) < 0;
+#else
+  /* This does not do the right thing for NaN, but this is irrelevant for
+     most use cases.  */
+  if (isnanl (arg))
+    return 0;
+  if (arg < 0.0L)
+    return 1;
+  else if (arg == 0.0L)
+    {
+      /* Distinguish 0.0L and -0.0L.  */
+      static long double plus_zero = 0.0L;
+      long double arg_mem = arg;
+      return (memcmp (&plus_zero, &arg_mem, SIZEOF_LDBL) != 0);
+    }
+  else
+    return 0;
+#endif
+}
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 3c0d76505..f549193b5 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -91,15 +91,6 @@ verify (FLT_RADIX == 2);
 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))
-
-/* On some platforms, isinf(x) returns 0, 1 or -1, indicating the sign
-   of the infinity, but other platforms return a boolean only. */
-#define DOUBLE_IS_POSITIVE_INFINITY(x) (isinf(x) && ((x) > 0))
-#define DOUBLE_IS_NEGATIVE_INFINITY(x) (isinf(x) && ((x) < 0))
-
 /* Test an inum to see if it can be converted to a double without loss
    of precision.  Note that this will sometimes return 0 even when 1
    could have been returned, e.g. for large powers of 2.  It is designed
@@ -654,12 +645,17 @@ scm_i_fraction2double (SCM z)
                               SCM_FRACTION_DENOMINATOR (z));
 }
 
-static int
-double_is_non_negative_zero (double x)
+static SCM
+scm_i_from_double (double val)
 {
-  static double zero = 0.0;
+  SCM z;
+
+  z = SCM_PACK_POINTER (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real"));
+
+  SCM_SET_CELL_TYPE (z, scm_tc16_real);
+  SCM_REAL_VALUE (z) = val;
 
-  return !memcmp (&x, &zero, sizeof(double));
+  return z;
 }
 
 SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0, 
@@ -724,7 +720,7 @@ SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0,
   else if (SCM_REALP (n))
     {
       double val = SCM_REAL_VALUE (n);
-      if (DOUBLE_IS_FINITE (val))
+      if (isfinite (val))
 	{
 	  double rem = fabs (fmod (val, 2.0));
 	  if (rem == 1.0)
@@ -758,7 +754,7 @@ SCM_PRIMITIVE_GENERIC (scm_even_p, "even?", 1, 0, 0,
   else if (SCM_REALP (n))
     {
       double val = SCM_REAL_VALUE (n);
-      if (DOUBLE_IS_FINITE (val))
+      if (isfinite (val))
 	{
 	  double rem = fabs (fmod (val, 2.0));
 	  if (rem == 1.0)
@@ -778,7 +774,7 @@ SCM_PRIMITIVE_GENERIC (scm_finite_p, "finite?", 1, 0, 0,
 #define FUNC_NAME s_scm_finite_p
 {
   if (SCM_REALP (x))
-    return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x)));
+    return scm_from_bool (isfinite (SCM_REAL_VALUE (x)));
   else if (scm_is_real (x))
     return SCM_BOOL_T;
   else
@@ -876,7 +872,7 @@ SCM_DEFINE (scm_inf, "inf", 0, 0, 0,
       guile_ieee_init ();
       initialized = 1;
     }
-  return scm_from_double (guile_Inf);
+  return scm_i_from_double (guile_Inf);
 }
 #undef FUNC_NAME
 
@@ -891,7 +887,7 @@ SCM_DEFINE (scm_nan, "nan", 0, 0, 0,
       guile_ieee_init ();
       initialized = 1;
     }
-  return scm_from_double (guile_NaN);
+  return scm_i_from_double (guile_NaN);
 }
 #undef FUNC_NAME
 
@@ -916,7 +912,7 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
       double xx = SCM_REAL_VALUE (x);
       /* If x is a NaN then xx<0 is false so we return x unchanged */
       if (xx < 0.0)
-        return scm_from_double (-xx);
+        return scm_i_from_double (-xx);
       /* Handle signed zeroes properly */
       else if (SCM_UNLIKELY (xx == 0.0))
 	return flo0;
@@ -1311,7 +1307,7 @@ scm_i_inexact_floor_quotient (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_floor_quotient);  /* or return a NaN? */
   else
-    return scm_from_double (floor (x / y));
+    return scm_i_from_double (floor (x / y));
 }
 
 static SCM
@@ -1474,7 +1470,7 @@ scm_i_inexact_floor_remainder (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_floor_remainder);  /* or return a NaN? */
   else
-    return scm_from_double (x - y * floor (x / y));
+    return scm_i_from_double (x - y * floor (x / y));
 }
 
 static SCM
@@ -1678,8 +1674,8 @@ scm_i_inexact_floor_divide (double x, double y, SCM *qp, SCM *rp)
     {
       double q = floor (x / y);
       double r = x - q * y;
-      *qp = scm_from_double (q);
-      *rp = scm_from_double (r);
+      *qp = scm_i_from_double (q);
+      *rp = scm_i_from_double (r);
     }
 }
 
@@ -1844,7 +1840,7 @@ scm_i_inexact_ceiling_quotient (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_ceiling_quotient);  /* or return a NaN? */
   else
-    return scm_from_double (ceil (x / y));
+    return scm_i_from_double (ceil (x / y));
 }
 
 static SCM
@@ -2017,7 +2013,7 @@ scm_i_inexact_ceiling_remainder (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_ceiling_remainder);  /* or return a NaN? */
   else
-    return scm_from_double (x - y * ceil (x / y));
+    return scm_i_from_double (x - y * ceil (x / y));
 }
 
 static SCM
@@ -2230,8 +2226,8 @@ scm_i_inexact_ceiling_divide (double x, double y, SCM *qp, SCM *rp)
     {
       double q = ceil (x / y);
       double r = x - q * y;
-      *qp = scm_from_double (q);
-      *rp = scm_from_double (r);
+      *qp = scm_i_from_double (q);
+      *rp = scm_i_from_double (r);
     }
 }
 
@@ -2376,7 +2372,7 @@ scm_i_inexact_truncate_quotient (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_truncate_quotient);  /* or return a NaN? */
   else
-    return scm_from_double (trunc (x / y));
+    return scm_i_from_double (trunc (x / y));
 }
 
 static SCM
@@ -2511,7 +2507,7 @@ scm_i_inexact_truncate_remainder (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_truncate_remainder);  /* or return a NaN? */
   else
-    return scm_from_double (x - y * trunc (x / y));
+    return scm_i_from_double (x - y * trunc (x / y));
 }
 
 static SCM
@@ -2689,8 +2685,8 @@ scm_i_inexact_truncate_divide (double x, double y, SCM *qp, SCM *rp)
     {
       double q = trunc (x / y);
       double r = x - q * y;
-      *qp = scm_from_double (q);
-      *rp = scm_from_double (r);
+      *qp = scm_i_from_double (q);
+      *rp = scm_i_from_double (r);
     }
 }
 
@@ -2864,9 +2860,9 @@ 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));
+    return scm_i_from_double (floor (x/y + 0.5));
   else if (SCM_LIKELY (y < 0))
-    return scm_from_double (ceil (x/y - 0.5));
+    return scm_i_from_double (ceil (x/y - 0.5));
   else if (y == 0)
     scm_num_overflow (s_scm_centered_quotient);  /* or return a NaN? */
   else
@@ -3086,7 +3082,7 @@ scm_i_inexact_centered_remainder (double x, double y)
     scm_num_overflow (s_scm_centered_remainder);  /* or return a NaN? */
   else
     return scm_nan ();
-  return scm_from_double (x - q * y);
+  return scm_i_from_double (x - q * y);
 }
 
 /* Assumes that both x and y are bigints, though
@@ -3335,8 +3331,8 @@ scm_i_inexact_centered_divide (double x, double y, SCM *qp, SCM *rp)
   else
     q = guile_NaN;
   r = x - q * y;
-  *qp = scm_from_double (q);
-  *rp = scm_from_double (r);
+  *qp = scm_i_from_double (q);
+  *rp = scm_i_from_double (r);
 }
 
 /* Assumes that both x and y are bigints, though
@@ -3564,7 +3560,7 @@ scm_i_inexact_round_quotient (double x, double y)
   if (SCM_UNLIKELY (y == 0))
     scm_num_overflow (s_scm_round_quotient);  /* or return a NaN? */
   else
-    return scm_from_double (scm_c_round (x / y));
+    return scm_i_from_double (scm_c_round (x / y));
 }
 
 /* Assumes that both x and y are bigints, though
@@ -3775,7 +3771,7 @@ scm_i_inexact_round_remainder (double x, double y)
   else
     {
       double q = scm_c_round (x / y);
-      return scm_from_double (x - q * y);
+      return scm_i_from_double (x - q * y);
     }
 }
 
@@ -4006,8 +4002,8 @@ scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp)
     {
       double q = scm_c_round (x / y);
       double r = x - q * y;
-      *qp = scm_from_double (q);
-      *rp = scm_from_double (r);
+      *qp = scm_i_from_double (q);
+      *rp = scm_i_from_double (r);
     }
 }
 
@@ -5354,7 +5350,7 @@ idbl2str (double dbl, char *a, int radix)
     }
   else if (dbl == 0.0)
     {
-      if (!double_is_non_negative_zero (dbl))
+      if (copysign (1.0, dbl) < 0.0)
         a[ch++] = '-';
       strcpy (a + ch, "0.0");
       return ch + 3;
@@ -5566,7 +5562,7 @@ icmplx2str (double real, double imag, char *str, int radix)
 #endif
   /* Don't output a '+' for negative numbers or for Inf and
      NaN.  They will provide their own sign. */
-  if (sgn >= 0 && DOUBLE_IS_FINITE (imag))
+  if (sgn >= 0 && isfinite (imag))
     str[i++] = '+';
   i += idbl2str (imag, &str[i], radix);
   str[i++] = 'i';
@@ -6506,7 +6502,7 @@ SCM_DEFINE (scm_rational_p, "rational?", 1, 0, 0,
   else if (SCM_REALP (x))
     /* 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)));
+    return scm_from_bool (isfinite (SCM_REAL_VALUE (x)));
   else
     return SCM_BOOL_F;
 }
@@ -7181,7 +7177,7 @@ scm_max (SCM x, SCM y)
 	  double yyd = SCM_REAL_VALUE (y);
 
 	  if (xxd > yyd)
-	    return scm_from_double (xxd);
+	    return scm_i_from_double (xxd);
 	  /* If y is a NaN, then "==" is false and we return the NaN */
 	  else if (SCM_LIKELY (!(xxd == yyd)))
 	    return y;
@@ -7220,7 +7216,7 @@ scm_max (SCM x, SCM y)
         big_real:
           xx = scm_i_big2dbl (x);
           yy = SCM_REAL_VALUE (y);
-	  return (xx > yy ? scm_from_double (xx) : y);
+	  return (xx > yy ? scm_i_from_double (xx) : y);
 	}
       else if (SCM_FRACTIONP (y))
 	{
@@ -7238,7 +7234,7 @@ scm_max (SCM x, SCM y)
 	  double yyd = yy;
 
 	  if (yyd > xxd)
-	    return scm_from_double (yyd);
+	    return scm_i_from_double (yyd);
 	  /* If x is a NaN, then "==" is false and we return the NaN */
 	  else if (SCM_LIKELY (!(xxd == yyd)))
 	    return x;
@@ -7269,16 +7265,16 @@ scm_max (SCM x, SCM y)
 	  else if (SCM_UNLIKELY (xx != yy))
 	    return (xx != xx) ? x : y;  /* Return the NaN */
 	  /* xx == yy, but handle signed zeroes properly */
-	  else if (double_is_non_negative_zero (yy))
-	    return y;
-	  else
+	  else if (copysign (1.0, yy) < 0.0)
 	    return x;
+	  else
+	    return y;
 	}
       else if (SCM_FRACTIONP (y))
 	{
 	  double yy = scm_i_fraction2double (y);
 	  double xx = SCM_REAL_VALUE (x);
-	  return (xx < yy) ? scm_from_double (yy) : x;
+	  return (xx < yy) ? scm_i_from_double (yy) : x;
 	}
       else
 	return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max);
@@ -7297,7 +7293,7 @@ scm_max (SCM x, SCM y)
 	{
 	  double xx = scm_i_fraction2double (x);
 	  /* if y==NaN then ">" is false, so we return the NaN y */
-	  return (xx > SCM_REAL_VALUE (y)) ? scm_from_double (xx) : y;
+	  return (xx > SCM_REAL_VALUE (y)) ? scm_i_from_double (xx) : y;
 	}
       else if (SCM_FRACTIONP (y))
 	{
@@ -7359,7 +7355,7 @@ scm_min (SCM x, SCM y)
 	{
 	  double z = xx;
 	  /* if y==NaN then "<" is false and we return NaN */
-	  return (z < SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
+	  return (z < SCM_REAL_VALUE (y)) ? scm_i_from_double (z) : y;
 	}
       else if (SCM_FRACTIONP (y))
 	{
@@ -7390,7 +7386,7 @@ scm_min (SCM x, SCM y)
         big_real:
           xx = scm_i_big2dbl (x);
           yy = SCM_REAL_VALUE (y);
-	  return (xx < yy ? scm_from_double (xx) : y);
+	  return (xx < yy ? scm_i_from_double (xx) : y);
 	}
       else if (SCM_FRACTIONP (y))
 	{
@@ -7405,7 +7401,7 @@ scm_min (SCM x, SCM y)
 	{
 	  double z = SCM_I_INUM (y);
 	  /* if x==NaN then "<" is false and we return NaN */
-	  return (z < SCM_REAL_VALUE (x)) ? scm_from_double (z) : x;
+	  return (z < SCM_REAL_VALUE (x)) ? scm_i_from_double (z) : x;
 	}
       else if (SCM_BIGP (y))
 	{
@@ -7428,16 +7424,16 @@ scm_min (SCM x, SCM y)
 	  else if (SCM_UNLIKELY (xx != yy))
 	    return (xx != xx) ? x : y;  /* Return the NaN */
 	  /* xx == yy, but handle signed zeroes properly */
-	  else if (double_is_non_negative_zero (xx))
-	    return y;
-	  else
+	  else if (copysign (1.0, xx) < 0.0)
 	    return x;
+	  else
+	    return y;
 	}
       else if (SCM_FRACTIONP (y))
 	{
 	  double yy = scm_i_fraction2double (y);
 	  double xx = SCM_REAL_VALUE (x);
-	  return (yy < xx) ? scm_from_double (yy) : x;
+	  return (yy < xx) ? scm_i_from_double (yy) : x;
 	}
       else
 	return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min);
@@ -7456,7 +7452,7 @@ scm_min (SCM x, SCM y)
 	{
 	  double xx = scm_i_fraction2double (x);
 	  /* if y==NaN then "<" is false, so we return the NaN y */
-	  return (xx < SCM_REAL_VALUE (y)) ? scm_from_double (xx) : y;
+	  return (xx < SCM_REAL_VALUE (y)) ? scm_i_from_double (xx) : y;
 	}
       else if (SCM_FRACTIONP (y))
 	{
@@ -7515,7 +7511,7 @@ scm_sum (SCM x, SCM y)
       else if (SCM_REALP (y))
         {
           scm_t_inum xx = SCM_I_INUM (x);
-          return scm_from_double (xx + SCM_REAL_VALUE (y));
+          return scm_i_from_double (xx + SCM_REAL_VALUE (y));
         }
       else if (SCM_COMPLEXP (y))
         {
@@ -7579,7 +7575,7 @@ scm_sum (SCM x, SCM y)
 	  {
 	    double result = mpz_get_d (SCM_I_BIG_MPZ (x)) + SCM_REAL_VALUE (y);
 	    scm_remember_upto_here_1 (x);
-	    return scm_from_double (result);
+	    return scm_i_from_double (result);
 	  }
 	else if (SCM_COMPLEXP (y))
 	  {
@@ -7598,20 +7594,20 @@ scm_sum (SCM x, SCM y)
   else if (SCM_REALP (x))
     {
       if (SCM_I_INUMP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
       else if (SCM_BIGP (y))
 	{
 	  double result = mpz_get_d (SCM_I_BIG_MPZ (y)) + SCM_REAL_VALUE (x);
 	  scm_remember_upto_here_1 (y);
-	  return scm_from_double (result);
+	  return scm_i_from_double (result);
 	}
       else if (SCM_REALP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (SCM_REAL_VALUE (x) + SCM_COMPLEX_REAL (y),
 				 SCM_COMPLEX_IMAG (y));
       else if (SCM_FRACTIONP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
       else
 	return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum);
     }
@@ -7650,7 +7646,7 @@ scm_sum (SCM x, SCM y)
 					scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
 			       SCM_FRACTION_DENOMINATOR (x));
       else if (SCM_REALP (y))
-	return scm_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
+	return scm_i_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (SCM_COMPLEX_REAL (y) + scm_i_fraction2double (x),
 				 SCM_COMPLEX_IMAG (y));
@@ -7718,7 +7714,7 @@ scm_difference (SCM x, SCM y)
              bignum, but negating that gives a fixnum.  */
           return scm_i_normbig (scm_i_clonebig (x, 0));
         else if (SCM_REALP (x))
-          return scm_from_double (-SCM_REAL_VALUE (x));
+          return scm_i_from_double (-SCM_REAL_VALUE (x));
         else if (SCM_COMPLEXP (x))
           return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
                                    -SCM_COMPLEX_IMAG (x));
@@ -7791,9 +7787,9 @@ scm_difference (SCM x, SCM y)
 	   * (0.0 - 0.0) ==> 0.0, but (- 0.0) ==> -0.0.
 	   */
 	  if (xx == 0)
-	    return scm_from_double (- SCM_REAL_VALUE (y));
+	    return scm_i_from_double (- SCM_REAL_VALUE (y));
 	  else
-	    return scm_from_double (xx - SCM_REAL_VALUE (y));
+	    return scm_i_from_double (xx - SCM_REAL_VALUE (y));
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -7865,7 +7861,7 @@ scm_difference (SCM x, SCM y)
 	{
 	  double result = mpz_get_d (SCM_I_BIG_MPZ (x)) - SCM_REAL_VALUE (y);
 	  scm_remember_upto_here_1 (x);
-	  return scm_from_double (result);
+	  return scm_i_from_double (result);
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -7884,20 +7880,20 @@ scm_difference (SCM x, SCM y)
   else if (SCM_REALP (x))
     {
       if (SCM_I_INUMP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
       else if (SCM_BIGP (y))
 	{
 	  double result = SCM_REAL_VALUE (x) - mpz_get_d (SCM_I_BIG_MPZ (y));
 	  scm_remember_upto_here_1 (x);
-	  return scm_from_double (result);      
+	  return scm_i_from_double (result);
 	}
       else if (SCM_REALP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (SCM_REAL_VALUE (x) - SCM_COMPLEX_REAL (y),
 				 -SCM_COMPLEX_IMAG (y));
       else if (SCM_FRACTIONP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
       else
 	return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference);
     }
@@ -7937,7 +7933,7 @@ scm_difference (SCM x, SCM y)
 					       scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
 			       SCM_FRACTION_DENOMINATOR (x));
       else if (SCM_REALP (y))
-	return scm_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
+	return scm_i_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (scm_i_fraction2double (x) - SCM_COMPLEX_REAL (y),
 				 -SCM_COMPLEX_IMAG (y));
@@ -8017,7 +8013,7 @@ scm_product (SCM x, SCM y)
 	     and we must do the multiplication in order to handle
 	     infinities and NaNs properly. */
 	  else if (SCM_REALP (y))
-	    return scm_from_double (0.0 * SCM_REAL_VALUE (y));
+	    return scm_i_from_double (0.0 * SCM_REAL_VALUE (y));
 	  else if (SCM_COMPLEXP (y))
 	    return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y),
 					   0.0 * SCM_COMPLEX_IMAG (y));
@@ -8069,7 +8065,7 @@ scm_product (SCM x, SCM y)
 	  return result;
 	}
       else if (SCM_REALP (y))
-	return scm_from_double (xx * SCM_REAL_VALUE (y));
+	return scm_i_from_double (xx * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
 				 xx * SCM_COMPLEX_IMAG (y));
@@ -8099,7 +8095,7 @@ scm_product (SCM x, SCM y)
 	{
 	  double result = mpz_get_d (SCM_I_BIG_MPZ (x)) * SCM_REAL_VALUE (y);
 	  scm_remember_upto_here_1 (x);
-	  return scm_from_double (result);
+	  return scm_i_from_double (result);
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -8125,15 +8121,15 @@ scm_product (SCM x, SCM y)
 	{
 	  double result = mpz_get_d (SCM_I_BIG_MPZ (y)) * SCM_REAL_VALUE (x);
 	  scm_remember_upto_here_1 (y);
-	  return scm_from_double (result);
+	  return scm_i_from_double (result);
 	}
       else if (SCM_REALP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	return scm_c_make_rectangular (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL (y),
 				 SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
       else if (SCM_FRACTIONP (y))
-	return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
+	return scm_i_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
       else
 	return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
     }
@@ -8179,7 +8175,7 @@ scm_product (SCM x, SCM y)
 	return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
 			       SCM_FRACTION_DENOMINATOR (x));
       else if (SCM_REALP (y))
-	return scm_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
+	return scm_i_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
 	{
 	  double xx = scm_i_fraction2double (x);
@@ -8283,7 +8279,7 @@ scm_divide (SCM x, SCM y)
 	    scm_num_overflow (s_divide);
 	  else
 #endif
-	    return scm_from_double (1.0 / xx);
+	    return scm_i_from_double (1.0 / xx);
 	}
       else if (SCM_COMPLEXP (x))
 	{
@@ -8320,7 +8316,7 @@ scm_divide (SCM x, SCM y)
 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
 	      scm_num_overflow (s_divide);
 #else
-	      return scm_from_double ((double) xx / (double) yy);
+	      return scm_i_from_double ((double) xx / (double) yy);
 #endif
 	    }
 	  else if (xx % yy != 0)
@@ -8347,7 +8343,7 @@ scm_divide (SCM x, SCM y)
             /* FIXME: Precision may be lost here due to:
                (1) The cast from 'scm_t_inum' to 'double'
                (2) Double rounding */
-	    return scm_from_double ((double) xx / yy);
+	    return scm_i_from_double ((double) xx / yy);
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -8446,7 +8442,7 @@ scm_divide (SCM x, SCM y)
 #endif
             /* FIXME: Precision may be lost here due to:
                (1) scm_i_big2dbl (2) Double rounding */
-	    return scm_from_double (scm_i_big2dbl (x) / yy);
+	    return scm_i_from_double (scm_i_big2dbl (x) / yy);
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -8473,7 +8469,7 @@ scm_divide (SCM x, SCM y)
             /* FIXME: Precision may be lost here due to:
                (1) The cast from 'scm_t_inum' to 'double'
                (2) Double rounding */
-	    return scm_from_double (rx / (double) yy);
+	    return scm_i_from_double (rx / (double) yy);
 	}
       else if (SCM_BIGP (y))
 	{
@@ -8482,7 +8478,7 @@ scm_divide (SCM x, SCM y)
              (2) Double rounding */
 	  double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
 	  scm_remember_upto_here_1 (y);
-	  return scm_from_double (rx / dby);
+	  return scm_i_from_double (rx / dby);
 	}
       else if (SCM_REALP (y))
 	{
@@ -8492,7 +8488,7 @@ scm_divide (SCM x, SCM y)
 	    scm_num_overflow (s_divide);
 	  else
 #endif
-	    return scm_from_double (rx / yy);
+	    return scm_i_from_double (rx / yy);
 	}
       else if (SCM_COMPLEXP (y))
 	{
@@ -8500,7 +8496,7 @@ scm_divide (SCM x, SCM y)
 	  goto complex_div;
 	}
       else if (SCM_FRACTIONP (y))
-	return scm_from_double (rx / scm_i_fraction2double (y));
+	return scm_i_from_double (rx / scm_i_fraction2double (y));
       else
 	return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
@@ -8600,7 +8596,7 @@ scm_divide (SCM x, SCM y)
             /* FIXME: Precision may be lost here due to:
                (1) The conversion from fraction to double
                (2) Double rounding */
-	    return scm_from_double (scm_i_fraction2double (x) / yy);
+	    return scm_i_from_double (scm_i_fraction2double (x) / yy);
 	}
       else if (SCM_COMPLEXP (y)) 
 	{
@@ -8678,7 +8674,7 @@ SCM_PRIMITIVE_GENERIC (scm_truncate_number, "truncate", 1, 0, 0,
   if (SCM_I_INUMP (x) || SCM_BIGP (x))
     return x;
   else if (SCM_REALP (x))
-    return scm_from_double (trunc (SCM_REAL_VALUE (x)));
+    return scm_i_from_double (trunc (SCM_REAL_VALUE (x)));
   else if (SCM_FRACTIONP (x))
     return scm_truncate_quotient (SCM_FRACTION_NUMERATOR (x),
 				  SCM_FRACTION_DENOMINATOR (x));
@@ -8698,7 +8694,7 @@ SCM_PRIMITIVE_GENERIC (scm_round_number, "round", 1, 0, 0,
   if (SCM_I_INUMP (x) || SCM_BIGP (x))
     return x;
   else if (SCM_REALP (x))
-    return scm_from_double (scm_c_round (SCM_REAL_VALUE (x)));
+    return scm_i_from_double (scm_c_round (SCM_REAL_VALUE (x)));
   else if (SCM_FRACTIONP (x))
     return scm_round_quotient (SCM_FRACTION_NUMERATOR (x),
 			       SCM_FRACTION_DENOMINATOR (x));
@@ -8716,7 +8712,7 @@ SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0,
   if (SCM_I_INUMP (x) || SCM_BIGP (x))
     return x;
   else if (SCM_REALP (x))
-    return scm_from_double (floor (SCM_REAL_VALUE (x)));
+    return scm_i_from_double (floor (SCM_REAL_VALUE (x)));
   else if (SCM_FRACTIONP (x))
     return scm_floor_quotient (SCM_FRACTION_NUMERATOR (x),
 			       SCM_FRACTION_DENOMINATOR (x));
@@ -8733,7 +8729,7 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
   if (SCM_I_INUMP (x) || SCM_BIGP (x))
     return x;
   else if (SCM_REALP (x))
-    return scm_from_double (ceil (SCM_REAL_VALUE (x)));
+    return scm_i_from_double (ceil (SCM_REAL_VALUE (x)));
   else if (SCM_FRACTIONP (x))
     return scm_ceiling_quotient (SCM_FRACTION_NUMERATOR (x),
 				 SCM_FRACTION_DENOMINATOR (x));
@@ -8772,7 +8768,7 @@ SCM_PRIMITIVE_GENERIC (scm_expt, "expt", 2, 0, 0,
     }
   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)));
+      return scm_i_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));
@@ -8797,7 +8793,7 @@ SCM_PRIMITIVE_GENERIC (scm_sin, "sin", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* sin(exact0) = exact0 */
   else if (scm_is_real (z))
-    return scm_from_double (sin (scm_to_double (z)));
+    return scm_i_from_double (sin (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y;
       x = SCM_COMPLEX_REAL (z);
@@ -8818,7 +8814,7 @@ SCM_PRIMITIVE_GENERIC (scm_cos, "cos", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return SCM_INUM1;  /* cos(exact0) = exact1 */
   else if (scm_is_real (z))
-    return scm_from_double (cos (scm_to_double (z)));
+    return scm_i_from_double (cos (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y;
       x = SCM_COMPLEX_REAL (z);
@@ -8839,7 +8835,7 @@ SCM_PRIMITIVE_GENERIC (scm_tan, "tan", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* tan(exact0) = exact0 */
   else if (scm_is_real (z))
-    return scm_from_double (tan (scm_to_double (z)));
+    return scm_i_from_double (tan (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y, w;
       x = 2.0 * SCM_COMPLEX_REAL (z);
@@ -8864,7 +8860,7 @@ SCM_PRIMITIVE_GENERIC (scm_sinh, "sinh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* sinh(exact0) = exact0 */
   else if (scm_is_real (z))
-    return scm_from_double (sinh (scm_to_double (z)));
+    return scm_i_from_double (sinh (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y;
       x = SCM_COMPLEX_REAL (z);
@@ -8885,7 +8881,7 @@ SCM_PRIMITIVE_GENERIC (scm_cosh, "cosh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return SCM_INUM1;  /* cosh(exact0) = exact1 */
   else if (scm_is_real (z))
-    return scm_from_double (cosh (scm_to_double (z)));
+    return scm_i_from_double (cosh (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y;
       x = SCM_COMPLEX_REAL (z);
@@ -8906,7 +8902,7 @@ SCM_PRIMITIVE_GENERIC (scm_tanh, "tanh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* tanh(exact0) = exact0 */
   else if (scm_is_real (z))
-    return scm_from_double (tanh (scm_to_double (z)));
+    return scm_i_from_double (tanh (scm_to_double (z)));
   else if (SCM_COMPLEXP (z))
     { double x, y, w;
       x = 2.0 * SCM_COMPLEX_REAL (z);
@@ -8934,7 +8930,7 @@ SCM_PRIMITIVE_GENERIC (scm_asin, "asin", 1, 0, 0,
     {
       double w = scm_to_double (z);
       if (w >= -1.0 && w <= 1.0)
-        return scm_from_double (asin (w));
+        return scm_i_from_double (asin (w));
       else
         return scm_product (scm_c_make_rectangular (0, -1),
                             scm_sys_asinh (scm_c_make_rectangular (0, w)));
@@ -8962,9 +8958,9 @@ SCM_PRIMITIVE_GENERIC (scm_acos, "acos", 1, 0, 0,
     {
       double w = scm_to_double (z);
       if (w >= -1.0 && w <= 1.0)
-        return scm_from_double (acos (w));
+        return scm_i_from_double (acos (w));
       else
-        return scm_sum (scm_from_double (acos (0.0)),
+        return scm_sum (scm_i_from_double (acos (0.0)),
                         scm_product (scm_c_make_rectangular (0, 1),
                                      scm_sys_asinh (scm_c_make_rectangular (0, w))));
     }
@@ -8972,7 +8968,7 @@ SCM_PRIMITIVE_GENERIC (scm_acos, "acos", 1, 0, 0,
     { double x, y;
       x = SCM_COMPLEX_REAL (z);
       y = SCM_COMPLEX_IMAG (z);
-      return scm_sum (scm_from_double (acos (0.0)),
+      return scm_sum (scm_i_from_double (acos (0.0)),
                       scm_product (scm_c_make_rectangular (0, 1),
                                    scm_sys_asinh (scm_c_make_rectangular (-y, x))));
     }
@@ -8993,7 +8989,7 @@ SCM_PRIMITIVE_GENERIC (scm_atan, "atan", 1, 1, 0,
       if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
 	return z;  /* atan(exact0) = exact0 */
       else if (scm_is_real (z))
-        return scm_from_double (atan (scm_to_double (z)));
+        return scm_i_from_double (atan (scm_to_double (z)));
       else if (SCM_COMPLEXP (z))
         {
           double v, w;
@@ -9009,7 +9005,7 @@ SCM_PRIMITIVE_GENERIC (scm_atan, "atan", 1, 1, 0,
   else if (scm_is_real (z))
     {
       if (scm_is_real (y))
-        return scm_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
+        return scm_i_from_double (atan2 (scm_to_double (z), scm_to_double (y)));
       else
         return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan);
     }
@@ -9026,7 +9022,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_asinh, "asinh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* asinh(exact0) = exact0 */
   else if (scm_is_real (z))
-    return scm_from_double (asinh (scm_to_double (z)));
+    return scm_i_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),
@@ -9044,7 +9040,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_acosh, "acosh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM1)))
     return SCM_INUM0;  /* acosh(exact1) = exact0 */
   else if (scm_is_real (z) && scm_to_double (z) >= 1.0)
-    return scm_from_double (acosh (scm_to_double (z)));
+    return scm_i_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),
@@ -9062,7 +9058,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_atanh, "atanh", 1, 0, 0,
   if (SCM_UNLIKELY (scm_is_eq (z, SCM_INUM0)))
     return z;  /* atanh(exact0) = exact0 */
   else 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)));
+    return scm_i_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))),
@@ -9165,7 +9161,7 @@ SCM_PRIMITIVE_GENERIC (scm_real_part, "real-part", 1, 0, 0,
 #define FUNC_NAME s_scm_real_part
 {
   if (SCM_COMPLEXP (z))
-    return scm_from_double (SCM_COMPLEX_REAL (z));
+    return scm_i_from_double (SCM_COMPLEX_REAL (z));
   else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z))
     return z;
   else
@@ -9180,7 +9176,7 @@ SCM_PRIMITIVE_GENERIC (scm_imag_part, "imag-part", 1, 0, 0,
 #define FUNC_NAME s_scm_imag_part
 {
   if (SCM_COMPLEXP (z))
-    return scm_from_double (SCM_COMPLEX_IMAG (z));
+    return scm_i_from_double (SCM_COMPLEX_IMAG (z));
   else if (SCM_I_INUMP (z) || SCM_REALP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z))
     return SCM_INUM0;
   else
@@ -9249,9 +9245,9 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0,
 	return z;
     }
   else if (SCM_REALP (z))
-    return scm_from_double (fabs (SCM_REAL_VALUE (z)));
+    return scm_i_from_double (fabs (SCM_REAL_VALUE (z)));
   else if (SCM_COMPLEXP (z))
-    return scm_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
+    return scm_i_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
   else if (SCM_FRACTIONP (z))
     {
       if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
@@ -9273,7 +9269,7 @@ SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0,
 #define FUNC_NAME s_scm_angle
 {
   /* atan(0,-1) is pi and it'd be possible to have that as a constant like
-     flo0 to save allocating a new flonum with scm_from_double each time.
+     flo0 to save allocating a new flonum with scm_i_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))
@@ -9281,32 +9277,32 @@ SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0,
       if (SCM_I_INUM (z) >= 0)
         return flo0;
       else
-	return scm_from_double (atan2 (0.0, -1.0));
+	return scm_i_from_double (atan2 (0.0, -1.0));
     }
   else if (SCM_BIGP (z))
     {
       int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
       scm_remember_upto_here_1 (z);
       if (sgn < 0)
-	return scm_from_double (atan2 (0.0, -1.0));
+	return scm_i_from_double (atan2 (0.0, -1.0));
       else
         return flo0;
     }
   else if (SCM_REALP (z))
     {
       double x = SCM_REAL_VALUE (z);
-      if (x > 0.0 || double_is_non_negative_zero (x))
+      if (copysign (1.0, x) > 0.0)
         return flo0;
       else
-        return scm_from_double (atan2 (0.0, -1.0));
+        return scm_i_from_double (atan2 (0.0, -1.0));
     }
   else if (SCM_COMPLEXP (z))
-    return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
+    return scm_i_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
   else if (SCM_FRACTIONP (z))
     {
       if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
 	return flo0;
-      else return scm_from_double (atan2 (0.0, -1.0));
+      else return scm_i_from_double (atan2 (0.0, -1.0));
     }
   else
     return scm_wta_dispatch_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle);
@@ -9320,11 +9316,11 @@ SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact, "exact->inexact", 1, 0, 0,
 #define FUNC_NAME s_scm_exact_to_inexact
 {
   if (SCM_I_INUMP (z))
-    return scm_from_double ((double) SCM_I_INUM (z));
+    return scm_i_from_double ((double) SCM_I_INUM (z));
   else if (SCM_BIGP (z))
-    return scm_from_double (scm_i_big2dbl (z));
+    return scm_i_from_double (scm_i_big2dbl (z));
   else if (SCM_FRACTIONP (z))
-    return scm_from_double (scm_i_fraction2double (z));
+    return scm_i_from_double (scm_i_fraction2double (z));
   else if (SCM_INEXACTP (z))
     return z;
   else
@@ -9353,7 +9349,7 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
 	return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1,
                                    s_scm_inexact_to_exact);
 
-      if (!SCM_LIKELY (DOUBLE_IS_FINITE (val)))
+      if (!SCM_LIKELY (isfinite (val)))
 	SCM_OUT_OF_RANGE (1, z);
       else if (val == 0.0)
         return SCM_INUM0;
@@ -9406,89 +9402,190 @@ SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
 {
   SCM_ASSERT_TYPE (scm_is_real (x), x, SCM_ARG1, FUNC_NAME, "real");
   SCM_ASSERT_TYPE (scm_is_real (eps), eps, SCM_ARG2, FUNC_NAME, "real");
-  eps = scm_abs (eps);
-  if (scm_is_false (scm_positive_p (eps)))
-    {
-      /* eps is either zero or a NaN */
-      if (scm_is_true (scm_nan_p (eps)))
-	return scm_nan ();
-      else if (SCM_INEXACTP (eps))
-	return scm_exact_to_inexact (x);
-      else
-	return x;
-    }
-  else if (scm_is_false (scm_finite_p (eps)))
-    {
-      if (scm_is_true (scm_finite_p (x)))
-	return flo0;
-      else
-	return scm_nan ();
-    }
-  else if (scm_is_false (scm_finite_p (x))) /* checks for both inf and nan */
-    return x;
-  else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x, eps)),
-				     scm_ceiling (scm_difference (x, eps)))))
+
+  if (SCM_UNLIKELY (!scm_is_exact (eps) || !scm_is_exact (x)))
     {
-      /* There's an integer within range; we want the one closest to zero */
-      if (scm_is_false (scm_less_p (eps, scm_abs (x))))
-	{
-	  /* zero is within range */
-	  if (SCM_INEXACTP (x) || SCM_INEXACTP (eps))
-	    return flo0;
-	  else
-	    return SCM_INUM0;
-	}
-      else if (scm_is_true (scm_positive_p (x)))
-	return scm_ceiling (scm_difference (x, eps));
+      if (SCM_UNLIKELY (scm_is_false (scm_finite_p (eps))))
+        {
+          if (scm_is_false (scm_nan_p (eps)) && scm_is_true (scm_finite_p (x)))
+            return flo0;
+          else
+            return scm_nan ();
+        }
+      else if (SCM_UNLIKELY (scm_is_false (scm_finite_p (x))))
+        return x;
       else
-	return scm_floor (scm_sum (x, eps));
-    }
-  else
-    {
-      /* Use continued fractions to find closest ratio.  All
-	 arithmetic is done with exact numbers.
+        return scm_exact_to_inexact
+          (scm_rationalize (scm_inexact_to_exact (x),
+                            scm_inexact_to_exact (eps)));
+    }
+  else
+    {
+      /* X and EPS are exact rationals.
+
+         The code that follows is equivalent to the following Scheme code:
+
+         (define (exact-rationalize x eps)
+           (let ((n1  (if (negative? x) -1 1))
+                 (x   (abs x))
+                 (eps (abs eps)))
+             (let ((lo (- x eps))
+                   (hi (+ x eps)))
+               (if (<= lo 0)
+                   0
+                   (let loop ((nlo (numerator lo)) (dlo (denominator lo))
+                              (nhi (numerator hi)) (dhi (denominator hi))
+                              (n1 n1) (d1 0) (n2 0) (d2 1))
+                     (let-values (((qlo rlo) (floor/ nlo dlo))
+                                  ((qhi rhi) (floor/ nhi dhi)))
+                       (let ((n0 (+ n2 (* n1 qlo)))
+                             (d0 (+ d2 (* d1 qlo))))
+                         (cond ((zero? rlo) (/ n0 d0))
+                               ((< qlo qhi) (/ (+ n0 n1) (+ d0 d1)))
+                               (else (loop dhi rhi dlo rlo n0 d0 n1 d1))))))))))
       */
 
-      SCM ex = scm_inexact_to_exact (x);
-      SCM int_part = scm_floor (ex);
-      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;
+      int n1_init = 1;
+      SCM lo, hi;
+
+      eps = scm_abs (eps);
+      if (scm_is_true (scm_negative_p (x)))
+        {
+          n1_init = -1;
+          x = scm_difference (x, SCM_UNDEFINED);
+        }
 
-      ex = scm_difference (ex, int_part);            /* x = x-int_part */
-      rx = scm_divide (ex, SCM_UNDEFINED); 	       /* rx = 1/x */
+      /* X and EPS are non-negative exact rationals. */
 
-      /* We stop after a million iterations just to be absolutely sure
-	 that we don't go into an infinite loop.  The process normally
-	 converges after less than a dozen iterations.
-      */
+      lo = scm_difference (x, eps);
+      hi = scm_sum (x, eps);
 
-      while (++i < 1000000)
-	{
-	  a = scm_sum (scm_product (a1, tt), a2);    /* a = a1*tt + a2 */
-	  b = scm_sum (scm_product (b1, tt), b2);    /* b = b1*tt + b2 */
-	  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))),
-			 eps)))                      /* abs(x-a/b) <= eps */
-	    {
-	      SCM res = scm_sum (int_part, scm_divide (a, b));
-	      if (SCM_INEXACTP (x) || SCM_INEXACTP (eps))
-		return scm_exact_to_inexact (res);
-	      else
-		return res;
-	    }
-	  rx = scm_divide (scm_difference (rx, tt),  /* rx = 1/(rx - tt) */
-			   SCM_UNDEFINED);
-	  tt = scm_floor (rx);                       /* tt = floor (rx) */
-	  a2 = a1;
-	  b2 = b1;
-	  a1 = a;
-	  b1 = b;
-	}
-      scm_num_overflow (s_scm_rationalize);
+      if (scm_is_false (scm_positive_p (lo)))
+        /* If zero is included in the interval, return it.
+           It is the simplest rational of all. */
+        return SCM_INUM0;
+      else
+        {
+          SCM result;
+          mpz_t n0, d0, n1, d1, n2, d2;
+          mpz_t nlo, dlo, nhi, dhi;
+          mpz_t qlo, rlo, qhi, rhi;
+
+          /* LO and HI are positive exact rationals. */
+
+          /* Our approach here follows the method described by Alan
+             Bawden in a message entitled "(rationalize x y)" on the
+             rrrs-authors mailing list, dated 16 Feb 1988 14:08:28 EST:
+
+             http://groups.csail.mit.edu/mac/ftpdir/scheme-mail/HTML/rrrs-1988/msg00063.html
+
+             In brief, we compute the continued fractions of the two
+             endpoints of the interval (LO and HI).  The continued
+             fraction of the result consists of the common prefix of the
+             continued fractions of LO and HI, plus one final term.  The
+             final term of the result is the smallest integer contained
+             in the interval between the remainders of LO and HI after
+             the common prefix has been removed.
+
+             The following code lazily computes the continued fraction
+             representations of LO and HI, and simultaneously converts
+             the continued fraction of the result into a rational
+             number.  We use MPZ functions directly to avoid type
+             dispatch and GC allocation during the loop. */
+
+          mpz_inits (n0, d0, n1, d1, n2, d2,
+                     nlo, dlo, nhi, dhi,
+                     qlo, rlo, qhi, rhi,
+                     NULL);
+
+          /* The variables N1, D1, N2 and D2 are used to compute the
+             resulting rational from its continued fraction.  At each
+             step, N2/D2 and N1/D1 are the last two convergents.  They
+             are normally initialized to 0/1 and 1/0, respectively.
+             However, if we negated X then we must negate the result as
+             well, and we do that by initializing N1/D1 to -1/0. */
+          mpz_set_si (n1, n1_init);
+          mpz_set_ui (d1, 0);
+          mpz_set_ui (n2, 0);
+          mpz_set_ui (d2, 1);
+
+          /* The variables NLO, DLO, NHI, and DHI are used to lazily
+             compute the continued fraction representations of LO and HI
+             using Euclid's algorithm.  Initially, NLO/DLO == LO and
+             NHI/DHI == HI. */
+          scm_to_mpz (scm_numerator   (lo), nlo);
+          scm_to_mpz (scm_denominator (lo), dlo);
+          scm_to_mpz (scm_numerator   (hi), nhi);
+          scm_to_mpz (scm_denominator (hi), dhi);
+
+          /* As long as we're using exact arithmetic, the following loop
+             is guaranteed to terminate. */
+          for (;;)
+            {
+              /* Compute the next terms (QLO and QHI) of the continued
+                 fractions of LO and HI. */
+              mpz_fdiv_qr (qlo, rlo, nlo, dlo);  /* QLO <-- floor (NLO/DLO), RLO <-- NLO - QLO * DLO */
+              mpz_fdiv_qr (qhi, rhi, nhi, dhi);  /* QHI <-- floor (NHI/DHI), RHI <-- NHI - QHI * DHI */
+
+              /* The next term of the result will be either QLO or
+                 QLO+1.  Here we compute the next convergent of the
+                 result based on the assumption that QLO is the next
+                 term.  If that turns out to be wrong, we'll adjust
+                 these later by adding N1 to N0 and D1 to D0. */
+              mpz_set (n0, n2); mpz_addmul (n0, n1, qlo);  /* N0 <-- N2 + (QLO * N1) */
+              mpz_set (d0, d2); mpz_addmul (d0, d1, qlo);  /* D0 <-- D2 + (QLO * D1) */
+
+              /* We stop iterating when an integer is contained in the
+                 interval between the remainders NLO/DLO and NHI/DHI.
+                 There are two cases to consider: either NLO/DLO == QLO
+                 is an integer (indicated by RLO == 0), or QLO < QHI. */
+              if (mpz_sgn (rlo) == 0 || mpz_cmp (qlo, qhi) != 0)
+                break;
+
+              /* Efficiently shuffle variables around for the next
+                 iteration.  First we shift the recent convergents. */
+              mpz_swap (n2, n1); mpz_swap (n1, n0);      /* N2 <-- N1 <-- N0 */
+              mpz_swap (d2, d1); mpz_swap (d1, d0);      /* D2 <-- D1 <-- D0 */
+
+              /* The following shuffling is a bit confusing, so some
+                 explanation is in order.  Conceptually, we're doing a
+                 couple of things here.  After substracting the floor of
+                 NLO/DLO, the remainder is RLO/DLO.  The rest of the
+                 continued fraction will represent the remainder's
+                 reciprocal DLO/RLO.  Similarly for the HI endpoint.
+                 So in the next iteration, the new endpoints will be
+                 DLO/RLO and DHI/RHI.  However, when we take the
+                 reciprocals of these endpoints, their order is
+                 switched.  So in summary, we want NLO/DLO <-- DHI/RHI
+                 and NHI/DHI <-- DLO/RLO. */
+              mpz_swap (nlo, dhi); mpz_swap (dhi, rlo); /* NLO <-- DHI <-- RLO */
+              mpz_swap (nhi, dlo); mpz_swap (dlo, rhi); /* NHI <-- DLO <-- RHI */
+            }
+
+          /* There is now an integer in the interval [NLO/DLO NHI/DHI].
+             The last term of the result will be the smallest integer in
+             that interval, which is ceiling(NLO/DLO).  We have already
+             computed floor(NLO/DLO) in QLO, so now we adjust QLO to be
+             equal to the ceiling.  */
+          if (mpz_sgn (rlo) != 0)
+            {
+              /* If RLO is non-zero, then NLO/DLO is not an integer and
+                 the next term will be QLO+1.  QLO was used in the
+                 computation of N0 and D0 above.  Here we adjust N0 and
+                 D0 to be based on QLO+1 instead of QLO.  */
+              mpz_add (n0, n0, n1);  /* N0 <-- N0 + N1 */
+              mpz_add (d0, d0, d1);  /* D0 <-- D0 + D1 */
+            }
+
+          /* The simplest rational in the interval is N0/D0 */
+          result = scm_i_make_ratio_already_reduced (scm_from_mpz (n0),
+                                                     scm_from_mpz (d0));
+          mpz_clears (n0, d0, n1, d1, n2, d2,
+                      nlo, dlo, nhi, dhi,
+                      qlo, rlo, qhi, rhi,
+                      NULL);
+          return result;
+        }
     }
 }
 #undef FUNC_NAME
@@ -9743,14 +9840,7 @@ scm_to_double (SCM val)
 SCM
 scm_from_double (double val)
 {
-  SCM z;
-
-  z = SCM_PACK_POINTER (scm_gc_malloc_pointerless (sizeof (scm_t_double), "real"));
-
-  SCM_SET_CELL_TYPE (z, scm_tc16_real);
-  SCM_REAL_VALUE (z) = val;
-
-  return z;
+  return scm_i_from_double (val);
 }
 
 int
@@ -9813,8 +9903,8 @@ log_of_shifted_double (double x, long shift)
 {
   double ans = log (fabs (x)) + shift * M_LN2;
 
-  if (x > 0.0 || double_is_non_negative_zero (x))
-    return scm_from_double (ans);
+  if (copysign (1.0, x) > 0.0)
+    return scm_i_from_double (ans);
   else
     return scm_c_make_rectangular (ans, M_PI);
 }
@@ -9846,7 +9936,7 @@ log_of_fraction (SCM n, SCM d)
     return (scm_difference (log_of_exact_integer (n),
 			    log_of_exact_integer (d)));
   else if (scm_is_false (scm_negative_p (n)))
-    return scm_from_double
+    return scm_i_from_double
       (log1p (scm_i_divide2double (scm_difference (n, d), d)));
   else
     return scm_c_make_rectangular
@@ -9929,8 +10019,8 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0,
       {
 	double re = scm_to_double (z);
 	double l = log10 (fabs (re));
-	if (re > 0.0 || double_is_non_negative_zero (re))
-	  return scm_from_double (l);
+	if (copysign (1.0, re) > 0.0)
+	  return scm_i_from_double (l);
 	else
 	  return scm_c_make_rectangular (l, M_LOG10E * M_PI);
       }
@@ -9967,7 +10057,7 @@ SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0,
     {
       /* 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)));
+      return scm_i_from_double (exp (scm_to_double (z)));
     }
   else
     return scm_wta_dispatch_1 (g_scm_exp, z, 1, s_scm_exp);
@@ -10126,7 +10216,7 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
                   if (root == floor (root))
                     return SCM_I_MAKINUM ((scm_t_inum) root);
                   else
-                    return scm_from_double (root);
+                    return scm_i_from_double (root);
                 }
               else
                 {
@@ -10170,7 +10260,7 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
                 return scm_c_make_rectangular
                   (0.0, ldexp (sqrt (-signif), expon / 2));
               else
-                return scm_from_double (ldexp (sqrt (signif), expon / 2));
+                return scm_i_from_double (ldexp (sqrt (signif), expon / 2));
             }
         }
       else if (SCM_FRACTIONP (z))
@@ -10203,7 +10293,7 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
               if (xx < 0)
                 return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift));
               else
-                return scm_from_double (ldexp (sqrt (xx), shift));
+                return scm_i_from_double (ldexp (sqrt (xx), shift));
             }
         }
 
@@ -10213,7 +10303,7 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
         if (xx < 0)
           return scm_c_make_rectangular (0.0, sqrt (-xx));
         else
-          return scm_from_double (sqrt (xx));
+          return scm_i_from_double (sqrt (xx));
       }
     }
   else
@@ -10244,8 +10334,8 @@ scm_init_numbers ()
 
   scm_add_feature ("complex");
   scm_add_feature ("inexact");
-  flo0 = scm_from_double (0.0);
-  flo_log10e = scm_from_double (M_LOG10E);
+  flo0 = scm_i_from_double (0.0);
+  flo_log10e = scm_i_from_double (M_LOG10E);
 
   exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
 
diff --git a/libguile/vm-i-scheme.c b/libguile/vm-i-scheme.c
index 2f1d5fe02..07bb6442d 100644
--- a/libguile/vm-i-scheme.c
+++ b/libguile/vm-i-scheme.c
@@ -219,8 +219,13 @@ VM_DEFINE_FUNCTION (151, ge, "ge?", 2)
  */
 
 /* The maximum/minimum tagged integers.  */
-#define INUM_MAX (INTPTR_MAX - 1)
-#define INUM_MIN (INTPTR_MIN + scm_tc2_int)
+#define INUM_MAX  \
+  ((scm_t_signed_bits) SCM_UNPACK (SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM)))
+#define INUM_MIN  \
+  ((scm_t_signed_bits) SCM_UNPACK (SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM)))
+#define INUM_STEP                                \
+  ((scm_t_signed_bits) SCM_UNPACK (SCM_INUM1)    \
+   - (scm_t_signed_bits) SCM_UNPACK (SCM_INUM0))
 
 #define FUNC2(CFUNC,SFUNC)				\
 {							\
@@ -238,28 +243,36 @@ VM_DEFINE_FUNCTION (151, ge, "ge?", 2)
 /* Assembly tagged integer arithmetic routines.  This code uses the
    `asm goto' feature introduced in GCC 4.5.  */
 
-#if defined __x86_64__ && SCM_GNUC_PREREQ (4, 5)
+#if SCM_GNUC_PREREQ (4, 5) && (defined __x86_64__ || defined __i386__)
+
+# undef _CX
+# ifdef __x86_64__
+#  define _CX "rcx"
+# else
+#  define _CX "ecx"
+# endif
 
 /* The macros below check the CPU's overflow flag to improve fixnum
-   arithmetic.  The %rcx register is explicitly clobbered because `asm
-   goto' can't have outputs, in which case the `r' constraint could be
-   used to let the register allocator choose a register.
+   arithmetic.  The _CX register (%rcx or %ecx) is explicitly
+   clobbered because `asm goto' can't have outputs, in which case the
+   `r' constraint could be used to let the register allocator choose a
+   register.
 
    TODO: Use `cold' label attribute in GCC 4.6.
    http://gcc.gnu.org/ml/gcc-patches/2010-10/msg01777.html  */
 
 # define ASM_ADD(x, y)							\
     {									\
-      asm volatile goto ("mov %1, %%rcx; "				\
-			 "test %[tag], %%cl; je %l[slow_add]; "		\
-			 "test %[tag], %0;   je %l[slow_add]; "		\
-			 "add %0, %%rcx;     jo %l[slow_add]; "		\
-			 "sub %[tag], %%rcx; "				\
-			 "mov %%rcx, (%[vsp])\n"			\
+      asm volatile goto ("mov %1, %%"_CX"; "				\
+			 "test %[tag], %%cl;   je %l[slow_add]; "	\
+			 "test %[tag], %0;     je %l[slow_add]; "	\
+			 "sub %[tag], %%"_CX"; "			\
+			 "add %0, %%"_CX";     jo %l[slow_add]; "	\
+			 "mov %%"_CX", (%[vsp])\n"			\
 			 : /* no outputs */				\
 			 : "r" (x), "r" (y),				\
 			   [vsp] "r" (sp), [tag] "i" (scm_tc2_int)	\
-			 : "rcx", "memory"				\
+			 : _CX, "memory", "cc"                          \
 			 : slow_add);					\
       NEXT;								\
     }									\
@@ -268,24 +281,106 @@ VM_DEFINE_FUNCTION (151, ge, "ge?", 2)
 
 # define ASM_SUB(x, y)							\
     {									\
-      asm volatile goto ("mov %0, %%rcx; "				\
-			 "test %[tag], %%cl; je %l[slow_sub]; "		\
-			 "test %[tag], %1;   je %l[slow_sub]; "		\
-			 "sub %1, %%rcx;     jo %l[slow_sub]; "		\
-			 "add %[tag], %%rcx; "				\
-			 "mov %%rcx, (%[vsp])\n"			\
+      asm volatile goto ("mov %0, %%"_CX"; "				\
+			 "test %[tag], %%cl;   je %l[slow_sub]; "	\
+			 "test %[tag], %1;     je %l[slow_sub]; "	\
+			 "sub %1, %%"_CX";     jo %l[slow_sub]; "	\
+			 "add %[tag], %%"_CX"; "			\
+			 "mov %%"_CX", (%[vsp])\n"			\
 			 : /* no outputs */				\
 			 : "r" (x), "r" (y),				\
 			   [vsp] "r" (sp), [tag] "i" (scm_tc2_int)	\
-			 : "rcx", "memory"				\
+			 : _CX, "memory", "cc"                          \
 			 : slow_sub);					\
       NEXT;								\
     }									\
   slow_sub:								\
     do { } while (0)
 
+# define ASM_MUL(x, y)							\
+    {									\
+      scm_t_signed_bits xx = SCM_I_INUM (x);				\
+      asm volatile goto ("mov %1, %%"_CX"; "				\
+			 "test %[tag], %%cl;   je %l[slow_mul]; "	\
+			 "sub %[tag], %%"_CX"; "			\
+			 "test %[tag], %0;     je %l[slow_mul]; "	\
+			 "imul %2, %%"_CX";    jo %l[slow_mul]; "	\
+			 "add %[tag], %%"_CX"; "			\
+			 "mov %%"_CX", (%[vsp])\n"			\
+			 : /* no outputs */				\
+			 : "r" (x), "r" (y), "r" (xx),			\
+			   [vsp] "r" (sp), [tag] "i" (scm_tc2_int)	\
+			 : _CX, "memory", "cc"				\
+			 : slow_mul);					\
+      NEXT;								\
+    }									\
+  slow_mul:								\
+    do { } while (0)
+
 #endif
 
+#if SCM_GNUC_PREREQ (4, 5) && defined __arm__
+
+# define ASM_ADD(x, y)							\
+    if (SCM_LIKELY (SCM_I_INUMP (x) && SCM_I_INUMP (y)))		\
+      {									\
+	asm volatile goto ("adds r0, %0, %1; bvs %l[slow_add]; "	\
+			   "str r0, [%[vsp]]\n"				\
+			   : /* no outputs */				\
+			   : "r" (x), "r" (y - scm_tc2_int),		\
+			     [vsp] "r" (sp)				\
+			   : "r0", "memory", "cc"			\
+			   : slow_add);					\
+	NEXT;								\
+      }									\
+  slow_add:								\
+    do { } while (0)
+
+# define ASM_SUB(x, y)							\
+    if (SCM_LIKELY (SCM_I_INUMP (x) && SCM_I_INUMP (y)))		\
+      {									\
+	asm volatile goto ("subs r0, %0, %1; bvs %l[slow_sub]; "	\
+			   "str r0, [%[vsp]]\n"				\
+			   : /* no outputs */				\
+			   : "r" (x), "r" (y - scm_tc2_int),		\
+			     [vsp] "r" (sp)				\
+			   : "r0", "memory", "cc"			\
+			   : slow_sub);					\
+	NEXT;								\
+      }									\
+  slow_sub:								\
+    do { } while (0)
+
+# if defined (__ARM_ARCH_3M__)  || defined (__ARM_ARCH_4__)		\
+  || defined (__ARM_ARCH_4T__)  || defined (__ARM_ARCH_5__)		\
+  || defined (__ARM_ARCH_5T__)  || defined (__ARM_ARCH_5E__)		\
+  || defined (__ARM_ARCH_5TE__) || defined (__ARM_ARCH_5TEJ__)		\
+  || defined (__ARM_ARCH_6__)   || defined (__ARM_ARCH_6J__)		\
+  || defined (__ARM_ARCH_6K__)  || defined (__ARM_ARCH_6Z__)		\
+  || defined (__ARM_ARCH_6ZK__) || defined (__ARM_ARCH_6T2__)		\
+  || defined (__ARM_ARCH_6M__)  || defined (__ARM_ARCH_7__)		\
+  || defined (__ARM_ARCH_7A__)  || defined (__ARM_ARCH_7R__)		\
+  || defined (__ARM_ARCH_7M__)  || defined (__ARM_ARCH_7EM__)		\
+  || defined (__ARM_ARCH_8A__)
+
+/* The ARM architectures listed above support the SMULL instruction */
+
+#  define ASM_MUL(x, y)							\
+    if (SCM_LIKELY (SCM_I_INUMP (x) && SCM_I_INUMP (y)))		\
+      {									\
+	scm_t_signed_bits rlo, rhi;					\
+	asm ("smull %0, %1, %2, %3\n"					\
+	     : "=r" (rlo), "=r" (rhi)					\
+	     : "r" (SCM_UNPACK (x) - scm_tc2_int),			\
+	       "r" (SCM_I_INUM (y)));					\
+	if (SCM_LIKELY (SCM_SRS (rlo, 31) == rhi))			\
+	  RETURN (SCM_PACK (rlo + scm_tc2_int));			\
+      }									\
+    do { } while (0)
+
+# endif
+
+#endif
 
 VM_DEFINE_FUNCTION (152, add, "add", 2)
 {
@@ -303,15 +398,14 @@ VM_DEFINE_FUNCTION (153, add1, "add1", 1)
 {
   ARGS1 (x);
 
-  /* Check for overflow.  */
-  if (SCM_LIKELY ((scm_t_intptr) SCM_UNPACK (x) < INUM_MAX))
+  /* Check for overflow.  We must avoid overflow in the signed
+     addition below, even if X is not an inum.  */
+  if (SCM_LIKELY ((scm_t_signed_bits) SCM_UNPACK (x) <= INUM_MAX - INUM_STEP))
     {
       SCM result;
 
-      /* Add the integers without untagging.  */
-      result = SCM_PACK ((scm_t_intptr) SCM_UNPACK (x)
-			 + (scm_t_intptr) SCM_UNPACK (SCM_I_MAKINUM (1))
-			 - scm_tc2_int);
+      /* Add 1 to the integer without untagging.  */
+      result = SCM_PACK ((scm_t_signed_bits) SCM_UNPACK (x) + INUM_STEP);
 
       if (SCM_LIKELY (SCM_I_INUMP (result)))
 	RETURN (result);
@@ -337,15 +431,14 @@ VM_DEFINE_FUNCTION (155, sub1, "sub1", 1)
 {
   ARGS1 (x);
 
-  /* Check for underflow.  */
-  if (SCM_LIKELY ((scm_t_intptr) SCM_UNPACK (x) > INUM_MIN))
+  /* Check for overflow.  We must avoid overflow in the signed
+     subtraction below, even if X is not an inum.  */
+  if (SCM_LIKELY ((scm_t_signed_bits) SCM_UNPACK (x) >= INUM_MIN + INUM_STEP))
     {
       SCM result;
 
-      /* Substract the integers without untagging.  */
-      result = SCM_PACK ((scm_t_intptr) SCM_UNPACK (x)
-			 - (scm_t_intptr) SCM_UNPACK (SCM_I_MAKINUM (1))
-			 + scm_tc2_int);
+      /* Substract 1 from the integer without untagging.  */
+      result = SCM_PACK ((scm_t_signed_bits) SCM_UNPACK (x) - INUM_STEP);
 
       if (SCM_LIKELY (SCM_I_INUMP (result)))
 	RETURN (result);
@@ -355,19 +448,24 @@ VM_DEFINE_FUNCTION (155, sub1, "sub1", 1)
   RETURN (scm_difference (x, SCM_I_MAKINUM (1)));
 }
 
-#undef ASM_ADD
-#undef ASM_SUB
-#undef FUNC2
-#undef INUM_MAX
-#undef INUM_MIN
-
 VM_DEFINE_FUNCTION (156, mul, "mul", 2)
 {
   ARGS2 (x, y);
+#ifdef ASM_MUL
+  ASM_MUL (x, y);
+#endif
   SYNC_REGISTER ();
   RETURN (scm_product (x, y));
 }
 
+#undef ASM_ADD
+#undef ASM_SUB
+#undef ASM_MUL
+#undef FUNC2
+#undef INUM_MAX
+#undef INUM_MIN
+#undef INUM_STEP
+
 VM_DEFINE_FUNCTION (157, div, "div", 2)
 {
   ARGS2 (x, y);
@@ -402,12 +500,11 @@ VM_DEFINE_FUNCTION (161, ash, "ash", 2)
   if (SCM_I_INUMP (x) && SCM_I_INUMP (y))
     {
       if (SCM_I_INUM (y) < 0)
-        {
-          /* Right shift, will be a fixnum. */
-          if (SCM_I_INUM (y) > -SCM_I_FIXNUM_BIT)
-            RETURN (SCM_I_MAKINUM (SCM_I_INUM (x) >> -SCM_I_INUM (y)));
-          /* fall through */
-        }
+        /* Right shift, will be a fixnum. */
+        RETURN (SCM_I_MAKINUM
+                (SCM_SRS (SCM_I_INUM (x),
+                          (-SCM_I_INUM (y) <= SCM_I_FIXNUM_BIT-1)
+                          ? -SCM_I_INUM (y) : SCM_I_FIXNUM_BIT-1)));
       else
         /* Left shift. See comments in scm_ash. */
         {
@@ -433,7 +530,8 @@ VM_DEFINE_FUNCTION (162, logand, "logand", 2)
 {
   ARGS2 (x, y);
   if (SCM_I_INUMP (x) && SCM_I_INUMP (y))
-    RETURN (SCM_I_MAKINUM (SCM_I_INUM (x) & SCM_I_INUM (y)));
+    /* Compute bitwise AND without untagging */
+    RETURN (SCM_PACK (SCM_UNPACK (x) & SCM_UNPACK (y)));
   SYNC_REGISTER ();
   RETURN (scm_logand (x, y));
 }
@@ -442,7 +540,8 @@ VM_DEFINE_FUNCTION (163, logior, "logior", 2)
 {
   ARGS2 (x, y);
   if (SCM_I_INUMP (x) && SCM_I_INUMP (y))
-    RETURN (SCM_I_MAKINUM (SCM_I_INUM (x) | SCM_I_INUM (y)));
+    /* Compute bitwise OR without untagging */
+    RETURN (SCM_PACK (SCM_UNPACK (x) | SCM_UNPACK (y)));
   SYNC_REGISTER ();
   RETURN (scm_logior (x, y));
 }
diff --git a/m4/copysign.m4 b/m4/copysign.m4
new file mode 100644
index 000000000..382a6c6bb
--- /dev/null
+++ b/m4/copysign.m4
@@ -0,0 +1,19 @@
+# copysign.m4 serial 1
+dnl Copyright (C) 2011-2013 Free Software Foundation, Inc.
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+AC_DEFUN([gl_FUNC_COPYSIGN],
+[
+  AC_REQUIRE([gl_MATH_H_DEFAULTS])
+
+  dnl Determine COPYSIGN_LIBM.
+  gl_MATHFUNC([copysign], [double], [(double, double)])
+  if test $gl_cv_func_copysign_no_libm = no \
+     && test $gl_cv_func_copysign_in_libm = no; then
+    HAVE_COPYSIGN=0
+    COPYSIGN_LIBM=
+  fi
+  AC_SUBST([COPYSIGN_LIBM])
+])
diff --git a/m4/gnulib-cache.m4 b/m4/gnulib-cache.m4
index 81e464632..d41ec7b22 100644
--- a/m4/gnulib-cache.m4
+++ b/m4/gnulib-cache.m4
@@ -27,7 +27,7 @@
 
 
 # Specification in the form of a command-line invocation:
-#   gnulib-tool --import --dir=. --local-dir=gnulib-local --lib=libgnu --source-base=lib --m4-base=m4 --doc-base=doc --tests-base=tests --aux-dir=build-aux --lgpl=3 --no-conditional-dependencies --libtool --macro-prefix=gl --no-vc-files accept alignof alloca-opt announce-gen autobuild bind byteswap c-strcase canonicalize-lgpl ceil clock-time close connect dirfd duplocale environ extensions flock floor fpieee frexp fstat full-read full-write func gendocs getaddrinfo getlogin getpeername getsockname getsockopt git-version-gen gitlog-to-changelog gnu-web-doc-update gnupload havelib iconv_open-utf inet_ntop inet_pton isinf isnan ldexp lib-symbol-versions lib-symbol-visibility libunistring listen localcharset locale log1p maintainer-makefile malloc-gnu malloca nl_langinfo nproc open pipe-posix pipe2 poll putenv recv recvfrom regex rename select send sendto setenv setsockopt shutdown socket stat-time stdlib strftime striconveh string sys_stat time times trunc verify vsnprintf warnings wchar
+#   gnulib-tool --import --dir=. --local-dir=gnulib-local --lib=libgnu --source-base=lib --m4-base=m4 --doc-base=doc --tests-base=tests --aux-dir=build-aux --lgpl=3 --no-conditional-dependencies --libtool --macro-prefix=gl --no-vc-files accept alignof alloca-opt announce-gen autobuild bind byteswap c-strcase canonicalize-lgpl ceil clock-time close connect copysign dirfd duplocale environ extensions flock floor fpieee frexp fstat full-read full-write func gendocs getaddrinfo getlogin getpeername getsockname getsockopt git-version-gen gitlog-to-changelog gnu-web-doc-update gnupload havelib iconv_open-utf inet_ntop inet_pton isfinite isinf isnan ldexp lib-symbol-versions lib-symbol-visibility libunistring listen localcharset locale log1p maintainer-makefile malloc-gnu malloca nl_langinfo nproc open pipe-posix pipe2 poll putenv recv recvfrom regex rename select send sendto setenv setsockopt shutdown socket stat-time stdlib strftime striconveh string sys_stat time times trunc verify vsnprintf warnings wchar
 
 # Specification in the form of a few gnulib-tool.m4 macro invocations:
 gl_LOCAL_DIR([gnulib-local])
@@ -45,6 +45,7 @@ gl_MODULES([
   clock-time
   close
   connect
+  copysign
   dirfd
   duplocale
   environ
@@ -71,6 +72,7 @@ gl_MODULES([
   iconv_open-utf
   inet_ntop
   inet_pton
+  isfinite
   isinf
   isnan
   largefile
diff --git a/m4/gnulib-comp.m4 b/m4/gnulib-comp.m4
index 55c003a58..8a1b80175 100644
--- a/m4/gnulib-comp.m4
+++ b/m4/gnulib-comp.m4
@@ -61,6 +61,7 @@ AC_DEFUN([gl_EARLY],
   # Code from module close:
   # Code from module configmake:
   # Code from module connect:
+  # Code from module copysign:
   # Code from module dirent:
   # Code from module dirfd:
   # Code from module dirname-lgpl:
@@ -108,12 +109,15 @@ AC_DEFUN([gl_EARLY],
   # Code from module inet_ntop:
   # Code from module inet_pton:
   # Code from module inline:
+  # Code from module isfinite:
   # Code from module isinf:
   # Code from module isnan:
   # Code from module isnand:
   # Code from module isnand-nolibm:
   # Code from module isnanf:
+  # Code from module isnanf-nolibm:
   # Code from module isnanl:
+  # Code from module isnanl-nolibm:
   # Code from module langinfo:
   # Code from module largefile:
   AC_REQUIRE([AC_SYS_LARGEFILE])
@@ -172,6 +176,7 @@ AC_DEFUN([gl_EARLY],
   # Code from module setsockopt:
   # Code from module shutdown:
   # Code from module signal-h:
+  # Code from module signbit:
   # Code from module size_max:
   # Code from module snippet/_Noreturn:
   # Code from module snippet/arg-nonnull:
@@ -292,6 +297,11 @@ AC_SUBST([LTALLOCA])
     AC_LIBOBJ([connect])
   fi
   gl_SYS_SOCKET_MODULE_INDICATOR([connect])
+  gl_FUNC_COPYSIGN
+  if test $HAVE_COPYSIGN = 0; then
+    AC_LIBOBJ([copysign])
+  fi
+  gl_MATH_MODULE_INDICATOR([copysign])
   gl_DIRENT_H
   gl_FUNC_DIRFD
   if test $ac_cv_func_dirfd = no && test $gl_cv_func_dirfd_macro = no; then
@@ -415,6 +425,11 @@ AC_SUBST([LTALLOCA])
   fi
   gl_ARPA_INET_MODULE_INDICATOR([inet_pton])
   gl_INLINE
+  gl_ISFINITE
+  if test $REPLACE_ISFINITE = 1; then
+    AC_LIBOBJ([isfinite])
+  fi
+  gl_MATH_MODULE_INDICATOR([isfinite])
   gl_ISINF
   if test $REPLACE_ISINF = 1; then
     AC_LIBOBJ([isinf])
@@ -445,6 +460,11 @@ AC_SUBST([LTALLOCA])
     gl_PREREQ_ISNANF
   fi
   gl_MATH_MODULE_INDICATOR([isnanf])
+  gl_FUNC_ISNANF_NO_LIBM
+  if test $gl_func_isnanf_no_libm != yes; then
+    AC_LIBOBJ([isnanf])
+    gl_PREREQ_ISNANF
+  fi
   gl_FUNC_ISNANL
   m4_ifdef([gl_ISNAN], [
     AC_REQUIRE([gl_ISNAN])
@@ -454,6 +474,11 @@ AC_SUBST([LTALLOCA])
     gl_PREREQ_ISNANL
   fi
   gl_MATH_MODULE_INDICATOR([isnanl])
+  gl_FUNC_ISNANL_NO_LIBM
+  if test $gl_func_isnanl_no_libm != yes; then
+    AC_LIBOBJ([isnanl])
+    gl_PREREQ_ISNANL
+  fi
   gl_LANGINFO_H
   AC_REQUIRE([gl_LARGEFILE])
   gl_FUNC_LDEXP
@@ -656,6 +681,13 @@ AC_SUBST([LTALLOCA])
   fi
   gl_SYS_SOCKET_MODULE_INDICATOR([shutdown])
   gl_SIGNAL_H
+  gl_SIGNBIT
+  if test $REPLACE_SIGNBIT = 1; then
+    AC_LIBOBJ([signbitf])
+    AC_LIBOBJ([signbitd])
+    AC_LIBOBJ([signbitl])
+  fi
+  gl_MATH_MODULE_INDICATOR([signbit])
   gl_SIZE_MAX
   gl_FUNC_SNPRINTF
   gl_STDIO_MODULE_INDICATOR([snprintf])
@@ -935,6 +967,7 @@ AC_DEFUN([gl_FILE_LIST], [
   lib/close.c
   lib/config.charset
   lib/connect.c
+  lib/copysign.c
   lib/dirent.in.h
   lib/dirfd.c
   lib/dirname-lgpl.c
@@ -976,11 +1009,14 @@ AC_DEFUN([gl_FILE_LIST], [
   lib/iconveh.h
   lib/inet_ntop.c
   lib/inet_pton.c
+  lib/isfinite.c
   lib/isinf.c
   lib/isnan.c
   lib/isnand-nolibm.h
   lib/isnand.c
+  lib/isnanf-nolibm.h
   lib/isnanf.c
+  lib/isnanl-nolibm.h
   lib/isnanl.c
   lib/itold.c
   lib/langinfo.in.h
@@ -1053,6 +1089,9 @@ AC_DEFUN([gl_FILE_LIST], [
   lib/setsockopt.c
   lib/shutdown.c
   lib/signal.in.h
+  lib/signbitd.c
+  lib/signbitf.c
+  lib/signbitl.c
   lib/size_max.h
   lib/snprintf.c
   lib/socket.c
@@ -1125,6 +1164,7 @@ AC_DEFUN([gl_FILE_LIST], [
   m4/close.m4
   m4/codeset.m4
   m4/configmake.m4
+  m4/copysign.m4
   m4/dirent_h.m4
   m4/dirfd.m4
   m4/dirname.m4
@@ -1163,6 +1203,7 @@ AC_DEFUN([gl_FILE_LIST], [
   m4/inline.m4
   m4/intmax_t.m4
   m4/inttypes_h.m4
+  m4/isfinite.m4
   m4/isinf.m4
   m4/isnan.m4
   m4/isnand.m4
@@ -1228,6 +1269,7 @@ AC_DEFUN([gl_FILE_LIST], [
   m4/servent.m4
   m4/setenv.m4
   m4/signal_h.m4
+  m4/signbit.m4
   m4/size_max.m4
   m4/snprintf.m4
   m4/socketlib.m4
diff --git a/m4/isfinite.m4 b/m4/isfinite.m4
new file mode 100644
index 000000000..b54b403b8
--- /dev/null
+++ b/m4/isfinite.m4
@@ -0,0 +1,165 @@
+# isfinite.m4 serial 13
+dnl Copyright (C) 2007-2013 Free Software Foundation, Inc.
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+AC_DEFUN([gl_ISFINITE],
+[
+  AC_REQUIRE([gl_MATH_H_DEFAULTS])
+  dnl Persuade glibc <math.h> to declare isfinite.
+  AC_REQUIRE([gl_USE_SYSTEM_EXTENSIONS])
+  AC_CHECK_DECLS([isfinite], , , [[#include <math.h>]])
+  if test "$ac_cv_have_decl_isfinite" = yes; then
+    gl_CHECK_MATH_LIB([ISFINITE_LIBM],
+     [x = isfinite (x) + isfinite ((float) x);])
+    if test "$ISFINITE_LIBM" != missing; then
+      dnl Test whether isfinite() on 'long double' works.
+      gl_ISFINITEL_WORKS
+      case "$gl_cv_func_isfinitel_works" in
+        *yes) ;;
+        *)    ISFINITE_LIBM=missing;;
+      esac
+      dnl Also, isfinite() on 'double' does not work on Linux/ia64 (because of
+      dnl signalling NaNs). But this does not have to be tested, since
+      dnl isfinite(long double) also does not work in this situation.
+    fi
+  fi
+  if test "$ac_cv_have_decl_isfinite" != yes ||
+     test "$ISFINITE_LIBM" = missing; then
+    REPLACE_ISFINITE=1
+    dnl No libraries are needed to link lib/isfinite.c.
+    ISFINITE_LIBM=
+  fi
+  AC_SUBST([ISFINITE_LIBM])
+])
+
+dnl Test whether isfinite() on 'long double' recognizes all numbers which are
+dnl neither finite nor infinite. This test fails e.g. on i686, x86_64, ia64,
+dnl because of
+dnl - pseudo-denormals on x86_64,
+dnl - pseudo-zeroes, unnormalized numbers, and pseudo-denormals on i686,
+dnl - pseudo-NaN, pseudo-Infinity, pseudo-zeroes, unnormalized numbers, and
+dnl   pseudo-denormals on ia64.
+AC_DEFUN([gl_ISFINITEL_WORKS],
+[
+  AC_REQUIRE([AC_PROG_CC])
+  AC_REQUIRE([gl_BIGENDIAN])
+  AC_REQUIRE([gl_LONG_DOUBLE_VS_DOUBLE])
+  AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
+  AC_CACHE_CHECK([whether isfinite(long double) works], [gl_cv_func_isfinitel_works],
+    [
+      AC_RUN_IFELSE([AC_LANG_SOURCE([[
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#define NWORDS \
+  ((sizeof (long double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
+typedef union { unsigned int word[NWORDS]; long double value; }
+        memory_long_double;
+/* On Irix 6.5, gcc 3.4.3 can't compute compile-time NaN, and needs the
+   runtime type conversion.  */
+#ifdef __sgi
+static long double NaNl ()
+{
+  double zero = 0.0;
+  return zero / zero;
+}
+#else
+# define NaNl() (0.0L / 0.0L)
+#endif
+int main ()
+{
+  int result = 0;
+
+  {
+    memory_long_double m;
+    unsigned int i;
+
+    /* The isfinite macro should be immune against changes in the sign bit and
+       in the mantissa bits.  The xor operation twiddles a bit that can only be
+       a sign bit or a mantissa bit (since the exponent never extends to
+       bit 31).  */
+    m.value = NaNl ();
+    m.word[NWORDS / 2] ^= (unsigned int) 1 << (sizeof (unsigned int) * CHAR_BIT - 1);
+    for (i = 0; i < NWORDS; i++)
+      m.word[i] |= 1;
+    if (isfinite (m.value))
+      result |= 1;
+  }
+
+#if ((defined __ia64 && LDBL_MANT_DIG == 64) || (defined __x86_64__ || defined __amd64__) || (defined __i386 || defined __i386__ || defined _I386 || defined _M_IX86 || defined _X86_)) && !HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
+/* Representation of an 80-bit 'long double' as an initializer for a sequence
+   of 'unsigned int' words.  */
+# ifdef WORDS_BIGENDIAN
+#  define LDBL80_WORDS(exponent,manthi,mantlo) \
+     { ((unsigned int) (exponent) << 16) | ((unsigned int) (manthi) >> 16), \
+       ((unsigned int) (manthi) << 16) | (unsigned int) (mantlo) >> 16),    \
+       (unsigned int) (mantlo) << 16                                        \
+     }
+# else
+#  define LDBL80_WORDS(exponent,manthi,mantlo) \
+     { mantlo, manthi, exponent }
+# endif
+  { /* Quiet NaN.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0xFFFF, 0xC3333333, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 2;
+  }
+  {
+    /* Signalling NaN.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0xFFFF, 0x83333333, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 2;
+  }
+  /* The isfinite macro should recognize Pseudo-NaNs, Pseudo-Infinities,
+     Pseudo-Zeroes, Unnormalized Numbers, and Pseudo-Denormals, as defined in
+       Intel IA-64 Architecture Software Developer's Manual, Volume 1:
+       Application Architecture.
+       Table 5-2 "Floating-Point Register Encodings"
+       Figure 5-6 "Memory to Floating-Point Register Data Translation"
+   */
+  { /* Pseudo-NaN.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0xFFFF, 0x40000001, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 4;
+  }
+  { /* Pseudo-Infinity.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0xFFFF, 0x00000000, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 8;
+  }
+  { /* Pseudo-Zero.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0x4004, 0x00000000, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 16;
+  }
+  { /* Unnormalized number.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0x4000, 0x63333333, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 32;
+  }
+  { /* Pseudo-Denormal.  */
+    static memory_long_double x =
+      { LDBL80_WORDS (0x0000, 0x83333333, 0x00000000) };
+    if (isfinite (x.value))
+      result |= 64;
+  }
+#endif
+
+  return result;
+}]])], [gl_cv_func_isfinitel_works=yes], [gl_cv_func_isfinitel_works=no],
+      [case "$host_cpu" in
+                               # Guess no on ia64, x86_64, i386.
+         ia64 | x86_64 | i*86) gl_cv_func_isfinitel_works="guessing no";;
+         *)                    gl_cv_func_isfinitel_works="guessing yes";;
+       esac
+      ])
+    ])
+])
diff --git a/m4/signbit.m4 b/m4/signbit.m4
new file mode 100644
index 000000000..d58caaf05
--- /dev/null
+++ b/m4/signbit.m4
@@ -0,0 +1,365 @@
+# signbit.m4 serial 13
+dnl Copyright (C) 2007-2013 Free Software Foundation, Inc.
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+AC_DEFUN([gl_SIGNBIT],
+[
+  AC_REQUIRE([gl_MATH_H_DEFAULTS])
+  AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
+  AC_CACHE_CHECK([for signbit macro], [gl_cv_func_signbit],
+    [
+      AC_RUN_IFELSE(
+        [AC_LANG_SOURCE([[
+#include <math.h>
+/* If signbit is defined as a function, don't use it, since calling it for
+   'float' or 'long double' arguments would involve conversions.
+   If signbit is not declared at all but exists as a library function, don't
+   use it, since the prototype may not match.
+   If signbit is not declared at all but exists as a compiler built-in, don't
+   use it, since it's preferable to use __builtin_signbit* (no warnings,
+   no conversions).  */
+#ifndef signbit
+# error "signbit should be a macro"
+#endif
+#include <string.h>
+]gl_SIGNBIT_TEST_PROGRAM
+])],
+        [gl_cv_func_signbit=yes],
+        [gl_cv_func_signbit=no],
+        [case "$host_os" in
+                   # Guess yes on glibc systems.
+           *-gnu*) gl_cv_func_signbit="guessing yes" ;;
+                   # If we don't know, assume the worst.
+           *)      gl_cv_func_signbit="guessing no" ;;
+         esac
+        ])
+    ])
+  dnl GCC 4.0 and newer provides three built-ins for signbit.
+  dnl They can be used without warnings, also in C++, regardless of <math.h>.
+  dnl But they may expand to calls to functions, which may or may not be in
+  dnl libc.
+  AC_CACHE_CHECK([for signbit compiler built-ins], [gl_cv_func_signbit_gcc],
+    [
+      AC_RUN_IFELSE(
+        [AC_LANG_SOURCE([[
+#if __GNUC__ >= 4
+# define signbit(x) \
+   (sizeof (x) == sizeof (long double) ? __builtin_signbitl (x) : \
+    sizeof (x) == sizeof (double) ? __builtin_signbit (x) : \
+    __builtin_signbitf (x))
+#else
+# error "signbit should be three compiler built-ins"
+#endif
+#include <string.h>
+]gl_SIGNBIT_TEST_PROGRAM
+])],
+        [gl_cv_func_signbit_gcc=yes],
+        [gl_cv_func_signbit_gcc=no],
+        [case "$host_os" in
+                   # Guess yes on glibc systems.
+           *-gnu*) gl_cv_func_signbit_gcc="guessing yes" ;;
+                   # If we don't know, assume the worst.
+           *)      gl_cv_func_signbit_gcc="guessing no" ;;
+         esac
+        ])
+    ])
+  dnl Use the compiler built-ins whenever possible, because they are more
+  dnl efficient than the system library functions (if they exist).
+  case "$gl_cv_func_signbit_gcc" in
+    *yes)
+      REPLACE_SIGNBIT_USING_GCC=1
+      ;;
+    *)
+      case "$gl_cv_func_signbit" in
+        *yes) ;;
+        *)
+          dnl REPLACE_SIGNBIT=1 makes sure the signbit[fdl] functions get built.
+          REPLACE_SIGNBIT=1
+          gl_FLOAT_SIGN_LOCATION
+          gl_DOUBLE_SIGN_LOCATION
+          gl_LONG_DOUBLE_SIGN_LOCATION
+          if test "$gl_cv_cc_float_signbit" = unknown; then
+            dnl Test whether copysignf() is declared.
+            AC_CHECK_DECLS([copysignf], , , [[#include <math.h>]])
+            if test "$ac_cv_have_decl_copysignf" = yes; then
+              dnl Test whether copysignf() can be used without libm.
+              AC_CACHE_CHECK([whether copysignf can be used without linking with libm],
+                [gl_cv_func_copysignf_no_libm],
+                [
+                  AC_LINK_IFELSE(
+                    [AC_LANG_PROGRAM(
+                       [[#include <math.h>
+                         float x, y;]],
+                       [[return copysignf (x, y) < 0;]])],
+                    [gl_cv_func_copysignf_no_libm=yes],
+                    [gl_cv_func_copysignf_no_libm=no])
+                ])
+              if test $gl_cv_func_copysignf_no_libm = yes; then
+                AC_DEFINE([HAVE_COPYSIGNF_IN_LIBC], [1],
+                  [Define if the copysignf function is declared in <math.h> and available in libc.])
+              fi
+            fi
+          fi
+          if test "$gl_cv_cc_double_signbit" = unknown; then
+            dnl Test whether copysign() is declared.
+            AC_CHECK_DECLS([copysign], , , [[#include <math.h>]])
+            if test "$ac_cv_have_decl_copysign" = yes; then
+              dnl Test whether copysign() can be used without libm.
+              AC_CACHE_CHECK([whether copysign can be used without linking with libm],
+                [gl_cv_func_copysign_no_libm],
+                [
+                  AC_LINK_IFELSE(
+                    [AC_LANG_PROGRAM(
+                       [[#include <math.h>
+                         double x, y;]],
+                       [[return copysign (x, y) < 0;]])],
+                    [gl_cv_func_copysign_no_libm=yes],
+                    [gl_cv_func_copysign_no_libm=no])
+                ])
+              if test $gl_cv_func_copysign_no_libm = yes; then
+                AC_DEFINE([HAVE_COPYSIGN_IN_LIBC], [1],
+                  [Define if the copysign function is declared in <math.h> and available in libc.])
+              fi
+            fi
+          fi
+          if test "$gl_cv_cc_long_double_signbit" = unknown; then
+            dnl Test whether copysignl() is declared.
+            AC_CHECK_DECLS([copysignl], , , [[#include <math.h>]])
+            if test "$ac_cv_have_decl_copysignl" = yes; then
+              dnl Test whether copysignl() can be used without libm.
+              AC_CACHE_CHECK([whether copysignl can be used without linking with libm],
+                [gl_cv_func_copysignl_no_libm],
+                [
+                  AC_LINK_IFELSE(
+                    [AC_LANG_PROGRAM(
+                       [[#include <math.h>
+                         long double x, y;]],
+                       [[return copysignl (x, y) < 0;]])],
+                    [gl_cv_func_copysignl_no_libm=yes],
+                    [gl_cv_func_copysignl_no_libm=no])
+                ])
+              if test $gl_cv_func_copysignl_no_libm = yes; then
+                AC_DEFINE([HAVE_COPYSIGNL_IN_LIBC], [1],
+                  [Define if the copysignl function is declared in <math.h> and available in libc.])
+              fi
+            fi
+          fi
+          ;;
+      esac
+      ;;
+  esac
+])
+
+AC_DEFUN([gl_SIGNBIT_TEST_PROGRAM], [[
+/* Global variables.
+   Needed because GCC 4 constant-folds __builtin_signbitl (literal)
+   but cannot constant-fold            __builtin_signbitl (variable).  */
+float vf;
+double vd;
+long double vl;
+int main ()
+{
+/* HP cc on HP-UX 10.20 has a bug with the constant expression -0.0.
+   So we use -p0f and -p0d instead.  */
+float p0f = 0.0f;
+float m0f = -p0f;
+double p0d = 0.0;
+double m0d = -p0d;
+/* On HP-UX 10.20, negating 0.0L does not yield -0.0L.
+   So we use another constant expression instead.
+   But that expression does not work on other platforms, such as when
+   cross-compiling to PowerPC on Mac OS X 10.5.  */
+long double p0l = 0.0L;
+#if defined __hpux || defined __sgi
+long double m0l = -LDBL_MIN * LDBL_MIN;
+#else
+long double m0l = -p0l;
+#endif
+  int result = 0;
+  if (signbit (vf)) /* link check */
+    vf++;
+  {
+    float plus_inf = 1.0f / p0f;
+    float minus_inf = -1.0f / p0f;
+    if (!(!signbit (255.0f)
+          && signbit (-255.0f)
+          && !signbit (p0f)
+          && (memcmp (&m0f, &p0f, sizeof (float)) == 0 || signbit (m0f))
+          && !signbit (plus_inf)
+          && signbit (minus_inf)))
+      result |= 1;
+  }
+  if (signbit (vd)) /* link check */
+    vd++;
+  {
+    double plus_inf = 1.0 / p0d;
+    double minus_inf = -1.0 / p0d;
+    if (!(!signbit (255.0)
+          && signbit (-255.0)
+          && !signbit (p0d)
+          && (memcmp (&m0d, &p0d, sizeof (double)) == 0 || signbit (m0d))
+          && !signbit (plus_inf)
+          && signbit (minus_inf)))
+      result |= 2;
+  }
+  if (signbit (vl)) /* link check */
+    vl++;
+  {
+    long double plus_inf = 1.0L / p0l;
+    long double minus_inf = -1.0L / p0l;
+    if (signbit (255.0L))
+      result |= 4;
+    if (!signbit (-255.0L))
+      result |= 4;
+    if (signbit (p0l))
+      result |= 8;
+    if (!(memcmp (&m0l, &p0l, sizeof (long double)) == 0 || signbit (m0l)))
+      result |= 16;
+    if (signbit (plus_inf))
+      result |= 32;
+    if (!signbit (minus_inf))
+      result |= 64;
+  }
+  return result;
+}
+]])
+
+AC_DEFUN([gl_FLOAT_SIGN_LOCATION],
+[
+  gl_FLOATTYPE_SIGN_LOCATION([float], [gl_cv_cc_float_signbit], [f], [FLT])
+])
+
+AC_DEFUN([gl_DOUBLE_SIGN_LOCATION],
+[
+  gl_FLOATTYPE_SIGN_LOCATION([double], [gl_cv_cc_double_signbit], [], [DBL])
+])
+
+AC_DEFUN([gl_LONG_DOUBLE_SIGN_LOCATION],
+[
+  gl_FLOATTYPE_SIGN_LOCATION([long double], [gl_cv_cc_long_double_signbit], [L], [LDBL])
+])
+
+AC_DEFUN([gl_FLOATTYPE_SIGN_LOCATION],
+[
+  AC_CACHE_CHECK([where to find the sign bit in a '$1'],
+    [$2],
+    [
+      AC_RUN_IFELSE(
+        [AC_LANG_SOURCE([[
+#include <stddef.h>
+#include <stdio.h>
+#define NWORDS \
+  ((sizeof ($1) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
+typedef union { $1 value; unsigned int word[NWORDS]; }
+        memory_float;
+static memory_float plus = { 1.0$3 };
+static memory_float minus = { -1.0$3 };
+int main ()
+{
+  size_t j, k, i;
+  unsigned int m;
+  FILE *fp = fopen ("conftest.out", "w");
+  if (fp == NULL)
+    return 1;
+  /* Find the different bit.  */
+  k = 0; m = 0;
+  for (j = 0; j < NWORDS; j++)
+    {
+      unsigned int x = plus.word[j] ^ minus.word[j];
+      if ((x & (x - 1)) || (x && m))
+        {
+          /* More than one bit difference.  */
+          fprintf (fp, "unknown");
+          return 2;
+        }
+      if (x)
+        {
+          k = j;
+          m = x;
+        }
+    }
+  if (m == 0)
+    {
+      /* No difference.  */
+      fprintf (fp, "unknown");
+      return 3;
+    }
+  /* Now m = plus.word[k] ^ ~minus.word[k].  */
+  if (plus.word[k] & ~minus.word[k])
+    {
+      /* Oh? The sign bit is set in the positive and cleared in the negative
+         numbers?  */
+      fprintf (fp, "unknown");
+      return 4;
+    }
+  for (i = 0; ; i++)
+    if ((m >> i) & 1)
+      break;
+  fprintf (fp, "word %d bit %d", (int) k, (int) i);
+  if (fclose (fp) != 0)
+    return 5;
+  return 0;
+}
+        ]])],
+        [$2=`cat conftest.out`],
+        [$2="unknown"],
+        [
+          dnl When cross-compiling, we don't know. It depends on the
+          dnl ABI and compiler version. There are too many cases.
+          $2="unknown"
+        ])
+      rm -f conftest.out
+    ])
+  case "$]$2[" in
+    word*bit*)
+      word=`echo "$]$2[" | sed -e 's/word //' -e 's/ bit.*//'`
+      bit=`echo "$]$2[" | sed -e 's/word.*bit //'`
+      AC_DEFINE_UNQUOTED([$4][_SIGNBIT_WORD], [$word],
+        [Define as the word index where to find the sign of '$1'.])
+      AC_DEFINE_UNQUOTED([$4][_SIGNBIT_BIT], [$bit],
+        [Define as the bit index in the word where to find the sign of '$1'.])
+      ;;
+  esac
+])
+
+# Expands to code that defines a function signbitf(float).
+# It extracts the sign bit of a non-NaN value.
+AC_DEFUN([gl_FLOAT_SIGNBIT_CODE],
+[
+  gl_FLOATTYPE_SIGNBIT_CODE([float], [f], [f])
+])
+
+# Expands to code that defines a function signbitd(double).
+# It extracts the sign bit of a non-NaN value.
+AC_DEFUN([gl_DOUBLE_SIGNBIT_CODE],
+[
+  gl_FLOATTYPE_SIGNBIT_CODE([double], [d], [])
+])
+
+# Expands to code that defines a function signbitl(long double).
+# It extracts the sign bit of a non-NaN value.
+AC_DEFUN([gl_LONG_DOUBLE_SIGNBIT_CODE],
+[
+  gl_FLOATTYPE_SIGNBIT_CODE([long double], [l], [L])
+])
+
+AC_DEFUN([gl_FLOATTYPE_SIGNBIT_CODE],
+[[
+static int
+signbit$2 ($1 value)
+{
+  typedef union { $1 f; unsigned char b[sizeof ($1)]; } float_union;
+  static float_union plus_one = { 1.0$3 };   /* unused bits are zero here */
+  static float_union minus_one = { -1.0$3 }; /* unused bits are zero here */
+  /* Compute the sign bit mask as the XOR of plus_one and minus_one.  */
+  float_union u;
+  unsigned int i;
+  u.f = value;
+  for (i = 0; i < sizeof ($1); i++)
+    if (u.b[i] & (plus_one.b[i] ^ minus_one.b[i]))
+      return 1;
+  return 0;
+}
+]])
diff --git a/module/rnrs/arithmetic/bitwise.scm b/module/rnrs/arithmetic/bitwise.scm
index 0acbc8cb6..5f66cf1c1 100644
--- a/module/rnrs/arithmetic/bitwise.scm
+++ b/module/rnrs/arithmetic/bitwise.scm
@@ -41,6 +41,18 @@
 	  bitwise-reverse-bit-field)
   (import (rnrs base (6))
 	  (rnrs control (6))
+          (rename (only (srfi srfi-60) bitwise-if
+                                       integer-length
+                                       first-set-bit
+                                       copy-bit
+                                       bit-field
+                                       copy-bit-field
+                                       rotate-bit-field
+                                       reverse-bit-field)
+                  (integer-length bitwise-length)
+                  (first-set-bit bitwise-first-bit-set)
+                  (bit-field bitwise-bit-field)
+                  (reverse-bit-field bitwise-reverse-bit-field))
 	  (rename (only (guile) lognot 
 			        logand 
 				logior
@@ -60,70 +72,21 @@
         (bitwise-not (logcount ei))
         (logcount ei)))
 
-  (define (bitwise-if ei1 ei2 ei3)
-    (bitwise-ior (bitwise-and ei1 ei2) (bitwise-and (bitwise-not ei1) ei3)))
-  
-  (define (bitwise-length ei)
-    (do ((result 0 (+ result 1))
-	 (bits (if (negative? ei) (bitwise-not ei) ei)
-	       (bitwise-arithmetic-shift bits -1)))
-	((zero? bits)
-	 result)))
-
-  (define (bitwise-first-bit-set ei)
-    (define (bitwise-first-bit-set-inner bits count)
-      (cond ((zero? bits) -1)
-	    ((logbit? 0 bits) count)
-	    (else (bitwise-first-bit-set-inner 
-		   (bitwise-arithmetic-shift bits -1) (+ count 1)))))
-    (bitwise-first-bit-set-inner ei 0))
-
   (define (bitwise-bit-set? ei1 ei2) (logbit? ei2 ei1))
 
   (define (bitwise-copy-bit ei1 ei2 ei3)
-    (bitwise-if (bitwise-arithmetic-shift-left 1 ei2) 
-		(bitwise-arithmetic-shift-left ei3 ei2)
-		ei1))
-
-  (define (bitwise-bit-field ei1 ei2 ei3)
-    (bitwise-arithmetic-shift-right 
-     (bitwise-and ei1 (bitwise-not (bitwise-arithmetic-shift-left -1 ei3)))
-     ei2))
+    ;; The specification states that ei3 should be either 0 or 1.
+    ;; However, other values have been tolerated by both Guile 2.0.x and
+    ;; the sample implementation given the R6RS library document, so for
+    ;; backward compatibility we continue to permit it.
+    (copy-bit ei2 ei1 (logbit? 0 ei3)))
 
   (define (bitwise-copy-bit-field ei1 ei2 ei3 ei4)
-    (bitwise-if (bitwise-and (bitwise-arithmetic-shift-left -1 ei2)
-			     (bitwise-not 
-			      (bitwise-arithmetic-shift-left -1 ei3)))
-		(bitwise-arithmetic-shift-left ei4 ei2)
-		ei1))
+    (copy-bit-field ei1 ei4 ei2 ei3))
 
-  (define bitwise-arithmetic-shift-left bitwise-arithmetic-shift)
-  (define (bitwise-arithmetic-shift-right ei1 ei2)
-    (bitwise-arithmetic-shift ei1 (- ei2)))
-  
   (define (bitwise-rotate-bit-field ei1 ei2 ei3 ei4)
-    (let ((width (- ei3 ei2)))
-      (if (positive? width)
-	  (let ((field (bitwise-bit-field ei1 ei2 ei3))
-		(count (modulo ei4 width)))
-	    (bitwise-copy-bit-field 
-	     ei1 ei2 ei3 
-	     (bitwise-ior (bitwise-arithmetic-shift-left field count)
-			  (bitwise-arithmetic-shift-right 
-			   field (- width count)))))
-	  ei1)))
+    (rotate-bit-field ei1 ei4 ei2 ei3))
 
-  (define (bitwise-reverse-bit-field ei1 ei2 ei3)
-    (define (reverse-bit-field-recursive n1 n2 len)
-      (if (> len 0)
-	  (reverse-bit-field-recursive
-	   (bitwise-arithmetic-shift-right n1 1) 
-	   (bitwise-copy-bit (bitwise-arithmetic-shift-left n2 1) 0 n1)
-	   (- len 1))
-	  n2))
-    (let ((width (- ei3 ei2)))
-      (if (positive? width)
-	  (let ((field (bitwise-bit-field ei1 ei2 ei3)))
-	    (bitwise-copy-bit-field
-	     ei1 ei2 ei3 (reverse-bit-field-recursive field 0 width)))
-	  ei1))))
+  (define bitwise-arithmetic-shift-left bitwise-arithmetic-shift)
+  (define (bitwise-arithmetic-shift-right ei1 ei2)
+    (bitwise-arithmetic-shift ei1 (- ei2))))
diff --git a/test-suite/tests/fractions.test b/test-suite/tests/fractions.test
index 3ee1347d8..90fad570c 100644
--- a/test-suite/tests/fractions.test
+++ b/test-suite/tests/fractions.test
@@ -150,7 +150,7 @@
   (testeqv (inexact->exact (exact->inexact 2135445/16777216)) 2135445/16777216)
   (testeq (< (- (exact->inexact 10197734562406803221/17452826108659293487)
 		.584302765576009) .0000001) #t)
-  (testeqv (rationalize #e0.76 1/10) 3/4)
+  (testeqv (rationalize #e0.76 1/10) 2/3)
   (testeqv (rationalize #e0.723 1/10) 2/3)
   (testeqv (rationalize #e0.723 1/100) 5/7)
   (testeqv (rationalize #e-0.723 1/100) -5/7)
@@ -351,7 +351,7 @@
   (test= (- 0+6i 1/4 0.5 7) -7.75+6.0i)
   (testeqv (rationalize #e2.5 1/1000) 5/2)
   (testeqv (rationalize 7/3 1/1000) 7/3)
-  (testeqv (rationalize #e3.14159265 1/10) 22/7)
+  (testeqv (rationalize #e3.14159265 1/10) 16/5)
   (testeqv (numerator (/ 8 -6)) -4)
   (testeqv (denominator (/ 8 -6)) 3)
   (testeqv (gcd (numerator 7/9) (denominator 7/9)) 1)
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index a36d49394..ffbbea26f 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -1431,6 +1431,35 @@
   (pass-if (eqv?  1/3   (rationalize  3/10 -1/10)))
   (pass-if (eqv? -1/3   (rationalize -3/10 -1/10)))
 
+  ;; Prior to Guile 2.0.10, rationalize used a faulty algorithm that
+  ;; incorrectly returned 2/3 and -2/3 in the following cases.
+  (pass-if (eqv?  1/2   (rationalize #e+0.67 1/4)))
+  (pass-if (eqv? -1/2   (rationalize #e-0.67 1/4)))
+
+  (pass-if (eqv?  1     (rationalize #e+0.67 1/3)))
+  (pass-if (eqv? -1     (rationalize #e-0.67 1/3)))
+
+  (pass-if (eqv?  1/2   (rationalize #e+0.66 1/3)))
+  (pass-if (eqv? -1/2   (rationalize #e-0.66 1/3)))
+
+  (pass-if (eqv?  1     (rationalize #e+0.67 2/3)))
+  (pass-if (eqv? -1     (rationalize #e-0.67 2/3)))
+
+  (pass-if (eqv?  0     (rationalize #e+0.66 2/3)))
+  (pass-if (eqv?  0     (rationalize #e-0.66 2/3)))
+
+  ;; Prior to Guile 2.0.10, rationalize used a faulty algorithm that
+  ;; incorrectly computed the following approximations of PI.
+  (with-test-prefix "pi"
+    (define *pi* #e3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679)
+    (pass-if (eqv? 16/5 (rationalize *pi* 1/10)))
+    (pass-if (eqv? 201/64 (rationalize *pi* 1/1000)))
+    (pass-if (eqv? 75948/24175 (rationalize *pi* (expt 10 -7))))
+    (pass-if (eqv? 100798/32085 (rationalize *pi* (expt 10 -8))))
+    (pass-if (eqv? 58466453/18610450 (rationalize *pi* (expt 10 -14))))
+    (pass-if (eqv? 2307954651196778721982809475299879198775111361078/734644782339796933783743757007944508986600750685
+                   (rationalize *pi* (expt 10 -95)))))
+
   (pass-if (test-eqv? (/  1.0 3) (rationalize  0.3  1/10)))
   (pass-if (test-eqv? (/ -1.0 3) (rationalize -0.3  1/10)))
   (pass-if (test-eqv? (/  1.0 3) (rationalize  0.3 -1/10)))
diff --git a/test-suite/tests/r6rs-arithmetic-bitwise.test b/test-suite/tests/r6rs-arithmetic-bitwise.test
index 3b358461c..3e23d81f0 100644
--- a/test-suite/tests/r6rs-arithmetic-bitwise.test
+++ b/test-suite/tests/r6rs-arithmetic-bitwise.test
@@ -62,7 +62,7 @@
 
 (with-test-prefix "bitwise-copy-bit"
   (pass-if "bitwise-copy-bit simple"
-    (eqv? (bitwise-copy-bit #b010 2 #b111) #b110)))
+    (eqv? (bitwise-copy-bit #b010 2 1) #b110)))
 
 (with-test-prefix "bitwise-bit-field"
   (pass-if "bitwise-bit-field simple"
diff --git a/test-suite/tests/r6rs-arithmetic-fixnums.test b/test-suite/tests/r6rs-arithmetic-fixnums.test
index 60c3b87e9..2d9b177f7 100644
--- a/test-suite/tests/r6rs-arithmetic-fixnums.test
+++ b/test-suite/tests/r6rs-arithmetic-fixnums.test
@@ -184,7 +184,7 @@
 
   (pass-if "fxbit-set? is #f on index of unset bit" (not (fxbit-set? 5 1))))
 
-(with-test-prefix "fxcopy-bit" (pass-if "simple" (fx=? (fxcopy-bit 2 2 7) 6)))
+(with-test-prefix "fxcopy-bit" (pass-if "simple" (fx=? (fxcopy-bit 2 2 1) 6)))
 
 (with-test-prefix "fxbit-field" 
   (pass-if "simple" (fx=? (fxbit-field 50 1 4) 1)))