[PATCH] Fix {,q}ecvt{,_r} on denormals

[Jakub Jelinek <jakub at redhat dot com>]

Hi!

As shown by the testcase below, ecvt/ecvt_r doesn't work on denormals
(well, it works for double denormals on IA-32 unless -mfpmath=sse,
as there if you are lucky f is computed in long double precision and
not spilled until done with it).  The problem is that a denormal
needs to be multiplied by a power of 10 bigger than what can fit into
the floating type (double resp. long double) if the result of the
multiplication should be bigger or equal than 1.0.
So current CVS code overflows f into +Inf and efcvt_r gives incorrect
result.

If ecvt/ecvt_r was some function where we would care a lot about speed,
it would be perhaps good to also use some precomputed powers
of ten (1e-4096, 1e-2048, ... 1e1, 1e2, 1e4, 1e8, 1e16, 1e32, 1e64, ... 1e4096),
so that we don't for close to smallest normalized doubles
multiply f by 10 ~ 300 times and for really small long doubles ~ 4900
times.  But given that it is a legacy function, I'm not sure if it
is ok to increase .rodata by 16 * sizeof (double) and 24 * sizeof (long
double) and code size because of that.  What do you think?

2004-12-21  Jakub Jelinek  <jakub@redhat.com>

	* misc/efgcvt_r.c (FLOAT_MIN_10_EXP, FLOAT_MIN_10_NORM): Define.
	(ecvt_r): Special case denormals.
	* misc/qefgcvt_r.c (FLOAT_MIN_10_EXP, FLOAT_MIN_10_NORM): Define.
	* misc/tst-efgcvt.c: Include float.h.
	(ecvt_tests): Add 2 new tests.

--- libc/misc/tst-efgcvt.c.jj	2001-07-06 06:55:36.000000000 +0200
+++ libc/misc/tst-efgcvt.c	2004-12-21 18:41:02.890224202 +0100
@@ -1,4 +1,4 @@
-/* Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
+/* Copyright (C) 1998, 1999, 2000, 2004 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -20,6 +20,7 @@
 # define _GNU_SOURCE	1
 #endif
 
+#include <float.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -59,6 +60,10 @@ static testcase ecvt_tests[] =
   { 123.01, -4, 3, "" },
   { 126.71, -4, 3, "" },
   { 0.0, 4, 1, "0000" },
+#if DBL_MANT_DIG == 53
+  { 0x1p-1074, 3, -323, "494" },
+  { -0x1p-1074, 3, -323, "494" },
+#endif
   /* -1.0 is end marker.  */
   { -1.0, 0, 0, "" }
 };
--- libc/misc/qefgcvt_r.c.jj	2004-05-07 14:32:16.000000000 +0200
+++ libc/misc/qefgcvt_r.c	2004-12-21 18:18:31.261313136 +0100
@@ -24,6 +24,7 @@
 #define FUNC_PREFIX q
 #define FLOAT_FMT_FLAG "L"
 #define FLOAT_NAME_EXT l
+#define FLOAT_MIN_10_EXP LDBL_MIN_10_EXP
 #if LDBL_MANT_DIG == 64
 # define NDIGIT_MAX 21
 #elif LDBL_MANT_DIG == 53
@@ -40,5 +41,16 @@
 # error "NDIGIT_MAX must be precomputed"
 # define NDIGIT_MAX (lrint (ceil (M_LN2 / M_LN10 * LDBL_MANT_DIG + 1.0)))
 #endif
+#if LDBL_MIN_10_EXP == -37
+# define FLOAT_MIN_10_NORM	1.0e-37L
+#elif LDBL_MIN_10_EXP == -307
+# define FLOAT_MIN_10_NORM	1.0e-307L
+#elif LDBL_MIN_10_EXP == -4931
+# define FLOAT_MIN_10_NORM	1.0e-4931L
+#else
+/* libc can't depend on libm.  */
+# error "FLOAT_MIN_10_NORM must be precomputed"
+# define FLOAT_MIN_10_NORM	exp10l (LDBL_MIN_10_EXP)
+#endif
 
 #include "efgcvt_r.c"
--- libc/misc/efgcvt_r.c.jj	2004-08-04 14:42:22.000000000 +0200
+++ libc/misc/efgcvt_r.c	2004-12-21 18:31:52.486396147 +0100
@@ -1,5 +1,5 @@
 /* Compatibility functions for floating point formatting, reentrant versions.
-   Copyright (C) 1995,96,97,98,99,2000,01,02 Free Software Foundation, Inc.
+   Copyright (C) 1995,96,97,98,99,2000,01,02,04 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -31,6 +31,7 @@
 # define FUNC_PREFIX
 # define FLOAT_FMT_FLAG
 # define FLOAT_NAME_EXT
+# define FLOAT_MIN_10_EXP DBL_MIN_10_EXP
 # if DBL_MANT_DIG == 53
 #  define NDIGIT_MAX 17
 # elif DBL_MANT_DIG == 24
@@ -43,6 +44,17 @@
 #  error "NDIGIT_MAX must be precomputed"
 #  define NDIGIT_MAX (lrint (ceil (M_LN2 / M_LN10 * DBL_MANT_DIG + 1.0)))
 # endif
+# if DBL_MIN_10_EXP == -37
+#  define FLOAT_MIN_10_NORM	1.0e-37
+# elif DBL_MIN_10_EXP == -307
+#  define FLOAT_MIN_10_NORM	1.0e-307
+# elif DBL_MIN_10_EXP == -4931
+#  define FLOAT_MIN_10_NORM	1.0e-4931
+# else
+/* libc can't depend on libm.  */
+#  error "FLOAT_MIN_10_NORM must be precomputed"
+#  define FLOAT_MIN_10_NORM	exp10 (DBL_MIN_10_EXP)
+# endif
 #endif
 
 #define APPEND(a, b) APPEND2 (a, b)
@@ -171,6 +183,17 @@ APPEND (FUNC_PREFIX, ecvt_r) (value, ndi
 	d = -value;
       else
 	d = value;
+      /* For denormalized numbers the d < 1.0 case below won't work,
+	 as f can overflow to +Inf.  */
+      if (d < FLOAT_MIN_10_NORM)
+	{
+	  value /= FLOAT_MIN_10_NORM;
+	  if (value < 0.0)
+	    d = -value;
+	  else
+	    d = value;
+	  exponent += FLOAT_MIN_10_EXP;
+	}
       if (d < 1.0)
 	{
 	  do

	Jakub