# Study into exact real arithmetic

London's Imperial College has a page for exact computation that holds plenty of nice things. We are going to look into the IC-Reals-Haskell.

Goal of this post is to obtain some understanding over how it's possible to do exact real arithmetic with a computer.

## Continued fractions

Continued fractions are fractional number representations that can extend to infinity.

Here's a generalized continued fraction that exactly represents π:

`pi = 2 + 2 / f 1`

`f n = 1 + n*(n+1) / f (n+1)`

You can open up the recursion and get the following sequence:

`z ↦ 2 + z`

`z ↦ 2 / (1 + z)`

`z ↦ 1*2 / (1 + z)`

`z ↦ 2*3 / (1 + z)`

`z ↦ 3*4 / (1 + z)`

`n ↦ z ↦ (n-1)*n / (1 + z)`

`...`

Now the program that exactly represents π
can be written as `π°(π¹(π²(π³(...))))`

.

`pi_lcf 0 z = 2 + z`

`pi_lcf 1 z = 2 / (1 + z)`

`pi_lcf n z = (n-1)*n / (1 + z)`

`pi_sequence :: (Enum a, Eq a, Fractional a) => [a -> a]`

`pi_sequence = [pi_lcf n | n <- [0..]]`

This decomposition allows writing a Haskell program that computes convergents for π.

`compute [] z = z`

`compute seq z = (head seq) (compute (tail seq) z)`

`convergent seq n = compute (take n seq) 0`

We get successively more accurate approximations.

`convergent pi_sequence 1 ≡ 2.0`

`convergent pi_sequence 2 ≡ 4.0`

`convergent pi_sequence 3 ≡ 2.6666666666666665`

`convergent pi_sequence 4 ≡ 3.555555555555556`

`convergent pi_sequence 5 ≡ 2.8444444444444446`

`convergent pi_sequence 6 ≡ 3.413333333333333`

`convergent pi_sequence 7 ≡ 2.9257142857142857`

`convergent pi_sequence 8 ≡ 3.3436734693877552`

`convergent pi_sequence 9 ≡ 2.972154195011338`

`convergent pi_sequence 10 ≡ 3.3023935500125976`

`convergent pi_sequence 11 ≡ 3.0021759545569067`

`convergent pi_sequence 12 ≡ 3.275101041334808`

`convergent pi_sequence 13 ≡ 3.023170192001361`

`convergent pi_sequence 14 ≡ 3.2557217452322345`

`convergent pi_sequence 15 ≡ 3.0386736288834193`

`convergent pi_sequence 16 ≡ 3.2412518708089806`

`convergent pi_sequence 17 ≡ 3.050589996055511`

`convergent pi_sequence 18 ≡ 3.2300364664117174`

`convergent pi_sequence 19 ≡ 3.0600345471268904`

`convergent pi_sequence 20 ≡ 3.221088996975674`

`convergent pi_sequence 200 ≡ 3.149456428081303`

`convergent pi_sequence 201 ≡ 3.133787490628162`

`convergent pi_sequence 2000 ≡ 3.1423781499034096`

`convergent pi_sequence 2001 ≡ 3.1408077460303945`

`convergent pi_sequence 20000 ≡ 3.1416711943878566`

`convergent pi_sequence 20001 ≡ 3.141514118681922`

`convergent pi_sequence 200000 ≡ 3.1416005075812445`

`convergent pi_sequence 200001 ≡ 3.141584799657246`

`convergent pi_sequence 2000000 ≡ 3.1415934389880547`

`convergent pi_sequence 2000001 ≡ 3.141591868192121`

`actual pi ≡ 3.141592653589793238462643383279..`

The difference between convergents give clues about how good approximation you have at any moment.

This method to compute π is not good for that purpose, but hopefully it revealed something about computable reals. The representations for computable reals are programs that give you increasingly better answers.

Here's the output from IC-Reals-Haskell -library.

`eshow epi 1 ≡ "unbounded"`

`eshow epi 10 ≡ "0.31e1"`

`eshow epi 100 ≡ "0.314159265358979323846264338328e1"`

`eshow epi 1000 ≡ "0.3141592653589793238462643383279502884197169399..e1"`

## Linear fractional transformations

Linear fractional transformation is a function of the following form, that has an inverse:

`az + b`

`z ↦ ────── {ad − bc ≠ 0}`

`cz + d`

LFTs are represented as matrices and they behave a bit like those, but you should remember that they are functions that return fractions.

Matrix multiplication can be used to compose LFTs:

`(w ↦ (aw + b) / (cw + d)) . (z ↦ (ez + f) / (gz + h))`

`ez + f`

`────── a + b`

`gz + h`

`z ↦ ────────────`

`ez + f`

`────── c + d`

`gz + h`

`multiply fraction by gz + h and reorder terms`

`(ae + bg)z + af + bh`

`z ↦ ────────────────────`

`(ce + dg)z + cf + dh`

If we choose a base interval, then the LFT becomes an interval representation.
The intervals commonly used for this purpose are `[-1 1]`

and `[0 ∞]`

.
The IC-Reals-Haskell appears to use the zero-to-infinity interval.

If we plug the interval into the LFT, we get the interval it represents:

`z+0`

`N = z ↦ ─── [0 ∞]`

`z+2`

`(0+0) / (0+2) = 0`

`(∞+0) / (∞+2) = (1/0+0) / (1/0+2) = 1/1 = 1`

`represents: [0 1]`

`3z+1`

`0 = z ↦ ──── [0 ∞]`

`1z+3`

`represents: [1/3 3/1]`

`2z+1`

`P = z ↦ ──── [0 ∞]`

`+1`

`represents: [1 ∞]`

The above LFTs form digits that can be used for representing any number on a real number line as a succession of intervals that eventually shrink to zero.

