How do you print the EXACT value of a floating point number?

There's been a lot of work on printing floating-point numbers. The gold standard is to print out a decimal equivalent of minimal length such that when the decimal equivalent is read back in, you get the same floating-point number that you started with, no matter what the rounding mode is during readback. You can read about the algorithm in the excellent paper by Burger and Dybvig.

This question has a bureaucratic part and an algorithmic part. A floating point number is stored internally as (2^e × m), where e is an exponent (itself in binary) and m is a mantissa. The bureaucratic part of the question is how to access this data, but R. seems more interested in the algorithmic part of the question, namely, converting (2^e × m) to a fraction (a/b) in decimal form. The answer to the bureaucratic question in several languages is frexp (which is an interesting detail that I didn’t know before today).

It is true that at first glance, it takes O(e²) work just to write 2^e in decimal, and more time still for the mantissa. But, thanks to the magic of the Schönhage–Strassen fast multiplication algorithm, you can do it in Õ(e) time, where the tilde means “up to log factors”. If you view Schönhage–Strassen as magic, then it’s not that hard to think of what to do. If e is even, you can recursively compute 2^e/2, and then square it using fast multiplication. On the other hand if e is odd, you can recursively compute 2^e−1 and then double it. You have to be careful to check that there is a version of Schönhage–Strassen in base 10. Although it is not widely documented, it can be done in any base.

Converting a very long mantissa from binary to base 10 is not exactly the same question, but it has a similar answer. You can divide the mantissa into two halves, m = a × 2^k + b. Then recursively convert a and b to base 10, convert 2^k to base 10, and do another fast multiplication to compute m in base 10.

The abstract result behind all of this is that you can convert integers from one base to another in Õ(N) time.

If the question is about standard 64-bit floating point numbers, then it’s too small for the fancy Schönhage–Strassen algorithm. In this range you can instead save work with various tricks. One approach is to store all 2048 values of 2^e in a lookup table, and then work in the mantissa with asymmetric multiplication (in between long multiplication and short multiplication). Another trick is to work in base 10000 (or a higher power of 10, depending on architecture) instead of base 10. But, as R. points out in the comments, 128-bit floating point numbers already allow large enough exponents to call into question both lookup tables and standard long multiplication. As a practical matter, long multiplication is the fastest up to a handful of digits, then in a significant medium range one can use Karatsuba multiplication or Toom–Cook multiplication, and then after that a variation of Schönhage–Strassen is best not just in theory but also in practice.

Actually, the big integer package GMP already has Õ(N)-time radix conversion, as well as good heuristics for which choice of multiplication algorithm. The only difference between their solution and mine is that instead of doing any big arithmetic in base 10, they compute large powers of 10 in base 2. In this solution, they also need fast division, but that can be obtained from fast multiplication in any of several ways.

Although it's C# and your question is tagged with C, Jon Skeet has code to convert a double to its exact representation as a string: http://www.yoda.arachsys.com/csharp/DoubleConverter.cs

From a quick glance, it does not appear to be too hard to port to C, and even easier to write in C++.

Upon further reflection, it appears that Jon's algorithm is also O(e^2), as it also loops over the exponent. However, that means the algorithm is O(log(n)^2) (where n is the floating-point number), and I'm not sure you can convert from base 2 to base 10 in better than log-squared time.

I see you’ve accepted an answer already, but here are a couple of open source implementations of this conversion you might want to look at:

1. David Gay’s dtoa() function in dtoa.c: https://www.netlib.org/fp/dtoa.c.

2. The function ___printf_fp() in the /stdio-common/printf_fp.c file in Glibc (https://ftp.gnu.org/gnu/glibc/glibc-2.11.2.tar.gz, for example).

Both will print as many digits as you ask for in a %f-type printf, as I’ve written about at:

• https://www.exploringbinary.com/print-precision-of-dyadic-fractions-varies-by-language/ and
• https://www.exploringbinary.com/print-precision-of-floating-point-integers-varies-too/.