[Bug libfortran/48511] New: Implement Steele-White algorithm for numeric output

[Bug libfortran/48511] New: Implement Steele-White algorithm for numeric output

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

           Summary: Implement Steele-White algorithm for numeric output
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: enhancement
          Priority: P3
         Component: libfortran
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: thenlich@users.sourceforge.net


Created attachment 23923
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23923
Demonstration of correct rounding a la Steele-White

In their paper "How to Print Floating-Point Numbers Accurately" [1], Steele &
White describe an algorithm which for a internal binary floating-point number
outputs the shortest external decimal representation of that number which
converts back to the original number.

This should be implemented in libfortran, possibly as a separate I/O rounding
mode.

It would eliminate the problem of getting "odd" numbers, like 0.177699999999...
when the value 0.1977 is much shorter and represents the same exact binary
value.

Basically, their algorithm does the same as the attached program demonstrates,
only more efficiently.

[1] http://grouper.ieee.org/groups/754/email/pdfq3pavhBfih.pdf

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #1 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-04-08 15:57:55 UTC ---
I have looked at this before and we have discussed it on the list.  Ok to flag
this as an enhancement.  It would be worth investigating one kind (say kind=8)
and run some timing tests to see if we can get some performance improvement.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #2 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-04-08 16:10:30 UTC ---
See http://www.cs.indiana.edu/~burger/FP-Printing-PLDI96.pdf

For improved method.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

Janne Blomqvist <jb at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |NEW
   Last reconfirmed|                            |2011.04.09 21:19:09
                 CC|                            |jb at gcc dot gnu.org
     Ever Confirmed|0                           |1

--- Comment #3 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-04-09 21:19:09 UTC ---
Does any of the Fortran edit descriptors require, or for that matter allow,
this kind of "shortest decimal representation" output? If not, the one place
where I see this could be useful if we decide to change list formatted output
to always use the shortest field width for variables, as some other compilers
do.

However, the output part is only half of the puzzle; if we do this we must make
sure that the input routines are able to convert the shortest decimal
representation into the correct binary representation.

FWIW, rather than the Steele & White algorithm or the Burger one, most actual
uses seems to use David Gay's implementation for performance reasons. AFAIK
glibc uses code based on Gay's, and also libquadmath in turn uses code based on
glibc. So maybe we can find something there. There's also
libjava/classpath/native/fdlibm/dtoa.c .

