Safe Haskell | None |
---|---|

Language | Haskell98 |

Exact Real Arithmetic - Computable reals. Inspired by ''The most unreliable technique for computing pi.'' See also http://www.haskell.org/haskellwiki/Exact_real_arithmetic .

- type T = (Exponent, Mantissa)
- type FixedPoint = (Integer, Mantissa)
- type Mantissa = [Digit]
- type Digit = Int
- type Exponent = Int
- type Basis = Int
- moveToZero :: Basis -> Digit -> (Digit, Digit)
- checkPosDigit :: Basis -> Digit -> Digit
- checkDigit :: Basis -> Digit -> Digit
- nonNegative :: Basis -> T -> T
- nonNegativeMant :: Basis -> Mantissa -> Mantissa
- compress :: Basis -> T -> T
- compressFirst :: Basis -> T -> T
- compressMant :: Basis -> Mantissa -> Mantissa
- compressSecondMant :: Basis -> Mantissa -> Mantissa
- prependDigit :: Basis -> T -> T
- trim :: T -> T
- trimUntil :: Exponent -> T -> T
- trimOnce :: T -> T
- decreaseExp :: Basis -> T -> T
- pumpFirst :: Basis -> Mantissa -> Mantissa
- decreaseExpFP :: Basis -> (Exponent, FixedPoint) -> (Exponent, FixedPoint)
- pumpFirstFP :: Basis -> FixedPoint -> FixedPoint
- negativeExp :: Basis -> T -> T
- fromBaseCardinal :: Basis -> Integer -> T
- fromBaseInteger :: Basis -> Integer -> T
- mantissaToNum :: C a => Basis -> Mantissa -> a
- mantissaFromCard :: C a => Basis -> a -> Mantissa
- mantissaFromInt :: C a => Basis -> a -> Mantissa
- mantissaFromFixedInt :: Basis -> Exponent -> Integer -> Mantissa
- fromBaseRational :: Basis -> Rational -> T
- toFixedPoint :: Basis -> T -> FixedPoint
- fromFixedPoint :: Basis -> FixedPoint -> T
- toDouble :: Basis -> T -> Double
- fromDouble :: Basis -> Double -> T
- fromDoubleApprox :: Basis -> Double -> T
- fromDoubleRough :: Basis -> Double -> T
- mantLengthDouble :: Basis -> Exponent
- liftDoubleApprox :: Basis -> (Double -> Double) -> T -> T
- liftDoubleRough :: Basis -> (Double -> Double) -> T -> T
- showDec :: Exponent -> T -> String
- showHex :: Exponent -> T -> String
- showBin :: Exponent -> T -> String
- showBasis :: Basis -> Exponent -> T -> String
- powerBasis :: Basis -> Exponent -> T -> T
- rootBasis :: Basis -> Exponent -> T -> T
- fromBasis :: Basis -> Basis -> T -> T
- fromBasisMant :: Basis -> Basis -> Mantissa -> Mantissa
- cmp :: Basis -> T -> T -> Ordering
- lessApprox :: Basis -> Exponent -> T -> T -> Bool
- trunc :: Exponent -> T -> T
- equalApprox :: Basis -> Exponent -> T -> T -> Bool
- ifLazy :: Basis -> Bool -> T -> T -> T
- mean2 :: Basis -> (Digit, Digit) -> (Digit, Digit) -> Digit
- withTwoMantissas :: Mantissa -> Mantissa -> a -> ((Digit, Mantissa) -> (Digit, Mantissa) -> a) -> a
- align :: Basis -> Exponent -> T -> T
- alignMant :: Basis -> Exponent -> T -> Mantissa
- absolute :: T -> T
- absMant :: Mantissa -> Mantissa
- fromLaurent :: T Digit -> T
- toLaurent :: T -> T Digit
- liftLaurent2 :: (T Digit -> T Digit -> T Digit) -> T -> T -> T
- liftLaurentMany :: ([T Digit] -> T Digit) -> [T] -> T
- add :: Basis -> T -> T -> T
- sub :: Basis -> T -> T -> T
- neg :: Basis -> T -> T
- addSome :: Basis -> [T] -> T
- addMany :: Basis -> [T] -> T
- type Series = [(Exponent, T)]
- series :: Basis -> Series -> T
- seriesPlain :: Basis -> Series -> T
- splitAtPadZero :: Int -> Mantissa -> (Mantissa, Mantissa)
- splitAtMatchPadZero :: [()] -> Mantissa -> (Mantissa, Mantissa)
- truncSeriesSummands :: Series -> Series
- scale :: Basis -> Digit -> T -> T
- scaleSimple :: Basis -> T -> T
- scaleMant :: Basis -> Digit -> Mantissa -> Mantissa
- mulSeries :: Basis -> T -> T -> Series
- mul :: Basis -> T -> T -> T
- divide :: Basis -> T -> T -> T
- divMant :: Basis -> Mantissa -> Mantissa -> Mantissa
- divMantSlow :: Basis -> Mantissa -> Mantissa -> Mantissa
- reciprocal :: Basis -> T -> T
- divIntMant :: Basis -> Digit -> Mantissa -> Mantissa
- divIntMantInf :: Basis -> Digit -> Mantissa -> Mantissa
- divInt :: Basis -> Digit -> T -> T
- sqrt :: Basis -> T -> T
- sqrtDriver :: Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T
- sqrtMant :: Basis -> Mantissa -> Mantissa
- sqrtFP :: Basis -> FixedPoint -> Mantissa
- sqrtNewton :: Basis -> T -> T
- sqrtFPNewton :: Basis -> FixedPoint -> Mantissa
- lazyInits :: [a] -> [[a]]
- expSeries :: Basis -> T -> Series
- expSmall :: Basis -> T -> T
- expSeriesLazy :: Basis -> T -> Series
- expSmallLazy :: Basis -> T -> T
- exp :: Basis -> T -> T
- intPower :: Basis -> Integer -> T -> T -> T -> T
- cardPower :: Basis -> Integer -> T -> T -> T
- powerSeries :: Basis -> Rational -> T -> Series
- powerSmall :: Basis -> Rational -> T -> T
- power :: Basis -> Rational -> T -> T
- root :: Basis -> Integer -> T -> T
- cosSinhSmall :: Basis -> T -> (T, T)
- cosSinSmall :: Basis -> T -> (T, T)
- cosSinFourth :: Basis -> T -> (T, T)
- cosSin :: Basis -> T -> (T, T)
- tan :: Basis -> T -> T
- cot :: Basis -> T -> T
- lnSeries :: Basis -> T -> Series
- lnSmall :: Basis -> T -> T
- lnNewton :: Basis -> T -> T
- lnNewton' :: Basis -> T -> T
- ln :: Basis -> T -> T
- angle :: Basis -> (T, T) -> T
- arctanSeries :: Basis -> T -> Series
- arctanSmall :: Basis -> T -> T
- arctanStem :: Basis -> Digit -> T
- arctan :: Basis -> T -> T
- arctanClassic :: Basis -> T -> T
- zero :: T
- one :: T
- minusOne :: T
- eConst :: Basis -> T
- recipEConst :: Basis -> T
- piConst :: Basis -> T
- sliceVertPair :: [a] -> [(a, a)]

