Study into exact real arithmetic Pt. 2

Reading the previous post. first is good idea because here we're going to go into some of the details.

I did modifications into the original IC-Reals-Haskell in the exact-computation page in order to understand it better. These modifications are available for you to download (ICReals.hs). The same license is retained as in the original copy.

"Exact Real Arithmetic using Möbius Transformations" by Peter John Potts explains what is going on in the library. I used both the paper and the code to understand the subject to the point that I can understand the structure of the library. Though I still do not entirely understand all the details so it's good to know where to read further on.

I'm leaving out the elementary functions from the explanation. For now it's sufficient to understand that the library can compute exactly something such as sin(1.4*π) + log(5).

Expression trees

The expression trees are in the heart of the system.

data Expression = ExpV Vector
                | ExpM Matrix Expression
                | ExpT Tensor Integer Expression Expression
                deriving (Eq, Show)

Note the only terminating element is the ExpV Vector. It is a rational number expressed with nominator/denominator pair and it seems necessary that the terminator identifies a rational number. Otherwise you end up with a different number domain with its own problems.

The matrices and tensors are like you'd expect based on the previous post, though there's a counter in the tensor. This counter is controlling the absorption strategy which is necessary to not get stuck during computation.

Square bracket application

Eventually we have to apply the LFT to something. The result tends to be some other kind of an LFT.

app  :: Matrix -> Expression -> Expression
tapp :: Tensor -> Integer -> Expression -> Expression -> Expression

I separated the application for matrices and tensors because they're quite different.

If we applied tensor into an expression that is itself a tensor then we'd get a higher-order element. To avoid them there's an absorption rule that converts the tensor into an exact floating point. We could also enforce this with types, using the structure:

data Expression = ExpV Vector
                | ExpM Matrix Expression
                deriving (Eq, Show)

data ExpressionT = ExpT Tensor Integer Expression Expression
                deriving (Eq, Show)

Then the user would themselves do utoe (tem t) to connect a tensor into an expression tree.

This kind of modification is probably a good idea, but may require some further rethinking of the library.

For emission I provided a modification of the bracket-application that checks that the refinement property holds for the resulting expression.

rapp :: Matrix -> Expression -> Maybe Expression

This avoids little bit of computation of inverses during the emission during emission. Basically it appears in these things:

fmap (Szer . dem) (rapp iszer e)

The fmap glues more stuff into the Maybe if there's a value. Effectively we calculate the Szer (dem (app iszer e)) if the refinement property holds, otherwise the expression evaluates to nothing. Success would mean that the expression emits a symbol. In this case it'd emit a Z-sign which means that the number drops into the [-1 1] -interval.

Refinement property

I'm still not sure what the refinement property means for a tensor. But we basically sum the signs in an expression and get the sign.

So say we have some constant pair a, b such as:

fa
──
fb

Then we would compute sign(a/b) = sign(sign(a) + sign(b)), examples:

sign(+a/+b) = sign(1 + 1) = +
sign(-a/+b) = sign(-1 + 1) = 0
sign(-a/-b) = sign(-1 + -1) = -

Then we could express LFT satisfying refinement property with:

xy(+) x(+) y(+) (+)
xy(-) x(-) y(-) (-)
x(+) (+)
x(-) (-)

In other words, if the columns in matrix are same sign and nonzero, then the refinement holds for the LFT represented by the matrix.

Refinement property would therefore hold for multiplication, division, addition, reciprocal LFT, but not for subtraction or negation LFT.

The goal with this thing is probably to ensure the interval is narrowing in the digit emission.

Normalization

The normalizer emits a symbol if it can. Otherwise it does an absorption step and tries again.

The sem emits a sign.

sem :: Expression -> Sefp
sem (ExpV v) = choices [
    fmap (Spos . dem) (rapp ispos (ExpV v)),
    fmap (Sneg . dem) (rapp isneg (ExpV v))] (error "sem (ExpV v)")
sem e = choices [
    fmap (Szer . dem) (rapp iszer e),
    fmap (Sinf . dem) (rapp isinf e),
    fmap (Spos . dem) (rapp ispos e),
    fmap (Sneg . dem) (rapp isneg e)] (sem (em e))

