Centred Dyadic Approximations

There are two constructors for approximations:

• Approx is encodes some finite interval with dyadic endpoints. A real number is approximated by the approximation is it belongs to the interval.
• Bottom is the trivial approximation that approximates all real numbers.

The four fields of an Approx m e s should be thought of as:

mb

the midpoint bound, ie maximum bits available for the midpoint

m

the midpoint

e

the error term

s

the exponent

Thus, a value Approx p m e s is to be interpreted as the interval [(m-e)*2^s, (m+e)*2^s] where |m| <= 2^p.

Centred intervals

We have opted for a centred representation of the intervals. It is also possible to represent the endpoints as Dyadic numbers. The rationale for a centred repersentation is that we often normalise an approximation Approx p m e s so that e is limited in size. This allows many operations to only work on one large number m.

Potential for overflow

Since the last field (the exponent) is only an Int it may overflow. This is an optimisation that was adopted since it seems unlikely that overflow in a 64 bit Int exponent would occur. In a 32 bit system, this is potentially an issue.

The Integer data type is unbonded, but is, of course, bounded by the available memory available in the computer. No attempt has been made to check for exhausted memory.

Approximations as a Domain

Ordered by reverse inclusion the dyadic intervals encoded by the Approx approximations (including Bottom) constitute the compact elements of a Scott domain D. (It is a substructure of the (algebraic) interval domain.) We will identify our approximations with the compact elements of D.

Increasing sequences in D have suprema. A sequence converges if the length of the approximations tend to zero. The supremum of a converging sequence is a singleton set containing a real number. Let ρ be the map taking a converging sequence to the unique real number in the supremum. The computations on (computable) real numbers is via this representation map ρ.

There is no check that the sequences we have are in fact increasing, but we are assuming that all sequences are pairwise consistent. We can thus create an increasing sequence by considering the sequence of finite suprema. For correctness, we have to ensure that all operations done on consistent sequences result in consistent sequences. If non-consistent sequences are somehow input we can make no guarantees at all about the computed value.

Note, that we cannot ensure that converging sequences are mapped to converging sequences because of properties of computable real arithmetic. In particular, at any discuntinuity, it is impossible to compute a converging sequence.

The Computable Real data type

Computable reals are realised as infinite sequences of centred dyadic representations.

All approximations in such a sequence should be pairwise consistent, i.e., have a non-empty intersection. However, there is no check that this is actually the case.

If the diameter of the approximations tend to zero we say that the sequences converges to the unique real number in the intersection of all intervals. Given the domain D of approximations described in Approx, we have a representation (a retraction) ρ from the converging sequences in D to ℝ. Some operations on computable reals are partial, notably equality and rnf m seq there is no guarantee that a rnf m seq rnf m seq rnf m seq there is a bound on how much effort is rnf m seq mation. For involved computations it is rnf m seq proximations are trivial, i.e., rnf m seq ually converge, it will generate proper rnf m seq initial trivial approximations. rnf m seq The amount of added effort in each iteration is rather substantial so the expected precision of approximations increase very quickly.

The actual data type

In fact, the type CR is a newtype of ZipList Approx in the implementation of infinite sequences of approximations, as that allows for applicative style. Hopefully, it is not needed to access the internal representation of CR directly.

A type synonym. Used to denote number of bits after binary point.

Gives a decimal representation of an approximation. It tries to give as many decimal digits as possible given the precision of the approximation. The representation may be wrong by 1 ulp (unit in last place). If the value is not exact the representation will be followed by ~.

The representation is not always intuitive:

>>> showA (Approx 1 1 0)
"1.~"

The meaning of the above is that it is 1, but then the added ~ (which must be after the decimal point) means that the last position may be off by 1, i.e., it could be down to 0 or up to 2. And [0,2] is indeed the range encoded by the above approximation.

Similar to showA but can generate representations in other bases (<= 16).

Give the "bound on midpoint bit-size" component of an Approx.

The midpoint coponent should always be bounded by this as follows: abs m <= 2^mb.

Construct a centred approximation from the end-points.

Gives the lower bound of an approximation as an Extended Dyadic number.

Gives the upper bound of an approximation as an Extended Dyadic number.

Gives the lower bound of an Approx as an exact Approx.

Gives the upper bound of an Approx as an exact Approx.

Gives the mid-point of an approximation as a Maybe Dyadic number.

Gives the centre of an Approx as an exact Approx.