# types

type FixedPoint = (Integer, Mantissa) Source

# basic helpers

checkPosDigit :: Basis -> Digit -> Digit Source

checkDigit :: Basis -> Digit -> Digit Source

nonNegative :: Basis -> T -> T Source

Converts all digits to non-negative digits, that is the usual positional representation. However the conversion will fail when the remaining digits are all zero. (This cannot be improved!)

nonNegativeMant :: Basis -> Mantissa -> Mantissa Source

Requires, that no digit is `(basis-1)`

or `(1-basis)`

.
The leading digit might be negative and might be `-basis`

or `basis`

.

compressFirst :: Basis -> T -> T Source

Compress first digit. May prepend a digit.

compressMant :: Basis -> Mantissa -> Mantissa Source

Does not prepend a digit.

compressSecondMant :: Basis -> Mantissa -> Mantissa Source

Compress second digit. Sometimes this is enough to keep the digits in the admissible range. Does not prepend a digit.

prependDigit :: Basis -> T -> T Source

decreaseExp :: Basis -> T -> T Source

Accept a high leading digit for the sake of a reduced exponent.
This eliminates one leading digit.
Like `pumpFirst`

but with exponent management.

pumpFirst :: Basis -> Mantissa -> Mantissa Source

Merge leading and second digit.
This is somehow an inverse of `compressMant`

