Tag: Exact arithmetic

Getting started with exact arithmetic and F#

In this blog post, I claimed that some exact arithmetic beyond rational numbers can be implemented on a computer. Today I want to show you how that might be done by showing you the beginning of my implementation. I chose F# for the task, since I have been waiting for an opportunity to check it out anyway. So this post is a more practical (first) follow up on the more theoretic one linked above with some of my F# developing experiences on the side.

F# turned out to be mostly pleasant to use, the only annoying thing that happened to me along the way was some weirdness of F# or of the otherwise very helpful IDE Rider: F# seems to need a compilation order of the source code files and I only found out by acts of desperation that this order is supposed to be controlled by drag & drop:

The code I want to (partially) explain is available on github:

https://github.com/felixwellen/ExactArithmetic

I will link to the current commit, when I discuss specifc sections below.

Prerequesite: Rational numbers and Polynomials

As explained in the ‘theory post’, polynomials will be the basic ingredient to cook more exact numbers from the rationals. The rationals themselves can be built from ‘BigInteger’s (source). The basic arithmetic operations follow the rules commonly tought in schools (here is addition):

static member (+) (l: Rational, r: Rational) =
    Rational(l.up * r.down + r.up * l.down,
             l.down * r.down)

‘up’ and ‘down’ are ‘BigInteger’s representing the nominator and denominator of the rational number. ‘-‘, ‘*’ and ‘/’ are defined in the same style and extended to polynomials with rational coefficients (source).

There are two things important for this post, that polynomials have and rationals do not have: Degrees and remainders. The degree of a polynomial is just the number of its coefficients minus one, unless it is constant zero. The zero-polynomial has degree -1 in my code, but that specific value is not too important – it just needs to be smaller than all the other degrees.

Remainders are a bit more work to calculate. For two polynomials P and Q where Q is not zero, there is always a unique polynomial R that has a smaller degree such that:

P = Q * D + R

For some polynomial D (the algorithm is here).

Numberfields and examples

The ingredients are put together in the type ‘NumberField’ which is the name used in algebra, so it is precisely what is described here. Yet it is far from obvious that this is the ‘same’ things as in my example code.

One source of confusion of this approach to exact arithmetic is that we do not know which solution of a polynomial equation we are using. In the example with the square root, the solutions only differ in the sign, but things can get more complicated. This ambiguity is also the reason that you will not find a function in my code, that approximates the elements of a numberfield by a decimal number. In order to do that, we would have to choose a particular solution first.

Now, in the form of unit tests unit tests, we can look at a very basic example of a number field: The one from the theory-post containing a solution of the equation X²=2:

let TwoAsPolynomial = Polynomial([|Rational(2,1)|])
let ModulusForSquareRootOfTwo = 
     Polynomial.Power(Polynomial.X,2) - TwoAsPolynomial
let E = NumberField(ModulusForSquareRootOfTwo)   
let TwoAsNumberFieldElement = NumberFieldElement(E, TwoAsPolynomial)

[<Fact>]
let ``the abstract solution is a solution of the given equation``() =
    let e = E.Solution in  (* e is a solution of the equation 'X^2-2=0' *)
    Assert.Equal(E.Zero, e * e - TwoAsNumberFieldElement)

There are applications of these numbers which have no obvious relation to square roots. For example, there are numberfields containing roots of unity, which would allow us to calculate with rotations in the plane by rational fraction of a full rotation. This might be the topic of a follow up post…

There are more exact numbers than you might know

There are lots of annoying problems a programmer can run into, when working with the inexact floating point numbers supported by common cpus. So the idea of using rational numbers

\frac{n}{d}

with arbitrary size integers n,d for floating point computation can be quite appealing. However, if you really try this in practice, you are likely to end up with unacceptable performance (for an example, see the post The Great Rational Explosion on this blog).
Another problem with the rationals is, that they don’t support lots of mathematical operations. For example, any prime will not have a rational square root and numbers like \pi and e are not rational.
This post is about the problem of the missing square roots. I will sketch how the rationals can be extended by some specific real numbers. The focus is on what can be done and not how it can be done.

