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.

Similar posts