guile-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Guile-commits] GNU Guile branch, master, updated. v2.1.0-131-gd8d7c7b


From: Mark H Weaver
Subject: [Guile-commits] GNU Guile branch, master, updated. v2.1.0-131-gd8d7c7b
Date: Tue, 06 Aug 2013 21:38:28 +0000

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU Guile".

http://git.savannah.gnu.org/cgit/guile.git/commit/?id=d8d7c7bf5706ce7873257eb88f0a5cc01b541858

The branch, master has been updated
       via  d8d7c7bf5706ce7873257eb88f0a5cc01b541858 (commit)
       via  524140436fc03ee439d5c358c8c7a4c2c559684a (commit)
       via  ca7b6f6869072a79175c2d2db1b63f706bdf9b25 (commit)
       via  19374ad2dec5529cbe4c82faf26247a57d09ce4d (commit)
       via  e1592f8a4068f116b88d210e85fc0676f001c61e (commit)
       via  8ba567480518fc333b662d12f52936866307aa90 (commit)
       via  48b6f151c71cdf22a2eda48429d9bb31143dc98b (commit)
       via  afa3c37ddc3682939a793a4798f3b55fc8d658ec (commit)
       via  f91a1864c365abef807714ed0b664849f099152c (commit)
       via  b2df1682df9edeb84acf7aacbb97d038aea7e501 (commit)
       via  4fa65b903bd0ad5ed62dca92df71325c0a110809 (commit)
       via  9f6211707b186e182aa1debfb52323bfa9fd26de (commit)
       via  ddf4ff2475515692f04f1e89256f1d1d993b5ef9 (commit)
       via  4350c15673a49ca1eacee5670b12d72e3272e3f5 (commit)
       via  cb1482e719a41182e3beec062ff6844c2ee19498 (commit)
       via  478fa0d53026f5420de5a1dab8b4f46e67138deb (commit)
       via  00472a22bbbbbeaf2c0e61520d4f155ace05e41c (commit)
       via  93da406f331a1849f05e63387442b9aaf33f9540 (commit)
       via  eba5c07715f556d3c27b85afb01baa8a189d7849 (commit)
       via  d9e7774fda909306ebdaab175684043946b0931b (commit)
       via  620c13e8fc02060e0af8fa38398ee4de745d41e9 (commit)
      from  e7f64971ed62a6b58f86b5d90a15b24733e36a8d (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit d8d7c7bf5706ce7873257eb88f0a5cc01b541858
Merge: e7f6497 5241404
Author: Mark H Weaver <address@hidden>
Date:   Tue Aug 6 17:37:34 2013 -0400

    Merge remote-tracking branch 'origin/stable-2.0'
    
    Conflicts:
        libguile/numbers.c
        libguile/vm-i-scheme.c

-----------------------------------------------------------------------

Summary of changes:
 lib/Makefile.am                               |   48 +++-
 lib/{mbtowc.c => copysign.c}                  |   16 +-
 lib/isfinite.c                                |   51 +++
 lib/{isnand-nolibm.h => isnanf-nolibm.h}      |   25 +-
 lib/{isnand-nolibm.h => isnanl-nolibm.h}      |   20 +-
 lib/signbitd.c                                |   64 +++
 lib/signbitf.c                                |   64 +++
 lib/signbitl.c                                |   64 +++
 libguile/numbers.c                            |  534 +++++++++++++++----------
 libguile/vm-i-scheme.c                        |  191 +++++++---
 m4/copysign.m4                                |   19 +
 m4/gnulib-cache.m4                            |    4 +-
 m4/gnulib-comp.m4                             |   42 ++
 m4/{isinf.m4 => isfinite.m4}                  |  104 +++---
 m4/signbit.m4                                 |  365 +++++++++++++++++
 module/rnrs/arithmetic/bitwise.scm            |   81 +---
 test-suite/tests/fractions.test               |    4 +-
 test-suite/tests/numbers.test                 |   29 ++
 test-suite/tests/r6rs-arithmetic-bitwise.test |    2 +-
 test-suite/tests/r6rs-arithmetic-fixnums.test |    2 +-
 20 files changed, 1312 insertions(+), 417 deletions(-)
 copy lib/{mbtowc.c => copysign.c} (78%)
 create mode 100644 lib/isfinite.c
 copy lib/{isnand-nolibm.h => isnanf-nolibm.h} (62%)
 copy lib/{isnand-nolibm.h => isnanl-nolibm.h} (74%)
 create mode 100644 lib/signbitd.c
 create mode 100644 lib/signbitf.c
 create mode 100644 lib/signbitl.c
 create mode 100644 m4/copysign.m4
 copy m4/{isinf.m4 => isfinite.m4} (61%)
 create mode 100644 m4/signbit.m4

diff --git a/lib/Makefile.am b/lib/Makefile.am
index 8857a90..2ba04b7 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/mbtowc.c b/lib/copysign.c
similarity index 78%
copy from lib/mbtowc.c
copy to lib/copysign.c
index 7777f0a..61efeb8 100644
--- a/lib/mbtowc.c
+++ b/lib/copysign.c
@@ -1,6 +1,5 @@
-/* Convert multibyte character to wide character.
+/* Copy sign into another 'double' number.
    Copyright (C) 2011-2013 Free Software Foundation, Inc.
-   Written by Bruno Haible <address@hidden>, 2011.
 
    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
@@ -17,10 +16,11 @@
 
 #include <config.h>
 
-#include <stdlib.h>
+/* Specification.  */
+#include <math.h>
 
-#include <errno.h>
-#include <string.h>
-#include <wchar.h>
-
-#include "mbtowc-impl.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 0000000..d9eddf5
--- /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 <address@hidden>, 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/isnand-nolibm.h b/lib/isnanf-nolibm.h
similarity index 62%
copy from lib/isnand-nolibm.h
copy to lib/isnanf-nolibm.h
index 3510202..56f4fde 100644
--- a/lib/isnand-nolibm.h
+++ b/lib/isnanf-nolibm.h
@@ -14,20 +14,27 @@
    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_ISNAND_IN_LIBC
-/* Get declaration of isnan macro.  */
+#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 isnand
-#  define isnand(x) __builtin_isnan ((double)(x))
+#  undef isnanf
+#  define isnanf(x) __builtin_isnanf ((float)(x))
+# elif defined isnan
+#  undef isnanf
+#  define isnanf(x) isnan ((float)(x))
 # else
-#  undef isnand
-#  define isnand(x) isnan ((double)(x))
+   /* 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 isnand
-# define isnand rpl_isnand
-extern int isnand (double x);
+# undef isnanf
+# define isnanf rpl_isnanf
+extern int isnanf (float x);
 #endif
diff --git a/lib/isnand-nolibm.h b/lib/isnanl-nolibm.h
similarity index 74%
copy from lib/isnand-nolibm.h
copy to lib/isnanl-nolibm.h
index 3510202..c5d0323 100644
--- a/lib/isnand-nolibm.h
+++ b/lib/isnanl-nolibm.h
@@ -14,20 +14,20 @@
    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_ISNAND_IN_LIBC
-/* Get declaration of isnan macro.  */
+#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 isnand
-#  define isnand(x) __builtin_isnan ((double)(x))
-# else
-#  undef isnand
-#  define isnand(x) isnan ((double)(x))
+#  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 isnand
-# define isnand rpl_isnand
-extern int isnand (double x);
+# 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 0000000..1c813da
--- /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 0000000..817484b
--- /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 0000000..159cfce
--- /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 3c0d765..f549193 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 2f1d5fe..07bb644 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 0000000..382a6c6
--- /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 81e4646..d41ec7b 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 55c003a..8a1b801 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/isinf.m4 b/m4/isfinite.m4
similarity index 61%
copy from m4/isinf.m4
copy to m4/isfinite.m4
index 513a1ba..b54b403 100644
--- a/m4/isinf.m4
+++ b/m4/isfinite.m4
@@ -1,59 +1,55 @@
-# isinf.m4 serial 9
+# 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_ISINF],
+AC_DEFUN([gl_ISFINITE],
 [
   AC_REQUIRE([gl_MATH_H_DEFAULTS])
-  dnl Persuade glibc <math.h> to declare isinf.
+  dnl Persuade glibc <math.h> to declare isfinite.
   AC_REQUIRE([gl_USE_SYSTEM_EXTENSIONS])
-  AC_CHECK_DECLS([isinf], , ,
-    [[#include <math.h>
-      #ifndef isinf
-      #error "isinf must be a macro, not a function"
-      #endif
-    ]])
-  if test "$ac_cv_have_decl_isinf" = yes; then
-    gl_CHECK_MATH_LIB([ISINF_LIBM], [x = isinf (x) + isinf ((float) x);])
-    if test "$ISINF_LIBM" != missing; then
-      dnl Test whether isinf() on 'long double' works.
-      gl_ISINFL_WORKS
-      case "$gl_cv_func_isinfl_works" in
+  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) ;;
-        *)    ISINF_LIBM=missing;;
+        *)    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_isinf" != yes ||
-     test "$ISINF_LIBM" = missing; then
-    REPLACE_ISINF=1
-    dnl No libraries are needed to link lib/isinf.c.
-    ISINF_LIBM=
+  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([ISINF_LIBM])
+  AC_SUBST([ISFINITE_LIBM])
 ])
 
