New module 'fpucw'.

20 files changed, 217 insertions, 1 deletions

diff --git a/ChangeLog b/ChangeLog
index 3409b07c79..861905555b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,35 @@
 2007-03-24  Bruno Haible  <bruno@clisp.org>
 
+	* modules/fpucw: New file.
+	* lib/fpucw.h: New file.
+	* lib/frexp.c: Include fpucw.h.
+	(DECL_ROUNDING, BEGIN_ROUNDING, END_ROUNDING): New macros.
+	(FUNC): Use them.
+	* lib/printf-frexp.c: Include fpucw.h.
+	(DECL_ROUNDING, BEGIN_ROUNDING, END_ROUNDING): New macros.
+	(FUNC): Use them.
+	* lib/vasnprintf.c: Include fpucw.h.
+	(VASNPRINTF): Invoke BEGIN/END_LONG_DOUBLE_ROUNDING around the
+	'long double' calculations.
+	* tests/test-frexpl.c: Include fpucw.h.
+	(main): Invoke BEGIN_LONG_DOUBLE_ROUNDING.
+	* tests/test-printf-frexpl.c: Include fpucw.h.
+	(main): Invoke BEGIN_LONG_DOUBLE_ROUNDING.
+	* modules/frexpl (Depends-on): Add fpucw.
+	* modules/printf-frexpl (Depends-on): Likewise.
+	* modules/fprintf-posix (Depends-on): Likewise.
+	* modules/snprintf-posix (Depends-on): Likewise.
+	* modules/sprintf-posix (Depends-on): Likewise.
+	* modules/vasnprintf-posix (Depends-on): Likewise.
+	* modules/vasprintf-posix (Depends-on): Likewise.
+	* modules/vfprintf-posix (Depends-on): Likewise.
+	* modules/vsnprintf-posix (Depends-on): Likewise.
+	* modules/vsprintf-posix (Depends-on): Likewise.
+	* modules/frexpl-tests (Depends-on): Likewise.
+	* modules/printf-frexpl-tests (Depends-on): Likewise.
+
+2007-03-24  Bruno Haible  <bruno@clisp.org>
+
 	* lib/float+.h: New file.
 	* lib/isnan.c: Include float+.h.
 	(SIZE): New macro.
