Constructive reals

As already mentioned the main idea behind the constructive reals is to represent a real number \(x\) not as a (possibily infinite) sequence of digits but as a rational function \(\tilde{f}_x(\delta)\). This function computes a rational approximation to \(x\) that agrees with it up to some arbitrary rational accuracy \(\delta\). For various reasons (for example because arithmetic on rational numbers is computationally expensive) it is desirable to use a simplified version in which (a) the accuracy \(\delta\) is of the form \(2^{p}\) for some integer \(p\) and (b) the rational function \(\tilde{f}_x(\delta)\) is replaced by an integer-valued function \(f_x(p)\) that is scaled by a factor of \(2^{-p}\) to get the rational approximation to \(x\). Concretely, the function \(f_x(p): \mathbb{Z} \to \mathbb{Z}\) satisfies

The basic building blocks of the construction are the integers: for \(x \in \mathbb{Z}\) we can take the function \(f_{x}(p) = 2^{-p} x\) as its representation. In the concrete implementation all integers are arbitrary-precision integers (JavaScript BigInt in this case).

Given two numbers \(x\) and \(y\) represented by \(f_{x}\) and \(f_{y}\) respectively and some basic arithmetic operation \(\circ\) such as addition or multiplication we need a way to construct a new function \(f_{x \circ y}\) representing \(x \circ y\). The concrete way to do this is different for each operation \(\circ\) but it always involves evaluating \(f_x\) and \(f_y\) to some (higher) precision and then combining the results with some operations on arbitrary-precision integers.

In the implementation there is an interface that all constructive numbers share:

interface CRBase {
    prec: number;
    appr: bigint;
    valid: boolean;
    id: number;
    eval: (prec: number) => bigint;
}

We store the highest precision prec and the corresponding approximation appr that were requested until now. The flag valid marks whether this approximation is valid and the id assigns an identifier to each number (this is not really needed but it was useful for quickly showing the computation as a tree above). Finally, the eval function implements \(f_{x}(p)\), i.e. it calculates an approximation to precision \(p\).

An integer is an extension of this interface that simply adds a field for the value:

interface CRInteger extends CRBase {
    type: "int";
    value: bigint;
}

The addition of two numbers on the other hand stores the two operands:

interface CRAdd extends CRBase {
    type: "add";
    x: CR;
    y: CR;
}

Here CR is any constructive real, i.e. type CR = CRInteger | CRAdd | .... The other operations work analogously.

We now need to describe how to construct a function \(f_{x \circ y}\) for some arithmetic operation \(\circ\). We also need the unary operations negation and inversion to handle subtraction and division.

For the addition, we have to compute \(f_{x}\) and \(f_{y}\) to two more levels of precision, add them and scale the result:

Here the \(\operatorname{div}\) operation is integer division.

The function for the negation \(-x\) of \(x\) is just the negation of the corresponding function:

Having access to addition and negation we can construct a function for the subtraction \(x - y\):

The function \(f_{x y}\) for the product \(x y\) is more complicated as the required precision for \(f_{x}\) depends on the value of \(y\) and vice versa. It turns out that the following formula works:

Here \(m(x)\) is the position of the most significant binary digit of \(x\). Writing \(x\) in binary and reading it left to right, this is the position of the first digit different from zero. The digit to the left of the decimal point is at position zero and the positions to the left of it are positive. For example, \(m(10) = 3\) because the binary representation of \(10\) is \(1010_{2}\).

The division operation \(x / y\) is built up from a multiplication of \(x\) with the inverse of \(y\). The function for the inverse operation is as follows:

This is scaled integer division. The denominator is an integer approximation to \(x\) with precision \(q(x) = 2 m(x) + p - 4\) which includes a factor of \(2^{q(x)}\). The first term in the numerator is \(2^{-p - q(x)}\) which has the extra factor of \(2^{-p}\) expected in \(f_{1 / x}(p)\). The second term in the numerator is there to round the integer division correctly to the nearest integer rather than simply flooring it. In his papers Boehm mentions that one could also construct the function \(f_{1 / x}(p)\) for the reciprocal using Newton-Raphson division but it does not seem like this is actually done in his implementation.

To convert a number \(x\) represented by its function \(f_{x}\) back to a string we first construct a scale factor \(r^{n}\) where \(r\) is the radix (i.e. \(r = 10\) for decimal, \(r = 2\) for binary, etc.) and \(n\) is the number of digits after the decimal point that we want to calculate. We (constructively) calculate \(x r^n\), evaluate the result with precision \(p = 0\) and convert it to a string in base \(r\). Then we insert a decimal point at the right location and add a minus sign in case the number is negative.

To convert a floating point number to a functional representation, we first extract the exponent \(e\) and mantissa \(m\). In IEE754 the exponent is stored with a bias which must be subtracted. We then store \(m\) as an integer and shift it by the exponent, \(m \ll e\). The function \(f_{x \ll s}(p)\) for a shift of \(x\) by \(s\) is

One (somewhat obvious) thing to note is that the number "0.1" converted in this way is not the same as the rational expression \(1 / 10\). The former inherits the imprecision of the floating point representation while the latter is exact.

More general functions of real numbers can be calculated from power series. For example, for the inverse tangent we have the following power series:

Some care must be taken to transform the argument to a region where the power series converges. For the \(\arctan\) one can for example use \(\arctan(x) + \arctan(x^{-1}) = \operatorname{sign}(x) \frac{\pi}{2}\). With the \(arctan\) function at hand we can calculate \(\pi\) using Machin's formula: