NetBSD. This platform not only lacks all 'long double' math functions, starting with frexpl() and ldexpl().
It also initializes the x86 FPU control word to a value that causes all floating-point operations to round to 'double' precision (53 mantissa bits). Nearly all 'long double' operations therefore lead to wrong results. To take a simple example, 1.01L * 0.5L evaluates to a value != 0.505L In theory the compiler should arrange to set the FPU control word at the beginning of a function that uses 'long double'. But gcc doesn't do it, and I cannot find a gcc option or __attribute__ that would do it. So we have to do it explicitly: Every function that does 'long float' computations needs to enable a control word that supports it. 2007-03-24 Bruno Haible <[EMAIL PROTECTED]> * 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. ========================= modules/fpucw ================================= 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 ========================== lib/fpucw.h ================================== /* Manipulating the FPU control word. Copyright (C) 2007 Free Software Foundation, Inc. Written by Bruno Haible <[EMAIL PROTECTED]>, 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 */ ========================================================================= *** lib/frexp.c 22 Mar 2007 11:52:34 -0000 1.2 --- lib/frexp.c 25 Mar 2007 02:26:41 -0000 *************** *** 28,33 **** --- 28,34 ---- # include <float.h> # ifdef USE_LONG_DOUBLE # include "isnanl-nolibm.h" + # include "fpucw.h" # else # include "isnan.h" # endif *************** *** 40,50 **** --- 41,57 ---- # 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,58 **** --- 60,66 ---- { int sign; int exponent; + DECL_ROUNDING /* Test for NaN, infinity, and zero. */ if (ISNAN (x) || x + x == x) *************** *** 68,73 **** --- 76,83 ---- 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,156 **** /* Here 0.5 <= x < 1.0. */ } *exp = exponent; ! return (sign < 0 ? - x : x); } #else --- 159,171 ---- /* Here 0.5 <= x < 1.0. */ } + if (sign < 0) + x = - x; + + END_ROUNDING (); + *exp = exponent; ! return x; } #else *** lib/printf-frexp.c 7 Mar 2007 01:39:07 -0000 1.5 --- lib/printf-frexp.c 25 Mar 2007 02:26:41 -0000 *************** *** 28,33 **** --- 28,36 ---- # 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,47 **** --- 45,53 ---- # 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,57 **** --- 58,66 ---- # define FREXP frexp # define LDEXP ldexp # endif + # define DECL_ROUNDING + # define BEGIN_ROUNDING() + # define END_ROUNDING() # define L_(literal) literal # endif *************** *** 59,64 **** --- 68,76 ---- 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,175 **** --- 182,189 ---- or 1.0 <= x < 2.0 and exponent >= MIN_EXP - 1. */ # endif + END_ROUNDING (); + *exp = exponent; return x; } *** lib/vasnprintf.c 25 Mar 2007 01:11:49 -0000 1.34 --- lib/vasnprintf.c 25 Mar 2007 02:26:42 -0000 *************** *** 57,62 **** --- 57,63 ---- # if HAVE_LONG_DOUBLE # include "isnanl-nolibm.h" # include "printf-frexpl.h" + # include "fpucw.h" # endif #endif *************** *** 415,420 **** --- 416,424 ---- else { int sign = 0; + DECL_LONG_DOUBLE_ROUNDING + + BEGIN_LONG_DOUBLE_ROUNDING (); if (arg < 0.0L) { *************** *** 542,547 **** --- 546,553 ---- while (*p != '\0') p++; } + + END_LONG_DOUBLE_ROUNDING (); } } else *** modules/fprintf-posix 9 Mar 2007 02:59:39 -0000 1.1 --- modules/fprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 14,19 **** --- 14,20 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_FPRINTF_POSIX *** modules/frexpl 23 Mar 2007 01:26:24 -0000 1.1 --- modules/frexpl 25 Mar 2007 02:26:42 -0000 *************** *** 9,14 **** --- 9,15 ---- Depends-on: math isnanl-nolibm + fpucw configure.ac: gl_FUNC_FREXPL *** modules/frexpl-tests 23 Mar 2007 01:27:57 -0000 1.1 --- modules/frexpl-tests 25 Mar 2007 02:26:42 -0000 *************** *** 2,7 **** --- 2,8 ---- tests/test-frexpl.c Depends-on: + fpucw configure.ac: *** modules/printf-frexpl 7 Mar 2007 01:14:04 -0000 1.3 --- modules/printf-frexpl 25 Mar 2007 02:26:42 -0000 *************** *** 11,16 **** --- 11,17 ---- Depends-on: math + fpucw configure.ac: gl_FUNC_PRINTF_FREXPL *** modules/printf-frexpl-tests 25 Feb 2007 14:29:34 -0000 1.1 --- modules/printf-frexpl-tests 25 Mar 2007 02:26:42 -0000 *************** *** 2,7 **** --- 2,8 ---- tests/test-printf-frexpl.c Depends-on: + fpucw configure.ac: *** modules/snprintf-posix 7 Mar 2007 01:59:05 -0000 1.1 --- modules/snprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 13,18 **** --- 13,19 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_SNPRINTF_POSIX *** modules/sprintf-posix 7 Mar 2007 03:47:50 -0000 1.1 --- modules/sprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 13,18 **** --- 13,19 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_SPRINTF_POSIX *** modules/vasnprintf-posix 4 Mar 2007 23:28:59 -0000 1.1 --- modules/vasnprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 12,17 **** --- 12,18 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_VASNPRINTF_POSIX *** modules/vasprintf-posix 5 Mar 2007 00:37:52 -0000 1.1 --- modules/vasprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 12,17 **** --- 12,18 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_VASPRINTF_POSIX *** modules/vfprintf-posix 9 Mar 2007 02:40:14 -0000 1.1 --- modules/vfprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 14,19 **** --- 14,20 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_VFPRINTF_POSIX *** modules/vsnprintf-posix 7 Mar 2007 01:51:44 -0000 1.1 --- modules/vsnprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 13,18 **** --- 13,19 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_VSNPRINTF_POSIX *** modules/vsprintf-posix 7 Mar 2007 03:27:10 -0000 1.1 --- modules/vsprintf-posix 25 Mar 2007 02:26:42 -0000 *************** *** 13,18 **** --- 13,19 ---- isnanl-nolibm printf-frexp printf-frexpl + fpucw configure.ac: gl_FUNC_VSPRINTF_POSIX *** tests/test-frexpl.c 23 Mar 2007 01:27:57 -0000 1.1 --- tests/test-frexpl.c 25 Mar 2007 02:26:42 -0000 *************** *** 24,29 **** --- 24,31 ---- #include <float.h> #include <stdlib.h> + #include "fpucw.h" + #define ASSERT(expr) if (!(expr)) abort (); static long double *************** *** 41,46 **** --- 43,51 ---- { int i; long double x; + DECL_LONG_DOUBLE_ROUNDING + + BEGIN_LONG_DOUBLE_ROUNDING (); { /* NaN. */ int exp = -9999; *** tests/test-printf-frexpl.c 25 Feb 2007 14:29:34 -0000 1.1 --- tests/test-printf-frexpl.c 25 Mar 2007 02:26:42 -0000 *************** *** 24,29 **** --- 24,31 ---- #include <float.h> #include <stdlib.h> + #include "fpucw.h" + #define ASSERT(expr) if (!(expr)) abort (); static long double *************** *** 41,46 **** --- 43,51 ---- { 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) {