.

decreaseExpFP :: Basis -> (Exponent, FixedPoint) -> (Exponent, FixedPoint) Source

pumpFirstFP :: Basis -> FixedPoint -> FixedPoint Source

negativeExp :: Basis -> T -> T Source

Make sure that a number with absolute value less than 1 has a (small) negative exponent. Also works with zero because it chooses an heuristic exponent for stopping.

# conversions

## integer

fromBaseCardinal :: Basis -> Integer -> T Source

fromBaseInteger :: Basis -> Integer -> T Source

mantissaToNum :: C a => Basis -> Mantissa -> a Source

mantissaFromCard :: C a => Basis -> a -> Mantissa Source

mantissaFromInt :: C a => Basis -> a -> Mantissa Source

## rational

fromBaseRational :: Basis -> Rational -> T Source

## fixed point

toFixedPoint :: Basis -> T -> FixedPoint Source

Split into integer and fractional part.

fromFixedPoint :: Basis -> FixedPoint -> T Source

## floating point

fromDouble :: Basis -> Double -> T Source

cf. `floatToDigits`

fromDoubleApprox :: Basis -> Double -> T Source

Only return as much digits as are contained in Double. This will speedup further computations.

fromDoubleRough :: Basis -> Double -> T Source

mantLengthDouble :: Basis -> Exponent Source

## text

## basis

powerBasis :: Basis -> Exponent -> T -> T Source

Convert from a `b`

basis representation to a `b^e`

basis.

Works well with every exponent.

rootBasis :: Basis -> Exponent -> T -> T Source

Convert from a `b^e`

basis representation to a `b`

basis.

Works well with every exponent.

fromBasis :: Basis -> Basis -> T -> T Source

Convert between arbitrary bases. This conversion is expensive (quadratic time).

# comparison

cmp :: Basis -> T -> T -> Ordering Source

The basis must be at least ***. Note: Equality cannot be asserted in finite time on infinite precise numbers. If you want to assert, that a number is below a certain threshold, you should not call this routine directly, because it will fail on equality. Better round the numbers before comparison.

ifLazy :: Basis -> Bool -> T -> T -> T Source

If all values are completely defined, then it holds

if b then x else y == ifLazy b x y

However if `b`

is undefined,
then it is at least known that the result is between `x`

and `y`

.

mean2 :: Basis -> (Digit, Digit) -> (Digit, Digit) -> Digit Source

mean2 b (x0,x1) (y0,y1)

computes ` round ((x0.x1 + y0.y1)/2) `

,
where `x0.x1`

and `y0.y1`

are positional rational numbers
with respect to basis `b`

withTwoMantissas :: Mantissa -> Mantissa -> a -> ((Digit, Mantissa) -> (Digit, Mantissa) -> a) -> a Source

alignMant :: Basis -> Exponent -> T -> Mantissa Source

Get the mantissa in such a form that it fits an expected exponent.

`x`

and `(e, alignMant b e x)`

represent the same number.

# arithmetic

fromLaurent :: T Digit -> T Source

addSome :: Basis -> [T] -> T Source

Add at most `basis`

summands.
More summands will violate the allowed digit range.

addMany :: Basis -> [T] -> T Source

Add many numbers efficiently by computing sums of sub lists with only little carry propagation.