For some discussion about issues with this kind of conversions, see e.g.
http://bugs.python.org/issue1580

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #4 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-04-10 08:36:24 UTC ---
(In reply to comment #3)
> Does any of the Fortran edit descriptors require, or for that matter allow,
> this kind of "shortest decimal representation" output?

Well, the obvious(?) answer is the various descriptors with 0 width.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #5 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-04-10 10:19:47 UTC ---
(In reply to comment #4)
> (In reply to comment #3)
> > Does any of the Fortran edit descriptors require, or for that matter allow,
> > this kind of "shortest decimal representation" output?
> 
> Well, the obvious(?) answer is the various descriptors with 0 width.

Not quite:

Since we're not talking about the shortest width (w) but the smallest number of
decimal digits in the fraction (d), only those descriptors where we can select
"a processor-dependent reasonable value for d" can be "shortened":

That would be list-directed and G0 (but not G0.d).

For all others the algorithm is still useful if we append zeros to fill the
required width: it is better to print 0.3000000 than 0.2999996

I'm still not sure where the Fortran standard says "a processor-dependent
reasonable value for d" that includes a "processor-dependent reasonable value
for d which also depends on the internal value to be printed" because that
would mess up tabular aligment for printing. This might be a feature users rely
on.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #6 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-04-10 12:24:54 UTC ---
(In reply to comment #5)
> (In reply to comment #4)
> > (In reply to comment #3)
> > > Does any of the Fortran edit descriptors require, or for that matter allow,
> > > this kind of "shortest decimal representation" output?
> > 
> > Well, the obvious(?) answer is the various descriptors with 0 width.
> 
> Not quite:
> 
> Since we're not talking about the shortest width (w) but the smallest number of
> decimal digits in the fraction (d), only those descriptors where we can select
> "a processor-dependent reasonable value for d" can be "shortened":
> 
> That would be list-directed and G0 (but not G0.d).

What I meant was edit descriptors with w=0 and where d is not present. Sorry if
that wasn't clear.

> For all others the algorithm is still useful if we append zeros to fill the
> required width: it is better to print 0.3000000 than 0.2999996

I don't this is useful, actually. If the user explicitly asks for d decimal
digits, we should use the existing simple algorithm ("simple" compared to the
shortest-width one) to round to d digits, rather than finding the shortest
representation and pad with zeros if necessary.

> I'm still not sure where the Fortran standard says "a processor-dependent
> reasonable value for d"

I don't think you'll find that anywhere. If d is present, d digits must be
printed, assuming it fits in the field width. 

> that includes a "processor-dependent reasonable value
> for d which also depends on the internal value to be printed" because that
> would mess up tabular aligment for printing. This might be a feature users rely
> on.

I think the intention with the w=0 edit descriptors is exactly that, the
compiler can choose the minimum width depending on the value (this being the
case where the shortest-width algorithm can be useful). If the user doesn't
want that, he can choose a non-zero width.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #7 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-04-27 14:11:10 UTC ---
Gay's routines (http://www.netlib.org/fp/) can handle double, float, extended,
quad; rounding directions may be specified. They can output a given number of
digits (compatible to all Fortran formats where this is required), thus
mimicking the "regular" printf behaviour; or find the shortest number of digits
that round-trip to the same value (good for G0 editing). This should be worth a
try.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #8 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-04-27 14:59:39 UTC ---
(In reply to comment #7)
> Gay's routines (http://www.netlib.org/fp/) can handle double, float, extended,
> quad

libquadmath initially used dtoa for __float128 [REAL(16)] I/O, but there were
several issues especially with NaN; thus, one now uses functions based on
GLIBC. (At least the code I saw in libquadmath/gdtoa was also a mess, maybe
using "dtoa.c" itself is better.)

In any case, one could think of using Fortran tailored I/O functions instead of
using the C ones; and Gay's routines might be a starting point. (Interestingly,
the netlib files are still maintained, the last changelog entry is from
2011-03-21!)

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #9 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-04-27 18:36:12 UTC ---
I've been thinking about this issue now and then over the past few weeks, and
funnily enough, while I initially thought this would be a nice improvement,
I've subsequently come to the opposite conclusion, namely that this shouldn't
be done:

- The use case is limited to making G0, and potentially list formatted, output
more compact in some cases. So the benefit isn't that big, IMHO.

- Even if we can start off with David Gay's code instead of reimplementing from
scratch, it means importing and maintaining potentially several thousand lines
of very tricky code.

- Getting something wrong may not only mean that we get some conversion wrong
in the last bit, but as can be seen in the recent issues with PHP and Java,
careless usage of Gay's code might lead to an infinite loop.

- As Tobias mentioned, there was enough issues with adapting Gay's code for
libquadmath that it was replaced with code from glibc.

- OTOH, currently we make use of {quadmath_}snprintf and strto{f,,ld,flt128},
which seem to be competently implemented at least on the widely used targets
(that is, they provide correctly rounded behavior), and are widely used also
outside of gfortran. This should definitely not be underestimated.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #10 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-04-28 09:30:10 UTC ---
As an addendum to the above, one thing we could do with modest effort and
without importing Gay's code would be to remove trailing zeros for G0 (and
maybe list directed) output.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #11 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-04-28 10:01:45 UTC ---
(In reply to comment #10)
> As an addendum to the above, one thing we could do with modest effort and
> without importing Gay's code would be to remove trailing zeros for G0 (and
> maybe list directed) output.

For this to be useful, it would require to choose a precision of not more than
P(bf) = floor(p×log10(2)); or floor(p×log10(2)) - 1, I'm not quite sure.

That would be P = {7, 15, 19, 34} or {6, 14, 18, 33}.

We then would have a loss of precision and no guarantee of round-trip
conversion to the same number.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #12 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-04-29 01:42:00 UTC ---
Created attachment 24137
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=24137
An implementation found via google

Found this while searching around. My impression is that this is what glibc
uses already. If my impression is true, there is no need for this PR.  If not,
maybe this is a starting point.

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

--- Comment #13 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-05-05 10:57:33 UTC ---
G95 is actually using this method for list-directed output.

print *, .3, .33, .333, .3333, .33333, .333333, .3333333
g95 =>  0.3 0.33 0.333 0.3333 0.33333 0.333333 0.3333333
gfortran => 0.300000012 0.330000013 0.333000004 0.333299994 0.333330005
0.333332986 0.333333313

From ftoa.c:


/* All software © 1996 Robert G. Burger.  Permission is hereby granted,
...
 * This code is described in an article by Robert G. Burger and
 * R. Kent Dybvig in "Proceedings of the ACM SIGPLAN '96 Conference on
 * Programming Language Design and Implementation", pages 108-116, May
 * 1996.
 *
 * Modified for use with g95 by Andy Vaught.  Thanks Bob!
 */

[Bug libfortran/48511] Implement Steele-White algorithm for numeric output

[fxcoudert at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=48511

Francois-Xavier Coudert <fxcoudert at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
                 CC|                            |fxcoudert at gcc dot gnu.org
         Resolution|---                         |WONTFIX

--- Comment #14 from Francois-Xavier Coudert <fxcoudert at gcc dot gnu.org> ---
Seeing the PR track, it seems both Janne and Jerry are against reimplementing
Steele-White. I am too: unless we can show very good features or speed
improvement, we'd best do exactly what we do now, i.e. use the system's libc
(and possibly libquadmath). So I am closing this PR as WONTFIX.


In any case, if we wanted to have faster floating-point I/O, we should move to
something more modern than Steele-White. The current "state of the art" (as
implemented in new language runtimes, such as Julia, rust, modern javascript,
etc.) is to combine grisu3
(http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf) with
dragon4 (Steele-White) as a fall-back for the cases where grisu3 doesn't round
exactly.