This implementation of strtod() is perfectly accurate in the following
sense: It always chooses the nearest representable floating point
number, breaking ties by choosing the floating point number whose least
significant fraction digit is zero (the IEEE round-to-even rule).
Coupled with a truly accurate print routine, this makes it possible to
print a floating point number out as human readable text, and later
read it back in to obtain the exact same binary floating point number.
(Note that the GNU printf_fp routine is not accurate enough to do
this.)

The algorithm is roughly as follows:  We divide the number into an
integer part and a fraction part.  The integer part is converted to a
multiple- precision integer, the upper bits of which become the
floating point fraction.  If this process leaves some fraction bits
unfilled, then the rest are filled from the decimal fraction part, as
follows:  The fraction digits are repeatedly multipled by a power of
two, in base 10 (actually base 10,000, to get 4 at a time), and the
integer part of the result (in base 2) is shifted into the floating
point fraction so far.  Finally, when we have filled all bits of the
floating point fraction, we look at the remaining bits of our
conversion to decide how to round it.

How does this algorithm compare to other algorithms?  The only
comparable algorithm I know of is Clinger's algorithm (and Gay's
variant of Clinger's algorithm), since other algorithms do not attempt
to guarantee perfect accuracy.

Compared to Clinger's algorithm, this algorithm has three advantages:
* It is easy to show that it is perfectly accurate.  (The algorithm,
  that is; the implementation may have bugs!)
* The algorithm does not require a full suite of bignum arithmetic
  operations.  In particular, it avoids bignum divides.
* This algorithm is relatively insensitive to the quality of the
  underlying floating point hardware.  It ought to be able to give
  perfectly accurate results even on machines that don't round correctly,
  like Crays, since it never performs a floating point operation that
  needs to be rounded.  (This turns out to be important on, e.g.,
  the 680x0, since compilers, including GCC, do not produce correct
  code for IEEE-conforming arithmetic.)

However, this implementation has poorer performance than the AT&T
netlib strtod (based on Clinger's algorithm).  Although it avoids
bignum divides, it pays for this by performing lots of short integer
divides.  It may be possible to replace these by table lookup, but in
order to get reasonable sized tables we will have to shrink the unit of
internal computation from 4 decimal digits to just 2 or perhaps 3, so
the cost may outweigh the benefit.

For very large values, we calculate powers of ten in the worst possible
way, and performance could be substantially improved by adding a real
bignum multiply, perhaps augmented by table lookup.

Some sample performance figures.  These were taken on a DECstation
5000, using the MIPS compiler with -O2 optimization.  "atof" is the
system C library's atof(), netlib is the AT&T Netlib strtod(), and
strtod() is my own strtod().  Times are in microseconds.  Note that the
system atof() makes no guarantee of perfect accuracy.

input:	"123"		"1.23e30"	"1.23"		"1.23e-30"
atof	12.1		18.5		17.3		 18.9
netlib	 6.7		10.9		 9.0		 90.7
strtod	14.0		57.0*		33.9#		142.5#

* This would be substantially improved by using a real bignum multiply
  to do large powers of ten.
# These *might* be improved by going to a scheme using table lookup
  instead of integer divide.

Internals:

strtod() parses the argument into an array of digits and an associated
integer exponent, and passes that to _strtod().  For example, given input
1.23e-30, what _strtod() gets is input "123" with exponent -31.  (That is,
there is an implicit 0 to the left of _strtod()'s input string.)

Overall performance might be improved slightly by combining the parsing
with the actual conversion, rather than doing them separately.

The algorithm uses a good deal of information from the ANSI C <float.h>
header, and ought to run correctly on a variety of floating point formats,
including at least IEEE, Vax, and probably Cray.  (It would not be
appropriate for the IBM 370, since that uses Radix 16 internally.)

The header "strtod.h" contains machine dependent definitions needed by
strtod(), and "bignum.h" contains machine dependent definitions needed
by the bignum routines.  Most of these should not change on common 32-bit
architectures.  Performance improvements might be realized by widening
the Suint and Duint types on true 64-bit machines.

To make best use of this strtod(), your library needs to provide an
implementation of ldexp() that works by directly scribbling on the bits
of a floating point number, rather than repeatedly multiplying by powers
of two in a loop.  The algorithm also relies on ldexp() to detect some
overflow and underflow conditions, although it detects some itself.

Mike Haertel
mike@cs.uoregon.edu
