'float' fixes for FreeBSD, AIX, IRIX

[Bruno Haible]

> > FreeBSD 6.4:
> > 
> > test-isinf.c:151: assertion failed
> > Abort trap (core dumped)
> > FAIL: test-isinf
> 
> This test failure occurs because lib/isinf.c relies on LDBL_MAX, and on
> FreeBSD/x86, LDBL_MAX is Infinity. Ouch.

We had known about this bug, see
<http://lists.gnu.org/archive/html/bug-gnulib/2008-07/msg00118.html>,
but not added a workaround.

The bug on FreeBSD/x86 is that the compiler thinks that 'long double's
have 53 bits of precision, but they really have 64 bits of precision.
When LDBL_MAX is defined to 0x0.ffffffffffffffffp16384L, gcc converts it
to Infinity. Whereas when LDBL_MAX is defined to 0x0.fffffffffffff8p16384L,
it is possible to multiply it with (1 + 2^-54) or so, and get finite
'long double' values that are greater than LDBL_MAX. What can we do?
The only workaround that I found (and that also disables GCC optimizations)
is to define LDBL_MAX as a variable in a separate file.

On AIX 7.1 with gcc 4.2 - a platform with a "double double" representation
for 'long double' - the situation is similar: The real LDBL_MAX value is
{ 0x7FEFFFFF, 0xFFFFFFFF, 0x7C8FFFFF, 0xFFFFFFFF }, but the largest value
that GCC allows us to define through a constant expression is
{ 0x7FEFFFFF, 0xFFFFFFFF, 0x7C8FFFFF, 0xFFFFFFFE }. Here too the fix is
to define LDBL_MAX externally.

On IRIX 6.5 - another platform with a "double double" representation for
'long double' - the situation is simpler. The cc compiler gets LDBL_MANT_DIG
wrong, and gcc gets LDBL_MIN_EXP, LDBL_MIN, LDBL_EPSILON wrong.

I'm applying these fixes. Unfortunately, it attaches an AC_LIBOBJ to a
header file, and unfortunately LDBL_MAX is not a constant expression any
more on all platforms. But that's the best we can get.


2011-06-20  Bruno Haible  <address@hidden>

        float: Work around <float.h> bugs on FreeBSD/x86, AIX with GCC, IRIX.
        * lib/float.in.h: Add workarounds for FreeBSD/x86, AIX with GCC, IRIX.
        * lib/float.c: New file.
        * m4/float_h.m4 (gl_FLOAT_H): Also handle FreeBSD, AIX, IRIX. Set
        REPLACE_FLOAT_LDBL.
        * modules/float (Files): Add lib/float.c.
        (configure.ac): Invoke AC_LIBOBJ.
        * doc/posix-headers/float.texi: Mention problems on FreeBSD, AIX, IRIX.

================================= lib/float.c =================================
/* Auxiliary definitions for <float.h>.
   Copyright (C) 2011 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 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 General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */

#include <config.h>

/* Specification.  */
#include <float.h>

#if (defined _ARCH_PPC || defined _POWER) && defined _AIX && (LDBL_MANT_DIG == 
106) && defined __GNUC__
const union gl_long_double_union gl_LDBL_MAX =
  { { DBL_MAX, DBL_MAX / (double)134217728UL / (double)134217728UL } };
#elif defined __i386__
const union gl_long_double_union gl_LDBL_MAX =
  { { 0xFFFFFFFF, 0xFFFFFFFF, 32766 } };
#else
/* This declaration is solely to ensure that after preprocessing
   this file is never empty.  */
typedef int dummy;
#endif
===============================================================================
--- doc/posix-headers/float.texi.orig   Mon Jun 20 12:33:40 2011
+++ doc/posix-headers/float.texi        Mon Jun 20 12:33:23 2011
@@ -12,6 +12,10 @@
 On OpenBSD 4.0, MirBSD 10, and BeOS, they are the same as the values of the
 @code{DBL_*} macros, although @samp{long double} is a larger type than
 @samp{double}.
+On FreeBSD/x86 6.4, they represent the incorrect 53-bit precision assumptions
+in the compiler, not the real 64-bit precision at runtime.
+On AIX 7.1 with GCC 4.2 and on IRIX 6.5, they don't reflect the
+``double double'' representation of @code{long double} correctly.
 @end itemize
 
 Portability problems not fixed by Gnulib:
--- lib/float.in.h.orig Mon Jun 20 12:33:40 2011
+++ lib/float.in.h      Mon Jun 20 12:11:55 2011
@@ -29,6 +29,7 @@
 #define address@hidden@_FLOAT_H
 
 /* 'long double' properties.  */
