Number.Positional
 Stability provisional Maintainer numericprelude@henning-thielemann.de
 Contents types basic helpers conversions integer rational fixed point floating point text basis comparison arithmetic algebraic functions transcendent functions exponential functions logarithmic functions constants elementary transcendental auxilary functions
Description
Exact Real Arithmetic - Computable reals. Inspired by ''The most unreliable technique for computing pi.'' See also http://www.haskell.org/haskellwiki/Exact_real_arithmetic .
Synopsis
 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 align :: Basis -> Exponent -> T -> T alignMant :: Basis -> Exponent -> T -> Mantissa absolute :: T -> T absMant :: Mantissa -> Mantissa fromLaurent :: T Int -> T toLaurent :: T -> T Int liftLaurent2 :: (T Int -> T Int -> T Int) -> T -> T -> T liftLaurentMany :: ([T Int] -> T Int) -> [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 -> Int -> Mantissa -> Mantissa divIntMantInf :: Basis -> Int -> 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 -> Int -> 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 T = (Exponent, Mantissa) Source
 type FixedPoint = (Integer, Mantissa) Source
 type Mantissa = [Digit] Source
 type Digit = Int Source
 type Exponent = Int Source
 type Basis = Int Source
basic helpers
 moveToZero :: Basis -> Digit -> (Digit, Digit) Source
 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.
 compress :: Basis -> T -> T Source
May prepend a digit.
 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
 trim :: T -> T Source
Eliminate leading zero digits. This will fail for zero.
 trimUntil :: Exponent -> T -> T Source
Trim until a minimum exponent is reached. Safe for zeros.
 trimOnce :: 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
 mantissaFromFixedInt :: Basis -> Exponent -> Integer -> 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
 toDouble :: Basis -> T -> Double Source
 fromDouble :: Basis -> Double -> T Source
cf. Numeric.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
 liftDoubleApprox :: Basis -> (Double -> Double) -> T -> T Source
 liftDoubleRough :: Basis -> (Double -> Double) -> T -> T Source
text
 showDec :: Exponent -> T -> String Source
Show a number with respect to basis 10^e.
 showHex :: Exponent -> T -> String Source
 showBin :: Exponent -> T -> String Source
 showBasis :: Basis -> Exponent -> T -> String Source
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).
 fromBasisMant :: Basis -> Basis -> Mantissa -> Mantissa Source
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.
 lessApprox :: Basis -> Exponent -> T -> T -> Bool Source
 trunc :: Exponent -> T -> T Source
 equalApprox :: Basis -> Exponent -> T -> T -> Bool Source
 align :: Basis -> Exponent -> T -> T 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.

 absolute :: T -> T Source
 absMant :: Mantissa -> Mantissa Source
arithmetic
 fromLaurent :: T Int -> T Source
 toLaurent :: T -> T Int Source
 liftLaurent2 :: (T Int -> T Int -> T Int) -> T -> T -> T Source
 liftLaurentMany :: ([T Int] -> T Int) -> [T] -> T Source
 add :: Basis -> T -> T -> T Source
Add two numbers but do not eliminate leading zeros.
 sub :: Basis -> T -> T -> T Source
 neg :: Basis -> T -> 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.
 type Series = [(Exponent, T)] Source
 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
 scale :: Basis -> Digit -> T -> T Source
 scaleSimple :: Basis -> T -> T Source
 scaleMant :: Basis -> Digit -> Mantissa -> Mantissa Source
 mulSeries :: Basis -> T -> T -> Series Source
 mul :: Basis -> T -> T -> T Source
For obtaining n result digits it is mathematically sufficient to know the first (n+1) digits of the operands. However this implementation needs (n+2) digits, because of calls to compress in both scale and series. We should fix that.
 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.

 divMant :: Basis -> Mantissa -> Mantissa -> Mantissa Source
 divMantSlow :: Basis -> Mantissa -> Mantissa -> Mantissa Source
 reciprocal :: Basis -> T -> T Source
 divIntMant :: Basis -> Int -> Mantissa -> Mantissa Source
Fast division for small integral divisors, which occur for instance in summands of power series.
 divIntMantInf :: Basis -> Int -> Mantissa -> Mantissa Source
 divInt :: Basis -> Digit -> T -> T Source
algebraic functions
 sqrt :: Basis -> T -> T Source
 sqrtDriver :: Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T Source
 sqrtMant :: Basis -> Mantissa -> Mantissa 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 nth digit of the square root depends roughly only on the first n digits of the input. This is because sqrt (1+eps) equalApprox 1 + eps/2. However this implementation requires 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 basisth 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,1]*** Exception: Prelude.undefined
```

The following routine is more lazy than inits and even lazier than Data.List.HT.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,1],[0,1,2],[0,1,2,*** Exception: Prelude.undefined
```
transcendent functions
exponential functions
 expSeries :: Basis -> T -> Series Source
 expSmall :: Basis -> T -> T Source
Absolute value of argument should be below 1.
 expSeriesLazy :: Basis -> T -> Series Source
 expSmallLazy :: Basis -> T -> T Source
 exp :: Basis -> T -> T Source
 intPower :: Basis -> Integer -> T -> T -> T -> T Source
 cardPower :: Basis -> Integer -> T -> 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.

 powerSmall :: Basis -> Rational -> T -> T Source
 power :: Basis -> Rational -> T -> T Source
 root :: Basis -> Integer -> T -> T Source
 cosSinhSmall :: Basis -> T -> (T, T) Source
Absolute value of argument should be below 1.
 cosSinSmall :: Basis -> T -> (T, T) Source
Absolute value of argument should be below 1.
 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.

 cosSin :: Basis -> T -> (T, T) Source
 tan :: Basis -> T -> T Source
 cot :: Basis -> T -> T Source
logarithmic functions
 lnSeries :: Basis -> T -> Series Source
 lnSmall :: Basis -> T -> T Source
 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 + 21.4140) 2 == 1.414213578500707 (1.4149 + 21.4149) 2 == 1.4142137288854335 That is, the digits 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.

 lnNewton' :: Basis -> T -> T Source
 ln :: Basis -> T -> T Source
 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 -> Int -> 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 Doubles.
 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
 zero :: T Source
 one :: T Source
 minusOne :: T Source
transcendental
 eConst :: Basis -> T Source
 recipEConst :: Basis -> T Source
 piConst :: Basis -> T Source
auxilary functions
 sliceVertPair :: [a] -> [(a, a)] Source
Produced by Haddock version 2.4.2