series :: Basis -> Series -> T Source

Add an infinite number of numbers. You must provide a list of estimate of the current remainders. The estimates must be given as exponents of the remainder. If such an exponent is too small, the summation will be aborted. If exponents are too big, computation will become inefficient.

seriesPlain :: Basis -> Series -> T Source

splitAtPadZero :: Int -> Mantissa -> (Mantissa, Mantissa) Source

Like `splitAt`

,
but it pads with zeros if the list is too short.
This way it preserves
` length (fst (splitAtPadZero n xs)) == n `

splitAtMatchPadZero :: [()] -> Mantissa -> (Mantissa, Mantissa) Source

truncSeriesSummands :: Series -> Series Source

help showing series summands

scaleSimple :: Basis -> T -> T Source

divide :: Basis -> T -> T -> T Source

Undefined if the divisor is zero - of course. Because it is impossible to assert that a real is zero, the routine will not throw an error in general.

ToDo: Rigorously derive the minimal required magnitude of the leading divisor digit.

reciprocal :: Basis -> T -> T Source

divIntMant :: Basis -> Digit -> Mantissa -> Mantissa Source

Fast division for small integral divisors, which occur for instance in summands of power series.

# algebraic functions

sqrtDriver :: Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T Source

sqrtFP :: Basis -> FixedPoint -> Mantissa Source

Square root.

We need a leading digit of type Integer,
because we have to collect up to 4 digits.
This presentation can also be considered as `FixedPoint`

.

ToDo: Rigorously derive the minimal required magnitude of the leading digit of the root.

Mathematically the `n`

th digit of the square root
depends roughly only on the first `n`

digits of the input.
This is because `sqrt (1+eps) `

.
However this implementation requires `equalApprox`

1 + eps/2`2*n`

input digits
for emitting `n`

digits.
This is due to the repeated use of `compressMant`

.
It would suffice to fully compress only every `basis`

th iteration (digit)
and compress only the second leading digit in each iteration.

Can the involved operations be made lazy enough to solve
`y = (x+frac)^2`

by
`frac = (y-x^2-frac^2) / (2*x)`

?

sqrtNewton :: Basis -> T -> T Source

sqrtFPNewton :: Basis -> FixedPoint -> Mantissa Source

Newton iteration doubles the number of correct digits in every step. Thus we process the data in chunks of sizes of powers of two. This way we get fastest computation possible with Newton but also more dependencies on input than necessary. The question arises whether this implementation still fits the needs of computational reals. The input is requested as larger and larger chunks, and the input itself might be computed this way, e.g. a repeated square root. Requesting one digit too much, requires the double amount of work for the input computation, which in turn multiplies time consumption by a factor of four, and so on.

Optimal fast implementation of one routine does not preserve fast computation of composed computations.

The routine assumes, that the integer parts is at least `b^2.`

lazyInits :: [a] -> [[a]] Source

List.inits is defined by
`inits = foldr (x ys -> [] : map (x:) ys) [[]]`

This is too strict for our application.

