### Wed, 21 Mar 2007

#### Xtreme Numerical Accuracy.

I'm working on a digital filter design program in Ocaml which was suffering from some numerical issues with Ocaml's native 64 bit floats. The problem was that the algorithm operates on both large floating point numbers and small floating point numbers. These numbers eventually end up in a matrix, and I then use Gaussian elimination to solve a set of simultaneous equations.

Anyone who has done any numerical computation will know that adding large floating point numbers to small floating numbers is a recipe for numerical inaccuracy. For me, the numerical issues were screwing things up badly.

When faced with a problem like this there are two possible solutions:

- Do all the computations symbolically, and only substitute numbers at the very last stage and then being careful to process the numerical parts in a way to minimize rounding and truncation error.
- Replace the floating point operations with operations on a number type which can provide higher (and maybe even arbitrary) precision.

The first option, doing all the computations symbolically was not practical due to the complexity of the computation. That left only the second option.

Looking around for what was available for Ocaml, I found the contfrac project on Sourceforge. As all the math geeks (hi Mark) have probably guessed by now, contfrac expresses numbers in terms of a really cool mathematical concept called continued fractions.

The idea is that any number can be represented by a (potentially infinite) list
of integers **[ a0 ; a1, a2, a3, ...]**.
Given the list of integers, the number itself can be calculated using:

All
rational numbers
have a finite length continued fraction expansion.
For example, the rational number 75/99 is expressed as **[ 0 ; 1, 3, 8 ]**.

Not surprisingly, all the
irrational numbers
have infinite length continued fraction expansions.
The surprising thing (for me at least) is that many of the irrational numbers
have CF expansions that are surprisingly regular.
The square root of two is expressed as **[ 1 ; 2, 2, 2, ...]** with an
infinitely repeating list of 2s.
The natural logarithm **e** is expressed as **[ 2 ; 1, 2, 1, 1, 4, 1, 1,
6, ...]** which again has a regular pattern, as does the golden ratio,
**[ 1 ; 1, 1, 1, ...]**.
While all the previous CF expansions have a degree of regularity, the expansion
of **pi**, is **[ 3 ; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2,
2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13,...]**, which looks completely random.

With numbers expressed as continued fractions, the Ocaml contfrac module then implements addition, subtraction, multiplication and division. Once the four arithmetic operations are defined, contfrac then implements a number of trigonometric and transcendental functions in terms of the same continued fractions.

Unfortunately, the module doesn't implement everything I need so I'm going to have to hack on some extra functionality. The actual Ocaml implementation uses Ocaml's lazy lists which is an aspect of Ocaml I hadn't played with yet. Time for some fiddling with lazy lists.

Posted at: 20:49 | Category: CodeHacking/Ocaml | Permalink