+
 #if defined __i386__ && (defined __BEOS__ || defined __OpenBSD__)
 /* Number of mantissa units, in base FLT_RADIX.  */
 # undef LDBL_MANT_DIG
@@ -59,5 +60,115 @@
 # define LDBL_MAX_10_EXP 4932
 #endif
 
+/* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
+   precision in the compiler but 64 bits of precision at runtime.  See
+   <http://lists.gnu.org/archive/html/bug-gnulib/2008-07/msg00063.html>.  */
+#if defined __i386__ && defined __FreeBSD__
+/* Number of mantissa units, in base FLT_RADIX.  */
+# undef LDBL_MANT_DIG
+# define LDBL_MANT_DIG   64
+/* Number of decimal digits that is sufficient for representing a number.  */
+# undef LDBL_DIG
+# define LDBL_DIG        18
+/* x-1 where x is the smallest representable number > 1.  */
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 1.084202172485504434007452800869941711426e-19L /* 2^-63 
*/
+/* Minimum e such that FLT_RADIX^(e-1) is a normalized number.  */
+# undef LDBL_MIN_EXP
+# define LDBL_MIN_EXP    (-16381)
+/* Maximum e such that FLT_RADIX^(e-1) is a representable finite number.  */
+# undef LDBL_MAX_EXP
+# define LDBL_MAX_EXP    16384
+/* Minimum positive normalized number.  */
+# undef LDBL_MIN
+# define LDBL_MIN        3.3621031431120935E-4932L /* = 0x1p-16382L */
+/* Maximum representable finite number.  */
+# undef LDBL_MAX
+/* LDBL_MAX is represented as { 0xFFFFFFFF, 0xFFFFFFFF, 32766 }.
+   But the largest literal that GCC allows us to write is
+   0x0.fffffffffffff8p16384L = { 0xFFFFF800, 0xFFFFFFFF, 32766 }.
+   So, define it like this through a reference to an external variable
+
+     const unsigned int LDBL_MAX[3] = { 0xFFFFFFFF, 0xFFFFFFFF, 32766 };
+     extern const long double LDBL_MAX;
+
+   Unfortunately, this is not a constant expression.  */
+union gl_long_double_union
+  {
+    struct { unsigned int lo; unsigned int hi; unsigned int exponent; } xd;
+    long double ld;
+  };
+extern const union gl_long_double_union gl_LDBL_MAX;
+# define LDBL_MAX (gl_LDBL_MAX.ld)
+/* Minimum e such that 10^e is in the range of normalized numbers.  */
+# undef LDBL_MIN_10_EXP
+# define LDBL_MIN_10_EXP (-4931)
+/* Maximum e such that 10^e is in the range of representable finite numbers.  
*/
+# undef LDBL_MAX_10_EXP
+# define LDBL_MAX_10_EXP 4932
+#endif
+
+/* On AIX 7.1 with gcc 4.2, the values of LDBL_MIN_EXP, LDBL_MIN, LDBL_MAX are
+   wrong.  */
+#if (defined _ARCH_PPC || defined _POWER) && defined _AIX && (LDBL_MANT_DIG == 
106) && defined __GNUC__
+# undef LDBL_MIN_EXP
+# define LDBL_MIN_EXP DBL_MIN_EXP
+# undef LDBL_MIN_10_EXP
+# define LDBL_MIN_10_EXP DBL_MIN_10_EXP
+# undef LDBL_MIN
+# define LDBL_MIN 2.22507385850720138309023271733240406422e-308L /* DBL_MIN = 
2^-1022 */
+# undef LDBL_MAX
+/* LDBL_MAX is represented as { 0x7FEFFFFF, 0xFFFFFFFF, 0x7C8FFFFF, 0xFFFFFFFF 
}.
+   It is not easy to define:
+     #define LDBL_MAX 1.79769313486231580793728971405302307166e308L
+   is too small, whereas
+     #define LDBL_MAX 1.79769313486231580793728971405302307167e308L
+   is too large.  Apparently a bug in GCC decimal-to-binary conversion.
+   Also, I can't get values larger than
+     #define LDBL63 ((long double) (1ULL << 63))
+     #define LDBL882 (LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * 
LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63)
+     #define LDBL945 (LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * 
LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63)
+     #define LDBL1008 (LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * 
LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 * LDBL63 
* LDBL63)
+     #define LDBL_MAX (LDBL1008 * 65535.0L + LDBL945 * (long double) 
9223372036821221375ULL + LDBL882 * (long double) 4611686018427387904ULL)
+   which is represented as { 0x7FEFFFFF, 0xFFFFFFFF, 0x7C8FFFFF, 0xF8000000 }.
+   So, define it like this through a reference to an external variable
+
+     const double LDBL_MAX[2] = { DBL_MAX, DBL_MAX / (double)134217728UL / 
(double)134217728UL };
+     extern const long double LDBL_MAX;
+
+   or through a pointer cast
+
+     #define LDBL_MAX \
+       (*(const long double *) (double[]) { DBL_MAX, DBL_MAX / 
(double)134217728UL / (double)134217728UL })
+
+   Unfortunately, this is not a constant expression, and the latter expression
+   does not work well when GCC is optimizing..  */
+union gl_long_double_union
+  {
+    struct { double hi; double lo; } dd;
+    long double ld;
+  };
+extern const union gl_long_double_union gl_LDBL_MAX;
+# define LDBL_MAX (gl_LDBL_MAX.ld)
+#endif
+
+/* On IRIX 6.5, with cc, the value of LDBL_MANT_DIG is wrong.
+   On IRIX 6.5, with gcc 4.2, the values of LDBL_MIN_EXP, LDBL_MIN, 
LDBL_EPSILON
+   are wrong.  */
+#if defined __sgi && (LDBL_MANT_DIG >= 106)
+# undef LDBL_MANT_DIG
+# define LDBL_MANT_DIG 106
+# if defined __GNUC__
+#  undef LDBL_MIN_EXP
+#  define LDBL_MIN_EXP DBL_MIN_EXP
+#  undef LDBL_MIN_10_EXP
+#  define LDBL_MIN_10_EXP DBL_MIN_10_EXP
+#  undef LDBL_MIN
+#  define LDBL_MIN 2.22507385850720138309023271733240406422e-308L /* DBL_MIN = 
2^-1022 */
+#  undef LDBL_EPSILON
+#  define LDBL_EPSILON 2.46519032881566189191165176650870696773e-32L /* 2^-105 
*/
+# endif
+#endif
+
 #endif /* address@hidden@_FLOAT_H */
 #endif /* address@hidden@_FLOAT_H */
