Floating point to decimal conversion is not so easy

Tagged:  •    •    •    •

Russ Cox of Plan-9 and Go fame posted a blog entry titled Floating Point to Decimal Conversion is Easy. While he is usually right, I believe this time he isn't.

Floating point to decimal conversion is easy if you are okay with ugly results. A good conversion routine will print the shortest decimal representation of the floating-point number, that is, the shortest decimal number whose closest floating-point representation equals the original number. You do not want 0.30000000001, you want 0.3, because the number right above 0.3 is 0.30000000003 and 0.30000000001 does not provide any extra precision.

Otherwise, your users will complain. (And you need to make sure you got it right, otherwise they will complain even more).

Here is how GNU Smalltalk does it. I'm pretty sure it is correct, too.

Unlike the Go example in Russ Cox's article, GNU Smalltalk starts from an exact rational representation and relies on LargeIntegers. This is not a very efficient algorithm, but it is correct and optimal. Throughout the code, num will contain the error between the number to be printed (self, supposed positive) and the current partial representation (digits / weight; these variables are defined below.

num := self asExactFraction.

The algorithm tracks the decimal representation of two adjacent floating point numbers, num and the immediately successive number:

"Smallest number such that self + eps ~= eps"
eps := 2 raisedToInteger: self exponent - self class precision + 1.

Note that eps is either an Integer, possibly large, or a Fraction. self exponent is the base-2 exponent of self, self class precision is the number of digits in the mantissa (53 for doubles).

The initial approximation is, well, 0. But the denominator, weight is computed from the beginning to be close to self:

digits := 0.
exponent := num floorLog: 10.
weight := 10 raisedToInteger: exponent.

As we add digits to the approximation, the denominator will be reduced while remaining a power of 10. exponent and weight could be computed with a table or with binary search.

We run the decimal conversion until we find a different digit in num versus each of num - eps and num + eps. Along the way, we remember whether we found a digit that is not 9:

allNines := true.
sameDown := true.
sameUp := true.

[digit := num // weight.
allNines := allNines and: [digit = 9].
sameDown := sameDown and: [(num - eps) // weight = digit].
sameUp := sameUp and: [(num + eps) // weight = digit].
num := num - (digit * weight).
prevWeight := weight.
weight := weight / 10.
sameDown or: [sameUp]] whileTrue.

Now we have a correct approximation, but perhaps not an optimal one. For simplicity, I'll accept a possible unoptimality in case of floating-point numbers that are integer and so big that num+0.5 cannot be represented exactly:

eps isInteger ifTrue: [eps := eps / 2].

Without this line, round-to-even behavior of the decimal-to-binary conversion may cause bugs.

At this point num is the error that remains in the decimal representation. So, the error by tweaking the lowest significant digits in the decimal representation will be something like num - (NN * prevWeight) and num + (NN * prevWeight), respectively if you add or subtract NN from the lowest significant digit in the decimal representations.

adjust will contain what I wrote above as NN:

The first improvement is to try rounding the last digit while not changing the meaning. Try rounding down:

(digit <= 5 and: [num + (digit * prevWeight) < (eps / 2)])
ifTrue: [adjust := digit negated].

so that e.g. 0.30000000001 will set adjust to -1 ("subtract 1", giving 0.3); and then up:

(digit > 5 and: [num + ((digit - 10) * prevWeight) > (eps / -2)])
ifTrue: [adjust := 10 - digit].

In this case 0.8999999996 will set adjust to 4 ("add 4", giving 0.9).

The second improvement is only done if the above didn't trigger. It tries adding 1, and see if that moves us closer to self. This is needed because digits are always found with //, which truncates the result rather than rounding it:

(adjust = 0 and: [digit > 0]) ifTrue: [
(num - prevWeight) abs < num ifTrue: [adjust := 1]].

The tweak is also skipped if the last digit is 0. In this case, the decimal representation we have may make the representation more precise but it would also make the result longer.

Now we can perform the adjustment (but be careful about moving the decimal point, if necessary, when rounding up!):

digits := digits + adjust.
(adjust > 0 and: [allNines])
ifTrue: [allNines := false. exponent := exponent + 1].

Now we have converted self asExactFraction (a fraction from one whose denominator is a power of 2) to another fraction, whose denominator is a power of 10. digits is the numerator of this new fraction, while the denominator is 10exponent. Let's convert the numerator to a string in order to print it:

digits := digits printString.

Keep the significant digits only, we convert it to a stream for easy iteration:

precision := digits findLast: [:ch | ch ~= \$0].
on: digits
from: 1
to: precision.

GNU Smalltalk now chooses whether to print in exponential notation or not. For simplicity I will cover the exponential notation case only.

The destination stream at last comes into play; it's called aStream. Remember that we stripped all zeros above, so we'll handle it specially. First print the first character, which is the integer part:

digitStream atEnd
ifTrue: [aStream nextPut: \$0].
ifFalse: [aStream nextPut: digitStream next].

Then the decimal point and everything else:

aStream nextPut: \$..
digitStream atEnd
ifTrue: [aStream nextPut: \$0]
ifFalse: [aStream nextPutAll: digitStream].

and finally the exponent letter:

aStream
nextPut: \$e;
print: exponent

That's it. It is not incredibly complex, but also not too optimized and full of off-by-one traps. So, good floating point to decimal conversion is not so easy!

Why would you not go from the binary format of the floating point number? part of the float is exp with a bias, part of the float is the mantissa or fraction and part is the sign. The float guys have busted a gut getting all the rounding right. if you just interpret their resulting bits and convert it to ascii, you'd be pretty good. (yeah, there are a couple of different floating formats but less than a handful are important. as a float vlsi designer, with the speed and memory of today's machines, why should we do float at all; just do ascii like we did in grammar school and manipulate the symbols.