diff --git a/lib/fpucw.h b/lib/fpucw.h
new file mode 100644
index 0000000000..16aa4bcb34
--- /dev/null
+++ b/lib/fpucw.h
@@ -0,0 +1,108 @@
+/* Manipulating the FPU control word.
+   Copyright (C) 2007 Free Software Foundation, Inc.
+   Written by Bruno Haible <bruno@clisp.org>, 2007.
+
+   This program is free software; you can redistribute it and/or modify
+   it under the terms of the GNU 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 General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program; if not, write to the Free Software Foundation,
+   Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
+
+#ifndef _FPUCW_H
+#define _FPUCW_H
+
+/* The i386 floating point hardware (the 387 compatible FPU, not the modern
+   SSE/SSE2 hardware) has a controllable rounding precision.  It is specified
+   through the 'PC' bits in the FPU control word ('fctrl' register).  (See
+   the GNU libc i386 <fpu_control.h> header for details.)
+
+   On some platforms, such as Linux or Solaris, the default precision setting
+   is set to "extended precision".  This means that 'long double' instructions
+   operate correctly, but 'double' computations often produce slightly
+   different results as on strictly IEEE 754 conforming systems.
+
+   On some platforms, such as NetBSD, the default precision is set to
+   "double precision".  This means that 'long double' instructions will operate
+   only as 'double', i.e. lead wrong results.
+
+   The FPU control word is under control of the application, i.e. it is
+   not required to be set either way by the ABI.  (In fact, the i386 ABI
+   http://refspecs.freestandards.org/elf/abi386-4.pdf page 3-12 = page 38
+   is not clear about it.  But in any case, gcc treats the control word
+   like a "preserved" register: it emits code that assumes that the control
+   word is preserved across calls, and it restores the control word at the
+   end of functions that modify it.)
+
+   See Vincent Lefèvre's page http://www.vinc17.org/research/extended.en.html
+   for a good explanation.
+   See http://www.uwsg.iu.edu/hypermail/linux/kernel/0103.0/0453.html for
+   some argumentation which setting should be the default.  */
+
+/* This header file provides the following facilities:
+     fpucw_t                        integral type holding the value of 'fctrl'
+     FPU_PC_MASK                    bit mask denoting the precision control
+     FPU_PC_DOUBLE                  precision control for 53 bits mantissa
+     FPU_PC_EXTENDED                precision control for 64 bits mantissa
+     GET_FPUCW ()                   yields the current FPU control word
+     SET_FPUCW (word)               sets the FPU control word
+     DECL_LONG_DOUBLE_ROUNDING      variable declaration for
+                                    BEGIN/END_LONG_DOUBLE_ROUNDING
+     BEGIN_LONG_DOUBLE_ROUNDING ()  starts a sequence of instructions with
+                                    'long double' safe operation precision
+     END_LONG_DOUBLE_ROUNDING ()    ends a sequence of instructions with
+                                    'long double' safe operation precision
+ */
+
+/* Inline assembler like this works only with GNU C.  */
+#if defined __i386__ && defined __GNUC__
+
+typedef unsigned short fpucw_t; /* glibc calls this fpu_control_t */
+
+# define FPU_PC_MASK 0x0300
+# define FPU_PC_DOUBLE 0x200    /* glibc calls this _FPU_DOUBLE */
+# define FPU_PC_EXTENDED 0x300  /* glibc calls this _FPU_EXTENDED */
+
+# define GET_FPUCW() \
+  ({ fpucw_t _cw;						\
+     __asm__ __volatile__ ("fnstcw %0" : "=m" (*&_cw));		\
+     _cw;							\
+   })
+# define SET_FPUCW(word) \
+  (void)({ fpucw_t _ncw = (word);				\
+           __asm__ __volatile__ ("fldcw %0" : : "m" (*&_ncw));	\
+         })
+
+# define DECL_LONG_DOUBLE_ROUNDING \
+  fpucw_t oldcw;
+# define BEGIN_LONG_DOUBLE_ROUNDING() \
+  (void)(oldcw = GET_FPUCW (),					\
+         SET_FPUCW ((oldcw & ~FPU_PC_MASK) | FPU_PC_EXTENDED))
+# define END_LONG_DOUBLE_ROUNDING() \
+  SET_FPUCW (oldcw)
+
+#else
+
+typedef unsigned int fpucw_t;
+
+# define FPU_PC_MASK 0
+# define FPU_PC_DOUBLE 0
+# define FPU_PC_EXTENDED 0
+
+# define GET_FPUCW() 0
+# define SET_FPUCW(word) (void)(word)
+
+# define DECL_LONG_DOUBLE_ROUNDING
+# define BEGIN_LONG_DOUBLE_ROUNDING()
+# define END_LONG_DOUBLE_ROUNDING()
+
+#endif
+
+#endif /* _FPUCW_H */
diff --git a/lib/frexp.c b/lib/frexp.c
index 3074f2e151..701f00eefb 100644
--- a/lib/frexp.c
+++ b/lib/frexp.c
@@ -28,6 +28,7 @@
 # include <float.h>
 # ifdef USE_LONG_DOUBLE
 #  include "isnanl-nolibm.h"
+#  include "fpucw.h"
 # else
 #  include "isnan.h"
 # endif
@@ -40,11 +41,17 @@
 #  define FUNC frexpl
 #  define DOUBLE long double
 #  define ISNAN isnanl
+#  define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
+#  define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
+#  define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
 #  define L_(literal) literal##L
 # else
 #  define FUNC frexp
 #  define DOUBLE double
 #  define ISNAN isnan
+#  define DECL_ROUNDING
+#  define BEGIN_ROUNDING()
+#  define END_ROUNDING()
 #  define L_(literal) literal
 # endif
 
@@ -53,6 +60,7 @@ FUNC (DOUBLE x, int *exp)
 {
   int sign;
   int exponent;
+  DECL_ROUNDING
 
   /* Test for NaN, infinity, and zero.  */
   if (ISNAN (x) || x + x == x)
@@ -68,6 +76,8 @@ FUNC (DOUBLE x, int *exp)
       sign = -1;
     }
 
+  BEGIN_ROUNDING ();
+
   {
     /* Since the exponent is an 'int', it fits in 64 bits.  Therefore the
        loops are executed no more than 64 times.  */
@@ -149,8 +159,13 @@ FUNC (DOUBLE x, int *exp)
     /* Here 0.5 <= x < 1.0.  */
   }
 
+  if (sign < 0)
+    x = - x;
+
+  END_ROUNDING ();
+
   *exp = exponent;
-  return (sign < 0 ? - x : x);
+  return x;
 }
 
 #else