--- m4/float_h.m4.orig  Mon Jun 20 12:33:40 2011
+++ m4/float_h.m4       Mon Jun 20 11:41:05 2011
@@ -1,4 +1,4 @@
-# float_h.m4 serial 6
+# float_h.m4 serial 7
 dnl Copyright (C) 2007, 2009-2011 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -9,12 +9,41 @@
   AC_REQUIRE([AC_PROG_CC])
   AC_REQUIRE([AC_CANONICAL_HOST])
   FLOAT_H=
+  REPLACE_FLOAT_LDBL=0
   case "$host_os" in
-    beos* | openbsd* | mirbsd*)
+    aix* | beos* | openbsd* | mirbsd* | irix*)
       FLOAT_H=float.h
-      gl_NEXT_HEADERS([float.h])
+      ;;
+    freebsd*)
+      case "$host_cpu" in
+changequote(,)dnl
+        i[34567]86 )
+changequote([,])dnl
+          FLOAT_H=float.h
+          ;;
+        x86_64 )
+          # On x86_64 systems, the C compiler may still be generating
+          # 32-bit code.
+          AC_EGREP_CPP([yes],
+            [#if defined __LP64__ || defined __x86_64__ || defined __amd64__
+             yes
+             #endif],
+            [],
+            [FLOAT_H=float.h])
+          ;;
+      esac
+      ;;
+  esac
+  case "$host_os" in
+    aix* | freebsd*)
+      if test -n "$FLOAT_H"; then
+        REPLACE_FLOAT_LDBL=1
+      fi
       ;;
   esac
+  if test -n "$FLOAT_H"; then
+    gl_NEXT_HEADERS([float.h])
+  fi
   AC_SUBST([FLOAT_H])
   AM_CONDITIONAL([GL_GENERATE_FLOAT_H], [test -n "$FLOAT_H"])
 ])
--- modules/float.orig  Mon Jun 20 12:33:40 2011
+++ modules/float       Mon Jun 20 02:23:42 2011
@@ -3,6 +3,7 @@
 
 Files:
 lib/float.in.h
+lib/float.c
 m4/float_h.m4
 
 Depends-on:
@@ -10,6 +11,9 @@
 
 configure.ac:
 gl_FLOAT_H
+if test $REPLACE_FLOAT_LDBL = 1; then
+  AC_LIBOBJ([float])
+fi
 
 Makefile.am:
 BUILT_SOURCES += $(FLOAT_H)

-- 
In memoriam Neda Agha-Soltan <http://en.wikipedia.org/wiki/Neda_Agha-Soltan>