Safe Haskell  None 

Language  Haskell2010 
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
 data Approx
 newtype CR = CR {}
 type Precision = Int
 showA :: Approx > String
 showInBaseA :: Int > Approx > String
 mBound :: Approx > Int
 mapMB :: (Int > Int) > Approx > Approx
 setMB :: Int > Approx > Approx
 enforceMB :: Approx > Approx
 approxAutoMB :: Integer > Integer > Int > Approx
 approxMB :: Int > Integer > Integer > Int > Approx
 approxMB2 :: Int > Int > Integer > Integer > Int > Approx
 endToApprox :: Int > Extended Dyadic > Extended Dyadic > Approx
 lowerBound :: Approx > Extended Dyadic
 upperBound :: Approx > Extended Dyadic
 lowerA :: Approx > Approx
 upperA :: Approx > Approx
 centre :: Approx > Maybe Dyadic
 centreA :: Approx > Approx
 radius :: Approx > Extended Dyadic
 diameter :: Approx > Extended Dyadic
 exact :: Approx > Bool
 approximatedBy :: Real a => a > Approx > Bool
 better :: Approx > Approx > Bool
 fromDyadic :: Dyadic > Approx
 fromDyadicMB :: Int > Dyadic > Approx
 toApprox :: Precision > Rational > Approx
 toApproxMB :: Int > Rational > Approx
 recipA :: Approx > Approx
 divAInteger :: Approx > Integer > Approx
 modA :: Approx > Approx > Approx
 divModA :: Approx > Approx > (Approx, Approx)
 toDoubleA :: Approx > Maybe Double
 precision :: Approx > Extended Precision
 significance :: Approx > Extended Int
 boundErrorTerm :: Approx > Approx
 boundErrorTermMB :: Approx > Approx
 limitSize :: Precision > Approx > Approx
 checkPrecisionLeft :: Approx > Approx
 limitAndBound :: Precision > Approx > Approx
 limitAndBoundMB :: Precision > Approx > Approx
 unionA :: Approx > Approx > Approx
 intersectionA :: Approx > Approx > Approx
 consistentA :: Approx > Approx > Bool
 poly :: [Approx] > Approx > Approx
 pow :: Num a => a > [a]
 powers :: Approx > [Approx]
 sqrtHeronA :: Precision > Approx > Approx
 sqrtA :: Approx > Approx
 sqrtRecA :: Precision > Approx > Approx
 findStartingValues :: (Double > Double) > [Double] > [Approx]
 sqrtD :: Int > Dyadic > Dyadic
 shiftD :: Int > Dyadic > Dyadic
 sqrA :: Approx > Approx
 log2Factorials :: [Int]
 taylorA :: Precision > [Approx] > Approx > Approx
 expA :: Approx > Approx
 expBinarySplittingA :: Precision > Approx > Approx
 expTaylorA :: Precision > Approx > Approx
 expTaylorA' :: Approx > Approx
 logA :: Approx > Approx
 logInternal :: Approx > Approx
 logBinarySplittingA :: Precision > Approx > Approx
 logTaylorA :: Precision > Approx > Approx
 sinTaylorA :: Approx > Approx
 sinTaylorRed1A :: Approx > (Approx, (Maybe Approx, Maybe Approx))
 sinTaylorRed2A :: Approx > Approx
 sinA :: Approx > Approx
 cosA :: Approx > Approx
 atanA :: Precision > Approx > Approx
 sinBinarySplittingA :: Precision > Approx > Approx
 cosBinarySplittingA :: Precision > Approx > Approx
 atanBinarySplittingA :: Precision > Approx > Approx
 atanTaylorA :: Precision > Approx > Approx
 piRaw :: [Approx]
 piA :: Precision > Approx
 piMachinA :: Precision > Approx
 piBorweinA :: Precision > Approx
 piAgmA :: Precision > Approx > Approx
 log2A :: Precision > Approx
 lnSuperSizeUnknownPi :: Precision > Approx > (Approx, Approx)
 logAgmA :: Precision > Approx > Approx
 agmLn :: Precision > Approx > Approx
 showCRN :: Int > CR > String
 showCR :: Int > CR > String
 ok :: Int > Approx > Approx
 require :: Int > CR > Approx
 toDouble :: CR > Maybe Double
 fromDouble :: Double > CR
 fromDoubleAsExactValue :: Double > CR
 polynomial :: [CR] > CR > CR
 taylorCR :: [CR] > CR > CR
 atanCR :: CR > CR
 piCRMachin :: CR
 piMachinCR :: CR
 piBorweinCR :: CR
 piBinSplitCR :: CR
 ln2 :: CR
 sinCRTaylor :: CR > CR
 sinCR :: CR > CR
 cosCR :: CR > CR
 sqrtCR :: CR > CR
 expCR :: CR > CR
Documentation
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
[(me)*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 nonconsistent 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
Enum Approx Source #  Not a sensible instance. Just used to allow to allow enumerating integers using '..' notation. 
Defined in Data.CDAR.Approx  
Eq Approx Source #  Two approximations are equal if they encode the same interval. 
Fractional Approx Source #  Not a proper Fractional type as Approx are intervals. 
Num Approx Source #  
Ord Approx Source #  Not a proper Ord type as Approx are intervals. 
Read Approx Source #  
Real Approx Source #  The 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. 
Defined in Data.CDAR.Approx toRational :: Approx > Rational #  
Show Approx Source #  
NFData Approx Source #  
Defined in Data.CDAR.Approx  
Scalable Approx 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 nonempty 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.
Instances
Eq CR Source #  
Floating CR Source #  
Defined in Data.CDAR.Approx  
Fractional CR Source #  
Num CR Source #  
Ord CR Source #  
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

Real CR Source #  
Defined in Data.CDAR.Approx toRational :: CR > Rational #  
Scalable CR Source #  
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 bitsize" component of an Approx
.
The midpoint coponent should always be bounded by this as follows:
abs m <= 2^mb
.
endToApprox :: Int > Extended Dyadic > Extended Dyadic > Approx Source #
Construct a centred approximation from the endpoints.
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 subinterval of the second.
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.
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 midpoint 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.
boundErrorTermMB :: Approx > Approx Source #
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
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.
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 NewtonRaphson 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' :: 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.
logInternal :: Approx > Approx Source #
logBinarySplittingA :: Precision > Approx > Approx Source #
Logarithm by binary splitting 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.
Computes a sequence of approximations of π using binary splitting summation of Ramanujan's series. See Haible and Papanikolaou 1998.
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 byproduct.
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.
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.
fromDouble :: Double > CR Source #
fromDoubleAsExactValue :: Double > CR Source #
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.
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.
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.