diff --git a/lib/printf-frexp.c b/lib/printf-frexp.c
index e6032a4fb0..978653d6c6 100644
--- a/lib/printf-frexp.c
+++ b/lib/printf-frexp.c
@@ -28,6 +28,9 @@
 
 # include <float.h>
 # include <math.h>
+# ifdef USE_LONG_DOUBLE
+#  include "fpucw.h"
+# endif
 
 /* This file assumes FLT_RADIX = 2.  If FLT_RADIX is a power of 2 greater
    than 2, or not even a power of 2, some rounding errors can occur, so that
@@ -42,6 +45,9 @@
 #   define FREXP frexpl
 #   define LDEXP ldexpl
 #  endif
+#  define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
+#  define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
+#  define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
 #  define L_(literal) literal##L
 # else
 #  define FUNC printf_frexp
@@ -52,6 +58,9 @@
 #   define FREXP frexp
 #   define LDEXP ldexp
 #  endif
+#  define DECL_ROUNDING
+#  define BEGIN_ROUNDING()
+#  define END_ROUNDING()
 #  define L_(literal) literal
 # endif
 
@@ -59,6 +68,9 @@ DOUBLE
 FUNC (DOUBLE x, int *exp)
 {
   int exponent;
+  DECL_ROUNDING
+
+  BEGIN_ROUNDING ();
 
 # ifdef USE_FREXP_LDEXP
   /* frexp and ldexp are usually faster than the loop below.  */
@@ -170,6 +182,8 @@ FUNC (DOUBLE x, int *exp)
      or 1.0 <= x < 2.0 and exponent >= MIN_EXP - 1.  */
 # endif
 
+  END_ROUNDING ();
+
   *exp = exponent;
   return x;
 }
diff --git a/lib/vasnprintf.c b/lib/vasnprintf.c
index 00419aaaaf..75d33b347e 100644
--- a/lib/vasnprintf.c
+++ b/lib/vasnprintf.c
@@ -57,6 +57,7 @@
 # if HAVE_LONG_DOUBLE
 #  include "isnanl-nolibm.h"
 #  include "printf-frexpl.h"
+#  include "fpucw.h"
 # endif
 #endif
 