The point of the following example is that we can extend the rational numbers by some new numbers such that we have the square root of 2 and are still able to perform the basic arithmetic operations, i.e. addition, subtraction, multiplication and division by non-zero numbers. And of course, we still want to have the absolute exactness of the rationals.
To achieve this goal, our new numbers can be represented as tuples (a,b) of rational numbers a,b. The idea is, that (a,b) stands for the real number

a+b\sqrt{2}

And now, a small miracle happens: If you have two of those new numbers then all of the basic arithmetic operations are again given as one of our new numbers! You can believe me or just check this for yourself now. So in total, we have extended the realm of exact calculation to real numbers which are of the form

a+b\sqrt{2}.

This is admittedly a bit lame, since we are still far away from calculating general square roots. But the miracle mentioned above is actually not so small and happens in way more generality: The real number \sqrt{2} that we added to the rationals is the solution of the polynomial equation

X^2=2 or equivalently X^2-2=0

and we can do the same thing with any irrational real number that is the solution of some polynomial equation

a_n X^n + a_{n-1} X^{n-1} + \dots + a_1 X^1 + a_0 = 0 or P(X)=0

(where a_n,\dots,a_0 are rational numbers).
For this general construction, the polynomial above has to be chosen such that n is as small as possible. Then the new numbers can be represented as tuples of length n. For the basic arithmetic operations, we think of a tuples (b_n,\dots,b_0) as the polynomial

b_n X^n+\dots b_1 X^1 + b_0

and define addition to be addition of those polynomials. To multiply two tuples, take the product of the polynomials and its remainder on division by P. You can check out how this recipie does the right thing in the example with P=X^2-2. To avoid confusion, let us put the polynomial defining in which extension of the rationals we are into the number. So know, a number is a tuple

(P,b_n,\dots,b_0).

In particular, there are extensions that contain a square root of an arbitrary rational. But if we look at, say, \sqrt{2} and \sqrt{3} we will get two extensions. So we have \sqrt{2} and \sqrt{3} but we can’t do anything that involves both of them. Or in other words, if we have two numbers, (P,b_n,\dots,b_0) and (Q,c_m,\dots,c_0) we don’t know how to perform basic arithmetic unless n=m and P=Q.

Here is an easy way to fix this: If we take a step back, we can realize that everything we did above doesn’t need specific features of the rational numbers. The basic arithmetic operations are enough. So we can apply the construction of new numbers including some root to an extension of the rationals. For the problem mentioned above, that would give us tuples of the form

(X^2-3, (X^2-2, a, b), (X^2-2, c, d))

– now with four rational numbers a,b,c,d.
The corresponding general pattern for an extension (of some extension …) is:

(P, a_0,\dots, a_n)

Where the numbers a_i are from some fixed extension and P is a polynomial with coefficients from this fixed extension. We also have to make sure, that P has a property called irreducible. It means that P has at least degree 1 and there is no pair of polynomials Q,R with lesser degree (and coefficients in the extension) such that their product P\cdot R is equal to P. If this is not the case, weird things will happen.
Note that if we keep extending our numbers on demand, we can now get square roots of everything, in fact, this even works for square roots of negative numbers. But there are some unmentioned problems, if we want to mix numbers from different towers of extensions. On the other hand, there are also lots of unmentioned possibilites: These construction are pretty general and they can also be applied to finite fields or some extension of the rationals which contains a trancendent element like \pi. But I will stop here since my goal of showing some exact representations of numbers that surprised myself once, is reached.

If you want to learn more about the mathematical background, this pdf might help. There is also a way to join two extensions to one of the form we discussed here, the key is this nice theorem: Primitive Element Theorem. There are lots of nice mathematical facts in this area. My former Algebra working group in Karlsruhe, in particular Fabian Ruoff, kindly reminded me of some of them over lunch, when I was discussing the topics of this post.
Then, there are of course implementations. Here is a class in the cgal library for square roots.