RFA: Support infinity, NaN, and denormalized numbers in floatformat.c

Ian Lance Taylor
2003-12-03T18:12:56Z

Here is a patch to improve the support for infinity, NaN, and
denormalized numbers in floatformat.c. It's not clear how to handle
these for non-IEEE formats. But the old code did entirely the wrong
thing even when using IEEE formats, which are used by pretty much all
modern processors. So I think we might as well do the right thing
here.

These patches cause the following compilation warnings on a Red Hat
Linux 7.3 system:

../../gcc/libiberty/floatformat.c: In function `floatformat_to_double':
../../gcc/libiberty/floatformat.c:319: warning: traditional C rejects the 'u' suffix
../../gcc/libiberty/floatformat.c:319: warning: traditional C rejects initialization of unions
../../gcc/libiberty/floatformat.c:321: warning: traditional C rejects the 'f' suffix

These are all cases in which the code is using macros which are
defined by the standard <math.h> header (namely INFINITY and NAN). So
these warnings do not indicate actual code problems. I don't see any
obvious way to avoid them.

After writing these patches, I decided to see what uses the
floatformat routines. I discovered only one use: the m68k
disassembler. gdb uses different routines which probably started as a
copy of the floatformat routines, but use the type DOUBLEST which can
be `long double' instead of `double'. The code has diverged, and
there is some support there for denormalized numbers. It doesn't look
right to me, but I haven't investigated much further.

It seems to me that the gdb code should move into libiberty. We can
easily add floatformat_to_long_double and
floatformat_from_long_double, which would only be available if the
long double type is available. gdb would then call the appropriate
function from libiberty. There is other gdb floatformat code which
could move over as well, such as floatformat_is_negative() and
floatformat_is_nan().

Any other opinions on the best way to merge the libiberty and gdb
implementations?

Ian

2003-12-03 Ian Lance Taylor <ian@wasabisystems.com>

* floatformat.c: Include "config.h" and <string.h> if available.
(INFINITY, NAN): Define if not defined by <math.h>.
(floatformat_to_double): Handle NaN, infinity, and denormalized
numbers.
(floatformat_from_double): Likewise.
(ieee_test): In debugging code, use little endian rather than big
endian. Correct tests to handle NaN and to check correct sign of
zero. Omit m68k extended test.
(main): Add more debugging cases.

Index: floatformat.c
===================================================================
RCS file: /cvs/gcc/gcc/libiberty/floatformat.c,v
retrieving revision 1.14
diff -u -p -r1.14 floatformat.c
--- floatformat.c 31 Oct 2003 05:29:37 -0000 1.14
+++ floatformat.c 3 Dec 2003 18:05:31 -0000
@@ -17,16 +17,33 @@ You should have received a copy of the G
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */

+/* This is needed to pick up the NAN macro on some systems. */
+#define _GNU_SOURCE
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <math.h>
+
+#ifdef HAVE_STRING_H
+#include <string.h>
+#endif
+
#include "ansidecl.h"
+#include "libiberty.h"
#include "floatformat.h"
-#include <math.h> /* ldexp */
-#ifdef ANSI_PROTOTYPES
-#include <stddef.h>
-extern void *memcpy (void *s1, const void *s2, size_t n);
-extern void *memset (void *s, int c, size_t n);
+
+#ifndef INFINITY
+#ifdef HUGE_VAL
+#define INFINITY HUGE_VAL
#else
-extern char *memcpy ();
-extern char *memset ();
+#define INFINITY (1.0 / 0.0)
+#endif
+#endif
+
+#ifndef NAN
+#define NAN (0.0 / 0.0)
#endif