Gives the radius of an approximation as a Dyadic number. Currently a partial function. Should be made to return an Extended Dyadic.

Gives the lower bound of an approximation as an Extended Dyadic number.

Returns True if the approximation is exact, i.e., it's diameter is 0.

Checks if a number is approximated by an approximation, i.e., if it belongs to the interval encoded by the approximation.

A partial order on approximations. The first approximation is better than the second if it is a sub-interval of the second.

Turns a Dyadic number into an exact approximation.

Turns a Dyadic number into an exact approximation.

Convert a rational number into an approximation of that number with Precision bits correct after the binary point.

Convert a rational number into an approximation of that number with mBound significant bits correct.

Compute the reciprocal of an approximation. The number of bits after the binary point is bounded by the 'midpoint bound' if the input is exact. Otherwise, a good approximation with essentially the same significance as the input will be computed.

Divide an approximation by an integer.

Compute the modulus of two approximations.

Compute the integer quotient (although returned as an Approx since it may be necessary to return Bottom if the integer quotient can't be determined) and the modulus as an approximation of two approximations.

Convert the centre of an approximation into a Maybe Double.

Computes the precision of an approximation. This is roughly the number of correct bits after the binary point.

Computes the significance of an approximation. This is roughly the number of significant bits.

This function bounds the error term of an Approx.

It is always the case that x `better` boundErrorTerm x.

Consider an approximation Approx mb m e s. If e has k bits then that essentially expresses that the last k bits of m are unknown or garbage. By scaling both m and e so that e has a small number of bits we save on memory space and computational effort to compute operations. On the other hand, if we remove too many bits in this way, the shift in the mid-point of the interval becomes noticable and it may adversely affect convergence speed of computations. The number of bits allowed for e after the operation is determined by the constant errorBits.

Domain interpretation and verification

For this implementation to be correct it is required that this function is below the identity on the domain D of Approx. For efficiency it is desirable to be as close to the identity as possible.

This function will map a converging sequence to a converging sequence.

Limits the size of an approximation by restricting how much precision an approximation can have.

It is always the case that x `better` limitSize x.

This is accomplished by restricting the exponent of the approximation from below. In other words, we limit the precision possible.

It is conceivable to limit the significance of an approximation rather than the precision. This would be an interesting research topic.

Domain interpretation and verification

For this implementation to be correct it is required that this function is below the identity on the domain D of Approx. For efficiency it is desirable to be as close to the identity as possible.

This function will NOT map a converging sequence to a converging sequence for a fixed precision argument. However, if the function is applied with increasing precision for a converging sequence, then this will give a converging sequence.

Throws an exception if the precision of an approximation is not larger than the deafult minimum.

Bounds the error term and limits the precision of an approximation.

It is always the case that x `better` limitAndBound x.

Find the hull of two approximations.

Find the intersection of two approximations.

Determine if two approximations overlap.

Given a list of polynom coefficients and a value this evaluates the polynomial at that value.

This works correctly only for exact coefficients.

Should give a tighter bound on the result since we reduce the dependency problem.

Gives a list of powers of a number, i.e., [1,x,x^2,...].

Computes powers of approximations. Should give tighter intervals than the general pow since take the dependency problem into account. However, so far benchmarking seems to indicate that the cost is too high, but this may depend on the application.

Old implementation of sqrt using Heron's method. Using the current method below proved to be more than twice as fast for small arguments (~50 bits) and many times faster for larger arguments.

Compute the square root of an approximation.

This and many other operations on approximations is just a reimplementation of interval arithmetic, with an extra argument limiting the effort put into the computation. This is done via the precision argument.

The resulting approximation should approximate the image of every point in the input approximation.

This uses Newton's method for computing the reciprocal of the square root.

The starting values for newton iterations can be found using the auxiliary function findStartingValues below.

For example, to generate the starting values for sqrtRecD above using three leading bits for both odd and even exponents the following was used:

findStartingValues ((1/) . sqrt) [1,1.125..2]

Approx 4172150648 1 (-32),Approx 3945434766 1 (-32),Approx 3752147976 1 (-32),Approx 3584793264 1 (-32),Approx 3438037830 1 (-32),Approx 3307969824 1 (-32),Approx 3191645366 1 (-32),Approx 3086800564 1 (-32)

> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it 0.111110~ 0.111011~ 0.111000~ 0.110101~ 0.110011~ 0.110001~ 0.110000~ 0.101110~ > findStartingValues ((1/) . sqrt) [2,2.25..4]

Approx 2950156016 1 (-32),Approx 2789843678 1 (-32),Approx 2653169278 1 (-32),Approx 2534831626 1 (-32),Approx 2431059864 1 (-32),Approx 2339087894 1 (-32),Approx 2256834080 1 (-32),Approx 2182697612 1 (-32)

> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it 0.101100~ 0.101010~ 0.101000~ 0.100110~ 0.100100~ 0.100011~ 0.100010~ 0.100001~

Computes the square root of a Dyadic number to the specified base. The Newton-Raphson method may overestimates the square root, but the overestimate is bounded by 1 ulp. For example, sqrtD 0 2 will give 2, whereas the closest integer to the square root is 1. Need double precision to guarantee correct rounding, which was not considered worthwhile.

This is actually Heron's method and is no longer used in Approx as it is faster to use sqrtRecD.

Shift a dyadic number to a given base and round in case of right shifts.

Square an approximation. Gives the exact image interval, as opposed to multiplicating a number with itself which will give a slightly larger interval due to the dependency problem.

Computes the list [lg 0!, lg 1!, lg 2!, ...].

Computes the sum of the form ∑ aₙxⁿ where aₙ and x are approximations.

Terms are added as long as they are larger than the current precision bound. The sum is adjusted for the tail of the series. For this to be correct we need the the terms to converge geometrically to 0 by a factor of at least 2.

The exponential of an approximation. There are three implementation using Taylor expansion here. This is just choosing between them.

More thorough benchmarking would be desirable.

Is faster for small approximations < ~2000 bits.

Exponential by binary splitting summation of Taylor series.

Exponential by summation of Taylor series.

Exponential by summation of Taylor series.

Computing the logarithm of an approximation. This chooses the fastest implementation.

More thorough benchmarking is desirable.

Binary splitting is faster than Taylor. AGM should be used over ~1000 bits.

Logarithm by binary splitting summation of Taylor series.

Logarithm by summation of Taylor series.

Computes sine by summation of Taylor series after two levels of range reductions.

First level of range reduction for sine. Brings it into the interval [-π2,π2].

Second level of range reduction for sine.

Computes the sine of an approximation. Chooses the best implementation.

Computes the cosine of an approximation. Chooses the best implementation.

Computes the arc tangent of an approximation. Chooses the best implementation.

Computes the sine of an approximation by binary splitting summation of Taylor series.

Begins by reducing the interval to [0,π/4].

Computes the cosine of an approximation by binary splitting summation of Taylor series.

Begins by reducing the interval to [0,π/4].

Computes the arc tangent of an approximation by binary splitting summation of Taylor series.

Computes a sequence of approximations of π using binary splitting summation of Ramanujan's series. See Haible and Papanikolaou 1998.

Computes an Approx of π of the given precision.

Compute π using Machin's formula. Lifted from computation on dyadic numbers.

Compute π using AGM as described in Borwein and Borwein's book 'Pi and the AGM'. Lifted from computation on dyadic numbers.

Compute π using AGM as described in Borwein and Borwein's book 'Pi and the AGM'.

Compute approximations of ln 2. Lifted from computation on dyadic numbers.

Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and the AGM'. An approximation of pi is produced as a by-product.

Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and the AGM'. TODO: adapt to mBound

Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and the AGM'.

Shows the internal representation of a CR. The first n approximations are shown on separate lines.

Shows a CR with the desired precision.

Check that an approximation has at least d bits of precision. This is used to bail out of computations where the size of approximation grow quickly.

Given a CR this functions finds an approximation of that number to within the precision required.

Gives a Double approximation of a CR number.

Evaluate a polynomial, given as a list of its coefficients, at the given point.

Computes the sum of a Taylor series, given as a list of its coefficients, at the given point.

π computed using Machin's formula. Computed directly on CR.

π computed using Machin's formula. Computed on Approx approximations.

π computed using Borwein's formula. Computed on Approx approximations.

π computed using binary splitting. Computed on Approx approximations.

The constant ln 2.

The sine function computed using Taylor's series. Computed directly on CR. Will have poor behaviour on larger inputs as no range reduction is performed.

The sine function computed using Taylor's series. Computed directly on CR.

The cosine function computed using Taylor's series. Computed directly on CR.

The square root function. Lifted from square root application on Approx approximations.

The exponential computed using Taylor's series. Computed directly on CR. Will have poor behaviour on larger inputs as no range reduction is performed.