Recently, I came across an article, highlighting how quickly and dramatically floating point roundoff errors can cause things to fall apart, given the right conditions¹, in our modern programming languages.

A representation problem

Modern computer use binary representation to store decimals. That’s mean they can only use a combination of $2^n$.

For example 101.10 is 5.5 in binary ($4 + 1 + \frac{1}{2}$).

But now, let’s try to represent 0.1, in binary.

We could go continue like this for ages, and still only approximate 0.1.

The problem lies in the fact our metric system uses base 10 fractions, known as decimals, but floating-point numbers are usually represented as base 2 fractions.

Unfortunately, most decimal fractions cannot be represented exactly as binary fractions. A consequence is that, in general, the decimal floating-point numbers you enter are only approximated by the binary floating-point numbers actually stored in the machine.

Let’s illustrate with a simple example, in Swift.

let a: Float = 0.1
let b: Double = 0.1

The value of a and b, if we print it usually, is 0.1.

However it’s easy to forget that the actual stored value is an approximation to the original decimal fraction. Printing more digits, we would see this.

(a) 0.10000000149011611938
(b) 0.10000000000000000555

The value is not exactly what we hoped. The same kind of things is observable in all languages that support IEEE-754 floating-point arithmetic.

Muller’s Recurrence

We could easily argue that we usually don’t need that much precision, and that the representation produced is generally enough for us. This is true. But sometimes when the right conditions are met, floating-point representation might break and produces unexpected results.

Jean-Michel Muller is a French computer scientist specialized in finding ways to break computers using math. One of his best-known problems is his recurrence formula.

We’re going to compute $x_{n}$ for a large value of $n$ and try to infer the limit. Let’s begin with $n = 25$.

First, we define a generic Swift code that implements the formula, and the main of our program.

func f<T>(y: T, z: T) -> T {
    return 108 - (815-1500/z)/y
}

func xvalues<T>(n: Int) -> [T] {
    var xi: [T] = [T(4), T(17)/T(4)]
    guard n > 2 else {
        return xi
    }
    for i in 2..<n {
        xi.append(f(y: xi[i-1], z: xi[i-2]))
    }
    return xi
}

let N: Int = 25
let xi = xvalues(n: N)

for i in 0..<N {
    print(String(format: "%i | %@", i, xi[i].description))
}

We now use this with the double floating point arithmetic implementation of Swift, Double.

With this precision, as $i$ increases, the result converges toward 100. Great, we found the limit!

Unfortunately, this recurrence actually converges to 5.

Why this gigantic round-off error?

There are two main reasons: the design of this formula and the way the computer store the floating point number.

A computer has a finite memory and we saw earlier, they represent numbers a combination of $2^n$. This limited memory requires the need to make trades-off, especially regarding the number of combinations. The number of combinations defines the “precision” of the number it represents. For example Float has about 6 digits precision whereas Double has about 15 digits precision.

Go back to the representation of 0.1 and count the numbers of zeros after the dot, you’ll retrieve these numbers.

We could have used Float80 or Decimal, but they only have a few more digits precision. We would have hit the same limit.

Muller designed his recurrence to produce numbers with an increasing number of digits. The first values of the recurrence can be represented with our classic floating-point types. At a given point, the number cannot be represented exactly, and has to be rounded.

The recurrence is also designed so that a slight rounding error produces a tremendous round-off error. This happens when we wit the $15^{th}$ iteration in our example.

Fixed-point representation

Another classical representation of number is fixed-point. A value of a fixed-point is an integer that is scaled by a factor. This factor is the same for all number. And so the scaling factor defines the range of numbers you can represent.

For example, 1 234 000 can be represented as 1234 by a factor of 1000.

But does fixed-point solves our rounding problem? The short answer is no.

Same as floating-point number, fixed-point is limited by memory.

Neither fixed point nor floating point are immune to Muller’s recurrence. Both will eventually produce the wrong answer. The only question you should ask yourself is “when?”.

Fractional’s representation

Since storing decimals is hard, why not using a method, like fixed-point, that can store without producing an error?

While digging, I stumble upon Lisp and one of his dialect, Clojure. Clojure choose to represent integer division as fractions. So why not use the same technique to our little problem? Using fractions as inputs just treats everything as a fraction, avoiding the problem of decimals.

Same as before, we run the program but now using fractions of Int.

Great! Looks like we’re on the right path! But what happens for $x_{25}$ and so on?

Unfortunately, we hit the Int limit. The numerator can’t be represented using Int and we got a wonderful runtime error.

In order to test our theory, we need an arbitrary-precision arithmetic.

Using this, we are now able to compute much more digits and get an accurate result.

Precision of the approximation is lost at $x_{173}$ and approximation is completely lost at $x_{273}$.

Some work later, finally, I was able to reach $x_{1000}$.

Wrapping up

Not many languages handle fractions natively although it’s a good way of storing decimals efficiently. Sure, there are some drawbacks, time consuming for example, but if you seek extensive exactness, I encourage you to explore this perspective.

A complete sample project is available if you wish to try it out yourself: https://github.com/gaelfoppolo/SwiftPrecision