The dem emits a digit.

dem :: Expression -> Uefp
dem (ExpV v) = Term v
dem e = choices [
    fmap (Dneg . dem) (rapp idneg e),
    fmap (Dpos . dem) (rapp idpos e),
    fmap (Dzer . dem) (rapp idzer e)] (dem (em e))

The em checks the absorption step is done according to the strategy. For matrices we just basically apply the matrix into the inner expression. In tensors we have to operate through a strategy.

em (ExpM m e)       = app m e
em (ExpT t i e1 e2) = tapp t i
    (ab e1 ((strategyo t i) == 1))
    (ab e2 ((strategyo t i) == 2))

The tensor absorption strategy may decide to delay the absorption of either side in an expression. Actually it'd seem that we never absorp the both sides at once, so this could be simplified further.

The second element in the absorption selection casts the 2-LFT into 1-LFT as mentioned before. Otherwise we'd end up with 3-LFTs that haven't been programmed into the system.

ab e b
    | not b     = ExpM identity e
    | tis e     = utoe (dem e)
    | otherwise = e

tis identifies an expression that is a tensor.

tis                  :: Expression -> Bool

The choices -expression is used during emission to select a digit if it can be emitted.

choices :: [Maybe a] -> a -> a
choices [] e = e
choices (Just a : rest) e = a
choices (Nothing : rest) e = choices rest e

Finally we only have to examine the tensor absorption strategy to finish the examination of the library for now.

LFT absorption strategy

It appears that the left/right absorption in a 2-LFT is chosen based on whether the refinement property holds.

strategyo t i
    | trefine t             = strategyr t i
    | otherwise             = mod i 2 + 1

If the simple strategy is used, then the absorption is done by alternating the side in absorption.

For refined elements the program decides the absorption by checking whether the columns in 2-LFT are disjoint.

strategyr t i
    | mdisjointm t1 t2      = 2
    | otherwise             = 1
    where
    t1 = fst (transpose t)
    t2 = snd (transpose t)

Pott's explains in the paper that this is avoiding some issues that would otherwise cause the computation to hang up.

The disjointness probably has something to do with the disjointness of matrices.

mdisjointm m n    = mlessm m n || mlessm n m
mlessm m (v,w)    = mlessv m v && mlessv m w
mlessv (v,w) x    = vlessv v x && vlessv w x
vlessv v w        = determinant (v,w) < 0

I haven't done much with determinants in matrices. I'd guess it may be obvious if I do a similar analysis as before.

    az + b          ew + f
z ↦ ──────      w ↦ ──────
    cz + d          gw + h

So we break the right-side matrix into ew/gw and f/h. Then we do the same for left-side matrix and check the determinants.

ag - ce < 0
ah - cf < 0
bg - de < 0
bh - df < 0

But.. This gets done to the transposed matrices, meaning we compare nominators, denominators.

Didn't figure it out.

Conclusion

If you wrap the whole thing into a continuation, you can create a calculator where user can input as many digits as he knows, and he gets a result where there's as many digits as can be produced.

The redundant number representations have potentially many other uses.

In the ICreals library they ensure that the computation can emit a digit that denotes where the result falls, without needing to know whether the result goes above or below it. The whole result produced by the library is sort of like emitting digits, each digit pointing out an interval of where the number falls.

Redundant binary representations can break carry chains by memoizing carry flags in the result when doing binary computations.

If we used RNR representations on price tags, we could prevent common form of carry-number-fraud common in modern commerce. Merchants could be forced to normalize numbers such that price tags always have maximal number of zero digits. eg. 499.9¤ would be denoted as 500.1̅¤, and you'd know the thing costs 500 credits minus 0.1. By simple examination we find out this rule gives surprisingly little leeway to giving numbers that are larger than they seem.

499.1¤ --> 500.9̅¤
498.0¤ --> 509̅.0¤
490.0¤ --> 490.0¤
490.9¤ --> 490.9¤

This price presentation would still be redundant. 15 and 25̅ are same value. The redundancy could be used in a price-slashing-fraud. It could be prevented by requiring that the representation with smallest absolute value is used if there are multiple maximal-zero presentations.

Finally we probably want to compute the sin(1.4*π) + log(5).

