cdar-mBound-0.1.0.1: Exact real arithmetic using Centred Dyadic Approximations
Safe HaskellNone
LanguageHaskell2010

Data.CDAR.Approx

Description

Computable Real Arithmetic

This module provides the data type CR that implements the real closed field of computable real numbers.

Centred Dyadic Approximations

The computable reals are realised as lists of rapidly shrinking intervals. The intervals used here are centred dyadic intervals, implemented here as the data type Approx.

For more information on the theoretical aspects see http://cs.swan.ac.uk/~csjens/pdf/centred.pdf.

Synopsis

Documentation

data Approx Source #

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.

Instances

Instances details
Enum Approx Source #

Not a sensible instance. Just used to allow to allow enumerating integers using '..' notation.

Instance details

Defined in Data.CDAR.Approx

Eq Approx Source #

Two approximations are equal if they encode the same interval.

Instance details

Defined in Data.CDAR.Approx

Methods

(==) :: Approx -> Approx -> Bool #

(/=) :: Approx -> Approx -> Bool #

Fractional Approx Source #

Not a proper Fractional type as Approx are intervals.

Instance details

Defined in Data.CDAR.Approx

Num Approx Source # 
Instance details

Defined in Data.CDAR.Approx

Ord Approx Source #

Not a proper Ord type as Approx are intervals.

Instance details

Defined in Data.CDAR.Approx

Read Approx Source # 
Instance details

Defined in Data.CDAR.Approx

Real Approx Source #

The toRational function is partial since there is no good rational number to return for the trivial approximation Bottom.

Note that converting to a rational number will give only a single rational point. Do not expect to be able to recover the interval from this value.

Instance details

Defined in Data.CDAR.Approx

Show Approx Source # 
Instance details

Defined in Data.CDAR.Approx

NFData Approx Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

rnf :: Approx -> () #

Scalable Approx Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

scale :: Approx -> Int -> Approx Source #

newtype CR Source #

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.

Constructors

CR 

Fields

Instances

Instances details
Eq CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

(==) :: CR -> CR -> Bool #

(/=) :: CR -> CR -> Bool #

Floating CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

pi :: CR #

exp :: CR -> CR #

log :: CR -> CR #

sqrt :: CR -> CR #

(**) :: CR -> CR -> CR #

logBase :: CR -> CR -> CR #

sin :: CR -> CR #

cos :: CR -> CR #

tan :: CR -> CR #

asin :: CR -> CR #

acos :: CR -> CR #

atan :: CR -> CR #

sinh :: CR -> CR #

cosh :: CR -> CR #

tanh :: CR -> CR #

asinh :: CR -> CR #

acosh :: CR -> CR #

atanh :: CR -> CR #

log1p :: CR -> CR #

expm1 :: CR -> CR #

log1pexp :: CR -> CR #

log1mexp :: CR -> CR #

Fractional CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

(/) :: CR -> CR -> CR #

recip :: CR -> CR #

fromRational :: Rational -> CR #

Num CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

(+) :: CR -> CR -> CR #

(-) :: CR -> CR -> CR #

(*) :: CR -> CR -> CR #

negate :: CR -> CR #

abs :: CR -> CR #

signum :: CR -> CR #

fromInteger :: Integer -> CR #

Ord CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

compare :: CR -> CR -> Ordering #

(<) :: CR -> CR -> Bool #

(<=) :: CR -> CR -> Bool #

(>) :: CR -> CR -> Bool #

(>=) :: CR -> CR -> Bool #

max :: CR -> CR -> CR #

min :: CR -> CR -> CR #

Read CR Source #

Reads a floating point representation of a real number and interprets that as a CR. Does not currently allow for the same format output by showCR.

Instance details

Defined in Data.CDAR.Approx

Real CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

toRational :: CR -> Rational #

Scalable CR Source # 
Instance details

Defined in Data.CDAR.Approx

Methods

scale :: CR -> Int -> CR Source #

type Precision = Int Source #

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

showA :: Approx -> String Source #

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.

showInBaseA :: Int -> Approx -> String Source #

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

mBound :: Approx -> Int Source #

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.

mapMB :: (Int -> Int) -> Approx -> Approx Source #

endToApprox :: Int -> Extended Dyadic -> Extended Dyadic -> Approx Source #

Construct a centred approximation from the end-points.

lowerBound :: Approx -> Extended Dyadic Source #

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

upperBound :: Approx -> Extended Dyadic Source #

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

lowerA :: Approx -> Approx Source #

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

upperA :: Approx -> Approx Source #

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

centre :: Approx -> Maybe Dyadic Source #

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

centreA :: Approx -> Approx Source #

Gives the centre of an Approx as an exact Approx.

radius :: Approx -> Extended Dyadic Source #

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

diameter :: Approx -> Extended Dyadic Source #

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

exact :: Approx -> Bool Source #

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

approximatedBy :: Real a => a -> Approx -> Bool Source #

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

better :: Approx -> Approx -> Bool Source #

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

fromDyadic :: Dyadic -> Approx Source #