@@ -415,6 +416,9 @@ VASNPRINTF (CHAR_T *resultbuf, size_t *lengthp, const CHAR_T *format, va_list ar
 		    else
 		      {
 			int sign = 0;
+			DECL_LONG_DOUBLE_ROUNDING
+
+			BEGIN_LONG_DOUBLE_ROUNDING ();
 
 			if (arg < 0.0L)
 			  {
@@ -542,6 +546,8 @@ VASNPRINTF (CHAR_T *resultbuf, size_t *lengthp, const CHAR_T *format, va_list ar
 			      while (*p != '\0')
 				p++;
 			  }
+
+			END_LONG_DOUBLE_ROUNDING ();
 		      }
 		  }
 		else
diff --git a/modules/fprintf-posix b/modules/fprintf-posix
index f79a8b428e..2be06d0a48 100644
--- a/modules/fprintf-posix
+++ b/modules/fprintf-posix
@@ -14,6 +14,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_FPRINTF_POSIX
diff --git a/modules/fpucw b/modules/fpucw
new file mode 100644
index 0000000000..d03bad3d4f
--- /dev/null
+++ b/modules/fpucw
@@ -0,0 +1,21 @@
+Description:
+Set the FPU control word, so as to allow correct 'long double' computations.
+
+Files:
+lib/fpucw.h
+
+Depends-on:
+
+configure.ac:
+
+Makefile.am:
+
+Include:
+"fpucw.h"
+
+License:
+LGPL
+
+Maintainer:
+Bruno Haible
+
diff --git a/modules/frexpl b/modules/frexpl
index 0bb8f586d8..465336e34e 100644
--- a/modules/frexpl
+++ b/modules/frexpl
@@ -9,6 +9,7 @@ m4/frexpl.m4
 Depends-on:
 math
 isnanl-nolibm
+fpucw
 
 configure.ac:
 gl_FUNC_FREXPL
diff --git a/modules/frexpl-tests b/modules/frexpl-tests
index 3dcaf11615..461f47cb1a 100644
--- a/modules/frexpl-tests
+++ b/modules/frexpl-tests
@@ -2,6 +2,7 @@ Files:
 tests/test-frexpl.c
 
 Depends-on:
+fpucw
 
 configure.ac:
 
diff --git a/modules/printf-frexpl b/modules/printf-frexpl
index 36d4f346cc..6d137621f1 100644
--- a/modules/printf-frexpl
+++ b/modules/printf-frexpl
@@ -11,6 +11,7 @@ m4/longdouble.m4
 
 Depends-on:
 math
+fpucw
 
 configure.ac:
 gl_FUNC_PRINTF_FREXPL
diff --git a/modules/printf-frexpl-tests b/modules/printf-frexpl-tests
index ab58971c0d..7eb69835a6 100644
--- a/modules/printf-frexpl-tests
+++ b/modules/printf-frexpl-tests
@@ -2,6 +2,7 @@ Files:
 tests/test-printf-frexpl.c
 
 Depends-on:
+fpucw
 
 configure.ac:
 
diff --git a/modules/snprintf-posix b/modules/snprintf-posix
index a5a6ac1fbf..65edc647ad 100644
--- a/modules/snprintf-posix
+++ b/modules/snprintf-posix
@@ -13,6 +13,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_SNPRINTF_POSIX
diff --git a/modules/sprintf-posix b/modules/sprintf-posix
index 498e3d54c1..05d7d4f613 100644
--- a/modules/sprintf-posix
+++ b/modules/sprintf-posix
@@ -13,6 +13,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_SPRINTF_POSIX
diff --git a/modules/vasnprintf-posix b/modules/vasnprintf-posix
index e72337ad37..0f2aeb7277 100644
--- a/modules/vasnprintf-posix
+++ b/modules/vasnprintf-posix
@@ -12,6 +12,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_VASNPRINTF_POSIX
diff --git a/modules/vasprintf-posix b/modules/vasprintf-posix
index e07e3d3f3a..d9bdf6d402 100644
--- a/modules/vasprintf-posix
+++ b/modules/vasprintf-posix
@@ -12,6 +12,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_VASPRINTF_POSIX
diff --git a/modules/vfprintf-posix b/modules/vfprintf-posix
index 1c0adb2062..ce4bde8f23 100644
--- a/modules/vfprintf-posix
+++ b/modules/vfprintf-posix
@@ -14,6 +14,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_VFPRINTF_POSIX
diff --git a/modules/vsnprintf-posix b/modules/vsnprintf-posix
index 4e42310f1a..cfb35a6611 100644
--- a/modules/vsnprintf-posix
+++ b/modules/vsnprintf-posix
@@ -13,6 +13,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_VSNPRINTF_POSIX
diff --git a/modules/vsprintf-posix b/modules/vsprintf-posix
index af073f517c..ef8b22c637 100644
--- a/modules/vsprintf-posix
+++ b/modules/vsprintf-posix
@@ -13,6 +13,7 @@ isnan-nolibm
 isnanl-nolibm
 printf-frexp
 printf-frexpl
+fpucw
 
 configure.ac:
 gl_FUNC_VSPRINTF_POSIX
diff --git a/tests/test-frexpl.c b/tests/test-frexpl.c
index f98fcf109d..447c3b9a7a 100644
--- a/tests/test-frexpl.c
+++ b/tests/test-frexpl.c
@@ -24,6 +24,8 @@
 #include <float.h>
 #include <stdlib.h>
 
+#include "fpucw.h"
+
 #define ASSERT(expr) if (!(expr)) abort ();
 
 static long double
@@ -41,6 +43,9 @@ main ()
 {
   int i;
   long double x;
+  DECL_LONG_DOUBLE_ROUNDING
+
+  BEGIN_LONG_DOUBLE_ROUNDING ();
 
   { /* NaN.  */
     int exp = -9999;
diff --git a/tests/test-printf-frexpl.c b/tests/test-printf-frexpl.c
index 32c92e8f30..15cd91d85b 100644
--- a/tests/test-printf-frexpl.c
+++ b/tests/test-printf-frexpl.c
@@ -24,6 +24,8 @@
 #include <float.h>
 #include <stdlib.h>
 
+#include "fpucw.h"
+
 #define ASSERT(expr) if (!(expr)) abort ();
 
 static long double
@@ -41,6 +43,9 @@ main ()
 {
   int i;
   long double x;
+  DECL_LONG_DOUBLE_ROUNDING
+
+  BEGIN_LONG_DOUBLE_ROUNDING ();
 
   for (i = 1, x = 1.0L; i <= LDBL_MAX_EXP; i++, x *= 2.0L)
     {