eshow
    (ExpT tadd 0
        (esin (ExpT tmul 0 (ExpV (14,10)) epi))
        (elog (ExpV (5,1))))
    10

Oh... that hangs up. It looks like I found a bug in the original source code. Lets see what it could be. I had no idea at this point.

Bugfixing

Examining at the eshow, it could either be the printing or normalization procedure.

*ICreals> sem (ExpT tadd 0 (esin (ExpT tmul 0 (ExpV (14,10)) epi)) (elog (ExpV (5,1))))
^CInterrupted.

It looks like the normalization expression hangs up. We already know that it's not getting to emission, so we know that the problem is in the em -function.

*ICreals> let foo = (ExpT tadd 0 (esin (ExpT tmul 0 (ExpV (14,10)) epi)) (elog (ExpV (5,1))))
*ICreals> em foo
^CInterrupted.

I was expecting em would print something, but no. This was interesting. On closer examination of the expression, it'd seem that utoe (dem e) is evaluated.

*ICreals> foo
ExpT (((0,0),(1,0)),((1,0),(0,1))) 0 ...
*ICreals> let ExpT t i a b = foo
*ICreals> a
ExpT (((0,1),(-1,0)),((-1,0),(0,1))) 2 ...
*ICreals> b
ExpT (((1,0),(1,1)),((-1,1),(-1,0))) 0 ...

When I evaluated the culprit, it hang up.

*ICreals> utoe (dem a)
^CInterrupted.

Oh right, this expression is not emitting the sign-symbol and it should since the sin returns a negative number.

*ICreals> stoe (sem a)
ExpM ((0,1),(-1,0)) ...

When I fix the issue, it seems to work, but I should check that the issue is out.

*ICreals> eshow (ExpT tadd 0 (esin (ExpT tmul 0 (ExpV (14,10)) epi)) (elog (ExpV (5,1)))) 100
"0.65838139613894680248431999984e0"
*ICreals> eshow (esin (ExpT tmul 0 (ExpV (14,10)) epi)) 100
"-0.95105651629515357211643933338e0"
*ICreals> eshow ((elog (ExpV (5,1)))) 100
"0.160943791243410037460075933323e1"

I guess it's time to do bit of counting by hand.

+1.60943791243410037460075933323e0
-0.95105651629515357211643933338e0
────────────────────────────────── +
 0.65838139613894680248431999985e0
────────────────────────────────── compare
"0.65838139613894680248431999984e0"

Ok. It seems to work after the fix. Anyway that last number is the approximate answer to the earlier question.

That bug could have had a driver's license as it had well-reached a legal drinking age. It was 21 years old and appeared in Pott's own Miranda code.

Some applications for exact real arithmetic

This library can be used for simple back-of-the-envelope calculations. It gives you exact approximate answers, aka. Most significant binary digits and a computation to get more of them.

Exactness in real arithmetic is surprisingly pleasing and the utilities provided so far provide enough functionality to apply it into clifford algebra.

I also remain wondering whether there were applications in this for symbolic computer algebra and constraint solving.

As a final example we can compute how far the horizon is on a sphere model that has the radius of the Earth.

sqrt((height + radius)² - radius²)

Here's the Haskell program for doing it with ICreals:

module EarthHorizon where

import ICreals

-- Note that no way of telling that this number is approximation
-- means we won't get entirely correct answers.
earth_radius = ExpV (6371000, 1)

horizon_distance view_height = esqrt $ sub
    (square (add view_height earth_radius))
    (square earth_radius)

add = ExpT tadd 0
sub = ExpT tsub 0
mul = ExpT tmul 0
square a = mul a a

This takes several minutes to compute.

*EarthHorizon> eshow (horizon_distance (ExpV (1,1))) 30
"0.356959e4"
(467.28 secs, 60,477,217,552 bytes)

The distance 3569.59m agrees with the result I got from elsewhere. I doubt the accuracy of 10 millimeters is warranted though.

Effort taken (60 gigabytes) is insane. :D But at least finally I figured out how to make use of all this computing capacity! In an useful sort of way, sort of.

My plans for this is to use it in a 3D-modeling environment. I guess I'd need better performing code&algorithms though.