Turns a Dyadic number into an exact approximation.

fromDyadicMB :: Int -> Dyadic -> Approx Source #

Turns a Dyadic number into an exact approximation.

toApprox :: Precision -> Rational -> Approx Source #

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

toApproxMB :: Int -> Rational -> Approx Source #

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

recipA :: Approx -> Approx Source #

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.

divAInteger :: Approx -> Integer -> Approx Source #

Divide an approximation by an integer.

modA :: Approx -> Approx -> Approx Source #

Compute the modulus of two approximations.

divModA :: Approx -> Approx -> (Approx, Approx) Source #

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.

toDoubleA :: Approx -> Maybe Double Source #

Convert the centre of an approximation into a Maybe Double.

precision :: Approx -> Extended Precision Source #

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

significance :: Approx -> Extended Int Source #

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

boundErrorTerm :: Approx -> Approx Source #

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.

limitSize :: Precision -> Approx -> Approx Source #

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.

checkPrecisionLeft :: Approx -> Approx Source #

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

limitAndBound :: Precision -> Approx -> Approx Source #

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

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

unionA :: Approx -> Approx -> Approx Source #

Find the hull of two approximations.

intersectionA :: Approx -> Approx -> Approx Source #

Find the intersection of two approximations.

consistentA :: Approx -> Approx -> Bool Source #

Determine if two approximations overlap.

poly :: [Approx] -> Approx -> Approx Source #

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.

pow :: Num a => a -> [a] Source #

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

powers :: Approx -> [Approx] Source #

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.

sqrtHeronA :: Precision -> Approx -> Approx Source #

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.

sqrtA :: Approx -> Approx Source #

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.

sqrtRecA :: Precision -> Approx -> Approx Source #

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

findStartingValues :: (Double -> Double) -> [Double] -> [Approx] Source #

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~

sqrtD :: Int -> Dyadic -> Dyadic Source #

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.

shiftD :: Int -> Dyadic -> Dyadic Source #

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

sqrA :: Approx -> Approx Source #

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.

log2Factorials :: [Int] Source #

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

taylorA :: Precision -> [Approx] -> Approx -> Approx Source #

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.

expA :: Approx -> Approx Source #

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.

expBinarySplittingA :: Precision -> Approx -> Approx Source #

Exponential by binary splitting summation of Taylor series.

expTaylorA :: Precision -> Approx -> Approx Source #

Exponential by summation of Taylor series.

expTaylorA' :: Approx -> Approx Source #

Exponential by summation of Taylor series.

logA :: Approx -> Approx Source #

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.

logBinarySplittingA :: Precision -> Approx -> Approx Source #

Logarithm by binary splitting summation of Taylor series.

logTaylorA :: Precision -> Approx -> Approx Source #

Logarithm by summation of Taylor series.

sinTaylorA :: Approx -> Approx Source #

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

sinTaylorRed1A :: Approx -> (Approx, (Maybe Approx, Maybe Approx)) Source #

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

sinTaylorRed2A :: Approx -> Approx Source #

Second level of range reduction for sine.

sinA :: Approx -> Approx Source #

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

cosA :: Approx -> Approx Source #

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

atanA :: Precision -> Approx -> Approx Source #

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

sinBinarySplittingA :: Precision -> Approx -> Approx Source #

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

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

cosBinarySplittingA :: Precision -> Approx -> Approx Source #

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

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

atanBinarySplittingA :: Precision -> Approx -> Approx Source #

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

piRaw :: [Approx] Source #

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

piA :: Precision -> Approx Source #

Computes an Approx of π of the given precision.

piMachinA :: Precision -> Approx Source #

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

piBorweinA :: Precision -> Approx Source #

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

piAgmA :: Precision -> Approx -> Approx Source #

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

log2A :: Precision -> Approx Source #

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

lnSuperSizeUnknownPi :: Precision -> Approx -> (Approx, Approx) Source #

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.

logAgmA :: Precision -> Approx -> Approx Source #

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

agmLn :: Precision -> Approx -> Approx Source #

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

showCRN :: Int -> CR -> String Source #

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

showCR :: Int -> CR -> String Source #

Shows a CR with the desired precision.

ok :: Int -> Approx -> Approx Source #

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.

require :: Int -> CR -> Approx Source #

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

toDouble :: CR -> Maybe Double Source #

Gives a Double approximation of a CR number.

polynomial :: [CR] -> CR -> CR Source #

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

taylorCR :: [CR] -> CR -> CR Source #

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

piCRMachin :: CR Source #

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

piMachinCR :: CR Source #

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

piBorweinCR :: CR Source #

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

piBinSplitCR :: CR Source #

π computed using binary splitting. Computed on Approx approximations.

ln2 :: CR Source #

The constant ln 2.

sinCRTaylor :: CR -> CR Source #

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.

sinCR :: CR -> CR Source #

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

cosCR :: CR -> CR Source #

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

sqrtCR :: CR -> CR Source #

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

expCR :: CR -> CR Source #

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