static unsigned long get_field PARAMS ((const unsigned char *,
@@ -271,9 +288,45 @@ floatformat_to_double (fmt, from, to)

exponent = get_field (ufrom, fmt->byteorder, fmt->totalsize,
fmt->exp_start, fmt->exp_len);
- /* Note that if exponent indicates a NaN, we can't really do anything useful
- (not knowing if the host has NaN's, or how to build one). So it will
- end up as an infinity or something close; that is OK. */
+
+ /* If the exponent indicates a NaN, we don't have information to
+ decide what to do. So we handle it like IEEE, except that we
+ don't try to preserve the type of NaN. FIXME. */
+ if ((unsigned long) exponent == fmt->exp_nan)
+ {
+ int nan;
+
+ mant_off = fmt->man_start;
+ mant_bits_left = fmt->man_len;
+ nan = 0;
+ while (mant_bits_left > 0)
+ {
+ mant_bits = min (mant_bits_left, 32);
+
+ if (get_field (ufrom, fmt->byteorder, fmt->totalsize,
+ mant_off, mant_bits) != 0)
+ {
+ /* This is a NaN. */
+ nan = 1;
+ break;
+ }
+
+ mant_off += mant_bits;
+ mant_bits_left -= mant_bits;
+ }
+
+ if (nan)
+ dto = NAN;
+ else
+ dto = INFINITY;
+
+ if (get_field (ufrom, fmt->byteorder, fmt->totalsize, fmt->sign_start, 1))
+ dto = -dto;
+
+ *to = dto;
+
+ return;
+ }

mant_bits_left = fmt->man_len;
mant_off = fmt->man_start;
@@ -306,8 +359,18 @@ floatformat_to_double (fmt, from, to)
mant = get_field (ufrom, fmt->byteorder, fmt->totalsize,
mant_off, mant_bits);

- dto += ldexp ((double)mant, exponent - mant_bits);
- exponent -= mant_bits;
+ /* Handle denormalized numbers. FIXME: What should we do for
+ non-IEEE formats? */
+ if (exponent == 0 && mant != 0)
+ dto += ldexp ((double)mant,
+ (- fmt->exp_bias
+ - mant_bits
+ - (mant_off - fmt->man_start)
+ + 1));
+ else
+ dto += ldexp ((double)mant, exponent - mant_bits);
+ if (exponent != 0)
+ exponent -= mant_bits;
mant_off += mant_bits;
mant_bits_left -= mant_bits;
}
@@ -392,33 +455,54 @@ floatformat_from_double (fmt, from, to)
int mant_bits_left;
unsigned char *uto = (unsigned char *)to;

- memcpy (&dfrom, from, sizeof (dfrom));
+ dfrom = *from;
memset (uto, 0, fmt->totalsize / FLOATFORMAT_CHAR_BIT);
+
+ /* If negative, set the sign bit. */
+ if (dfrom < 0)
+ {
+ put_field (uto, fmt->byteorder, fmt->totalsize, fmt->sign_start, 1, 1);
+ dfrom = -dfrom;
+ }
+
if (dfrom == 0)
- return; /* Result is zero */
+ {
+ /* 0.0. */
+ return;
+ }
+
if (dfrom != dfrom)
{
- /* From is NaN */
+ /* NaN. */
put_field (uto, fmt->byteorder, fmt->totalsize, fmt->exp_start,
fmt->exp_len, fmt->exp_nan);
- /* Be sure it's not infinity, but NaN value is irrel */
+ /* Be sure it's not infinity, but NaN value is irrelevant. */
put_field (uto, fmt->byteorder, fmt->totalsize, fmt->man_start,
32, 1);
return;
}

- /* If negative, set the sign bit. */
- if (dfrom < 0)
+ if (dfrom + dfrom == dfrom)
{
- put_field (uto, fmt->byteorder, fmt->totalsize, fmt->sign_start, 1, 1);
- dfrom = -dfrom;
+ /* This can only happen for an infinite value (or zero, which we
+ already handled above). */
+ put_field (uto, fmt->byteorder, fmt->totalsize, fmt->exp_start,
+ fmt->exp_len, fmt->exp_nan);
+ return;
}

- /* How to tell an infinity from an ordinary number? FIXME-someday */
-
mant = frexp (dfrom, &exponent);
- put_field (uto, fmt->byteorder, fmt->totalsize, fmt->exp_start, fmt->exp_len,
- exponent + fmt->exp_bias - 1);
+ if (exponent + fmt->exp_bias - 1 > 0)
+ put_field (uto, fmt->byteorder, fmt->totalsize, fmt->exp_start,
+ fmt->exp_len, exponent + fmt->exp_bias - 1);
+ else
+ {
+ /* Handle a denormalized number. FIXME: What should we do for
+ non-IEEE formats? */
+ put_field (uto, fmt->byteorder, fmt->totalsize, fmt->exp_start,
+ fmt->exp_len, 0);
+ mant = ldexp (mant, exponent + fmt->exp_bias - 1);
+ }

mant_bits_left = fmt->man_len;
mant_off = fmt->man_start;
@@ -431,12 +515,11 @@ floatformat_from_double (fmt, from, to)
mant_long = (unsigned long)mant;
mant -= mant_long;

- /* If the integer bit is implicit, then we need to discard it.
- If we are discarding a zero, we should be (but are not) creating
- a denormalized number which means adjusting the exponent
- (I think). */
+ /* If the integer bit is implicit, and we are not creating a
+ denormalized number, then we need to discard it. */
if ((unsigned int) mant_bits_left == fmt->man_len
- && fmt->intbit == floatformat_intbit_no)
+ && fmt->intbit == floatformat_intbit_no
+ && exponent + fmt->exp_bias - 1 > 0)
{
mant_long &= 0x7fffffff;
mant_bits -= 1;
@@ -468,6 +551,8 @@ floatformat_is_valid (fmt, from)

#ifdef IEEE_DEBUG

+#include <stdio.h>
+
/* This is to be run on a host which uses IEEE floating point. */

void
@@ -475,19 +560,31 @@ ieee_test (n)
double n;
{
double result;
- char exten[16];

- floatformat_to_double (&floatformat_ieee_double_big, &n, &result);
- if (n != result)
+ floatformat_to_double (&floatformat_ieee_double_little, (char *) &n,
+ &result);
+ if ((n != result && (! isnan (n) || ! isnan (result)))
+ || (n < 0 && result >= 0)
+ || (n >= 0 && result < 0))
printf ("Differ(to): %.20g -> %.20g\n", n, result);
- floatformat_from_double (&floatformat_ieee_double_big, &n, &result);
- if (n != result)
+
+ floatformat_from_double (&floatformat_ieee_double_little, &n,
+ (char *) &result);
+ if ((n != result && (! isnan (n) || ! isnan (result)))
+ || (n < 0 && result >= 0)
+ || (n >= 0 && result < 0))
printf ("Differ(from): %.20g -> %.20g\n", n, result);

- floatformat_from_double (&floatformat_m68881_ext, &n, exten);
- floatformat_to_double (&floatformat_m68881_ext, exten, &result);
- if (n != result)
- printf ("Differ(to+from): %.20g -> %.20g\n", n, result);
+#if 0
+ {
+ char exten[16];
+
+ floatformat_from_double (&floatformat_m68881_ext, &n, exten);
+ floatformat_to_double (&floatformat_m68881_ext, exten, &result);
+ if (n != result)
+ printf ("Differ(to+from): %.20g -> %.20g\n", n, result);
+ }
+#endif

#if IEEE_DEBUG > 1
/* This is to be run on a host which uses 68881 format. */
@@ -502,12 +599,22 @@ ieee_test (n)
int
main ()
{
+ ieee_test (0.0);
ieee_test (0.5);
ieee_test (256.0);
ieee_test (0.12345);
ieee_test (234235.78907234);
ieee_test (-512.0);
ieee_test (-0.004321);
+ ieee_test (1.2E-70);
+ ieee_test (1.2E-316);
+ ieee_test (4.9406564584124654E-324);
+ ieee_test (- 4.9406564584124654E-324);
+ ieee_test (- 0.0);
+ ieee_test (- INFINITY);
+ ieee_test (- NAN);
+ ieee_test (INFINITY);
+ ieee_test (NAN);
return 0;
}
#endif