`0N lfc={z ↦ (2z + 1) / (2z + 3)} interval=[1/3 1/1]`

`00 lfc={z ↦ (5z + 3) / (3z + 5)} interval=[3/5 5/3]`

`0P lfc={z ↦ (3z + 2) / (1z + 2)} interval=[1/1 1/3]`

The full notation comes with signs to cover the whole number line:

`z`

`POS = z ↦ ─ [0 ∞]`

`1`

`z+1`

`INF = z ↦ ──── [1 -1]`

`-z+1`

`1`

`NEG = z ↦ ─ [∞ 0]`

`z`

`z-1`

`ZER = z ↦ ─── [-1 1]`

`z+1`

This notation for real numbers has redundancy which means there are many representations for the same real as intervals overlap. The redundancy is there to allow intervals represented by the digits to sway.

Additionally the library is using a larger, "tensor" LFT to represent binary arithmetic.

`xya + xb + yc + d`

`(x, y) ↦ ─────────────────`

`xye + xf + yg + h`

It can represent multiplication, division, addition, subtraction and their combinations to two numbers this way.

## Decimal output

The `eshow`

seem to be the best point to start examining how the library works.
Haskell representation is extraordinarily clean and would be even better if it
had been documented a little bit more.

So lets resume to the pi -example. What is this thing? They say that you get about third of digits for the iteration number that appears in the expression.

`eshow epi 20 ≡ "0.31416e1"`

Why does the number 20 translate to about 5 significant digits? That's because the iteration number is number of digits in that exact binary float representation just mentioned.

The implementation of eshow is:

`eshow e i = mshow (stom (sem e) i)`

The `sem epi`

-function normalizes the digit
and produces a signed, exact floating point.

`Spos (Dpos (Dzer (Dzer (Dzer (Dzer (Dpos (Dzer (Dzer`

`(Dzer (Dpos (Dpos (Dzer (Dzer (Dzer (Dzer (Dzer`

`(Dpos (Dzer (Dneg (Dneg ...`

`+1000 0100 0110 0000 10NN..`

The floating point is collapsed into a matrix, with matrix multiplication rules.
eg. `stom {+1}`

is equal to `((2,0),(1,1))`

.

`stom (sem epi) 0 ≡ ((1,0),(0,1))`

`stom (sem epi) 1 ≡ ((2,0),(1,1))`

`stom (sem epi) 2 ≡ ((7,1),(5,3))`

`stom (sem epi) 3 ≡ ((13,3),(11,5))`

`stom (sem epi) 4 ≡ ((25,7),(23,9))`

`stom (sem epi) 5 ≡ ((49,15),(47,17))`

`stom (sem epi) 6 ≡ ((49,15),(48,16))`

`stom (sem epi) 7 ≡ ((195,61),(193,63))`

`stom (sem epi) 8 ≡ ((389,123),(387,125))`

`stom (sem epi) 9 ≡ ((777,247),(775,249))`

`stom (sem epi) 10 ≡ ((777,247),(776,248))`

`stom (sem epi) 11 ≡ ((1554,494),(1553,495))`

`stom (sem epi) 12 ≡ ((6215,1977),(6213,1979))`

`stom (sem epi) 13 ≡ ((12429,3955),(12427,3957))`

`stom (sem epi) 14 ≡ ((24857,7911),(24855,7913))`

`stom (sem epi) 15 ≡ ((49713,15823),(49711,15825))`

`stom (sem epi) 20 ≡ ((795395,253181),(795394,253182))`

The `mshow`

converts the matrix into a decimal string.
Matrices with zero determinant appear to represent exact rationals.

`mshow ((3,4), (3,4)) ≡ "3/4"`

`3z + 3 3(z+1) 3`

`────── = ────── = ─`

`4z + 4 4(z+1) 4`

The function `scientific`

contains the core logic
for extracting a decimal mantissa from a matrix.

The next explanation is going to be lacking because I don't understand the refinement property that the program is using.

If the interval of the matrix is unbounded, the range in denominator crosses zero. The function checks the condition and represents it with an empty sequence. Otherwise it starts outputting an exponent.

If the number in the matrix is too large, it increments the exponent and divides denominator in matrix by 10. Otherwise the mantissa is built.

When a number is selected to the mantissa, the LFT-matrix gets multiplied by 10 and the digit is subtracted away from it. Those are done by multiplying it with the following LFT:

`10z - n`

`z ↦ ───────`

`0z + 1`

The scientific number produced seem to itself be redundant as well. To confirm this you can give it a LFT such as:

`4z +6`

`z ↦ ───────`

`10z +10`

The program outputs `0.5e0`

. Basically you can mark it as 0.5e0±1:
The number is something between 0.4 and 0.6.

Despite operating with intervals, the library doesn't let you give in an interval. Every-non-infinite expression is terminated with a rational number.

`eshow (elog (ExpV (4, 10))) 20 ≡ "-0.91629e0"`

`eshow (elog (ExpV (5, 10))) 20 ≡ "-0.69315e0"`

`eshow (elog (ExpV (6, 10))) 20 ≡ "-0.51083e0"`

It's something to look for. Sometimes it could really help if you knew how much accuracy you can get with your input.

`eshow (elog (ExpV (501, 1000))) 20 ≡ "-0.69115e0"`

`eshow (elog (ExpV (499, 1000))) 20 ≡ "-0.69515e0"`

Here you see you need three digits in 0.5 before you get a result from logarithm function that is accurate to two digits.

So far I've only understood how the number format ticks. There's still the evaluator and the actual implementations of trigonometric functions, pi, and such left for me to explore. I really have to find the correct papers to examine because some of the clauses are just magic numbers all over it.

Well I'll get to that at some point. Now it's time to resume on the other projects.