-dnl Test whether isinf() works:
-dnl 1) Whether it correctly returns false for LDBL_MAX.
-dnl 2) Whether on 'long double' recognizes all numbers which are neither
-dnl    finite nor infinite. This test fails on OpenBSD/x86, but could also
-dnl    fail e.g. on i686, x86_64, ia64, 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_ISINFL_WORKS],
+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 isinf(long double) works], [gl_cv_func_isinfl_works],
+  AC_CACHE_CHECK([whether isfinite(long double) works], 
[gl_cv_func_isfinitel_works],
     [
-      AC_RUN_IFELSE(
-        [AC_LANG_SOURCE([[
+      AC_RUN_IFELSE([AC_LANG_SOURCE([[
 #include <float.h>
 #include <limits.h>
 #include <math.h>
@@ -76,14 +72,11 @@ int main ()
 {
   int result = 0;
 
-  if (isinf (LDBL_MAX))
-    result |= 1;
-
   {
     memory_long_double m;
     unsigned int i;
 
-    /* The isinf macro should be immune against changes in the sign bit and
+    /* 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).  */
@@ -91,8 +84,8 @@ int main ()
     m.word[NWORDS / 2] ^= (unsigned int) 1 << (sizeof (unsigned int) * 
CHAR_BIT - 1);
     for (i = 0; i < NWORDS; i++)
       m.word[i] |= 1;
-    if (isinf (m.value))
-      result |= 2;
+    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
@@ -111,17 +104,17 @@ int main ()
   { /* Quiet NaN.  */
     static memory_long_double x =
       { LDBL80_WORDS (0xFFFF, 0xC3333333, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 2;
   }
   {
     /* Signalling NaN.  */
     static memory_long_double x =
       { LDBL80_WORDS (0xFFFF, 0x83333333, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 2;
   }
-  /* The isinf macro should recognize Pseudo-NaNs, Pseudo-Infinities,
+  /* 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.
@@ -131,44 +124,41 @@ int main ()
   { /* Pseudo-NaN.  */
     static memory_long_double x =
       { LDBL80_WORDS (0xFFFF, 0x40000001, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 4;
   }
   { /* Pseudo-Infinity.  */
     static memory_long_double x =
       { LDBL80_WORDS (0xFFFF, 0x00000000, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 8;
   }
   { /* Pseudo-Zero.  */
     static memory_long_double x =
       { LDBL80_WORDS (0x4004, 0x00000000, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 16;
   }
   { /* Unnormalized number.  */
     static memory_long_double x =
       { LDBL80_WORDS (0x4000, 0x63333333, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 32;
   }
   { /* Pseudo-Denormal.  */
     static memory_long_double x =
       { LDBL80_WORDS (0x0000, 0x83333333, 0x00000000) };
-    if (isinf (x.value))
+    if (isfinite (x.value))
       result |= 64;
   }
 #endif
 
   return result;
-}]])], [gl_cv_func_isinfl_works=yes], [gl_cv_func_isinfl_works=no],
-      [
-       case "$host" in
-           # Guess no on OpenBSD ia64, x86_64, i386.
-         ia64-*-openbsd* | x86_64-*-openbsd* | i*86-*-openbsd*)
-            gl_cv_func_isinfl_works="guessing no";;
-         *)
-            gl_cv_func_isinfl_works="guessing yes";;
+}]])], [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 0000000..d58caaf
--- /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 0acbc8c..5f66cf1 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 3ee1347..90fad57 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 a36d493..ffbbea2 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 3b35846..3e23d81 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 60c3b87..2d9b177 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)))


hooks/post-receive
-- 
GNU Guile



reply via email to

[Prev in Thread] Current Thread [Next in Thread]