Prelude> List.inits (0:1:2:undefined) [[],[0],[0,1]*** Exception: Prelude.undefined

The following routine is more lazy than `inits`

and even lazier than `inits`

from `utility-ht`

package,
but it is restricted to infinite lists.
This degree of laziness is needed for `sqrtFP`

.

Prelude> lazyInits (0:1:2:undefined) [[],[0],[0,1],[0,1,2],[0,1,2,*** Exception: Prelude.undefined

# transcendent functions

## exponential functions

expSeriesLazy :: Basis -> T -> Series Source

expSmallLazy :: Basis -> T -> T Source

powerSeries :: Basis -> Rational -> T -> Series Source

Residue estimates will only hold for exponents with absolute value below one.

The computation is based on `Int`

,
thus the denominator should not be too big.
(Say, at most 1000 for 1000000 digits.)

It is not optimal to split the power into pure root and pure power (that means, with integer exponents). The root series can nicely handle all exponents, but for exponents above 1 the series summands rises at the beginning and thus make the residue estimate complicated. For powers with integer exponents the root series turns into the binomial formula, which is just a complicated way to compute a power which can also be determined by simple multiplication.

cosSinFourth :: Basis -> T -> (T, T) Source

Like `cosSinSmall`

but converges faster.
It calls `cosSinSmall`

with reduced arguments
using the trigonometric identities
```
cos (4*x) = 8 * cos x ^ 2 * (cos x ^ 2 - 1) + 1
sin (4*x) = 4 * sin x * cos x * (1 - 2 * sin x ^ 2)
```

Note that the faster convergence is hidden by the overhead.

The same could be achieved with a fourth power of a complex number.

## logarithmic functions

lnNewton :: Basis -> T -> T Source

x' = x - (exp x - y) / exp x = x + (y * exp (-x) - 1)

First, the dependencies on low-significant places are currently
much more than mathematically necessary.
Check
```
*Number.Positional> expSmall 1000 (-1,100 : replicate 16 0 ++ [undefined])
(0,[1,105,171,-82,76*** Exception: Prelude.undefined
```

Every multiplication cut off two trailing digits.
```
*Number.Positional> nest 8 (mul 1000 (-1,repeat 1)) (-1,100 : replicate 16 0 ++ [undefined])
(-9,[101,*** Exception: Prelude.undefined
```

Possibly the dependencies of expSmall
could be resolved by not computing `mul`

immediately
but computing `mul`

series which are merged and subsequently added.
But this would lead to an explosion of series.

Second, even if the dependencies of all atomic operations
are reduced to a minimum,
the mathematical dependencies of the whole iteration function
are less than the sums of the parts.
Lets demonstrate this with the square root iteration.
It is
```
(1.4140 + 2
```

That is, the digits *1.4140) * 2 == 1.414213578500707
(1.4149 + 2*1.4149) * 2 == 1.4142137288854335
`213`

do not depend mathematically on `x`

of `1.414x`

,
but their computation depends.
Maybe there is a glorious trick to reduce the computational dependencies
to the mathematical ones.

angle :: Basis -> (T, T) -> T Source

This is an inverse of `cosSin`

,
also known as `atan2`

with flipped arguments.
It's very slow because of the computation of sinus and cosinus.
However, because it uses the `atan2`

implementation as estimator,
the final application of arctan series should converge rapidly.

It could be certainly accelerated by not using cosSin and its fiddling with pi. Instead we could analyse quadrants before calling atan2, then calling cosSinSmall immediately.

arctanSeries :: Basis -> T -> Series Source

Arcus tangens of arguments with absolute value less than `1 / sqrt 3`

.

arctanSmall :: Basis -> T -> T Source

arctanStem :: Basis -> Digit -> T Source

Efficient computation of Arcus tangens of an argument of the form `1/n`

.

arctan :: Basis -> T -> T Source

This implementation gets the first decimal place for free
by calling the arcus tangens implementation for `Double`

s.

arctanClassic :: Basis -> T -> T Source

A classic implementation without ''cheating'' with floating point implementations.

For `x < 1 / sqrt 3`

(`1 / sqrt 3 == tan (pi/6)`

)
use `arctan`

power series.
(`sqrt 3`

is approximately `19/11`

.)

For `x > sqrt 3`

use
`arctan x = pi/2 - arctan (1/x)`

For other `x`

use

`arctan x = pi/4 - 0.5*arctan ((1-x^2)/2*x)`

(which follows from
`arctan x + arctan y == arctan ((x+y) / (1-x*y))`

which in turn follows from complex multiplication
`(1:+x)*(1:+y) == ((1-x*y):+(x+y))`

If `x`

is close to `sqrt 3`

or `1 / sqrt 3`

the computation is quite inefficient.

# constants

## elementary

## transcendental

recipEConst :: Basis -> T Source

# auxilary functions

sliceVertPair :: [a] -> [(a, a)] Source