{-# LANGUAGE RebindableSyntax #-}
{- |
Exact Real Arithmetic - Computable reals.
Inspired by ''The most unreliable technique for computing pi.''
See also <http://www.haskell.org/haskellwiki/Exact_real_arithmetic> .
-}
module Number.Positional where

import qualified MathObj.LaurentPolynomial as LPoly
import qualified MathObj.Polynomial.Core   as Poly

import qualified Algebra.IntegralDomain as Integral
import qualified Algebra.Ring           as Ring
import qualified Algebra.ToInteger      as ToInteger

import qualified Prelude as P98
import qualified NumericPrelude.Numeric as NP

import NumericPrelude.Base
import NumericPrelude.Numeric hiding (sqrt, tan, one, zero, )

import qualified Data.List as List
import Data.Char (intToDigit)

import qualified Data.List.Match as Match
import Data.Function.HT (powerAssociative, nest, )
import Data.Tuple.HT (swap, )
import Data.Maybe.HT (toMaybe, )
import Data.Bool.HT (select, if', )
import NumericPrelude.List (mapLast, )
import Data.List.HT
          (sliceVertical, mapAdjacent,
           padLeft, padRight, )


{-
FIXME:

defltBase = 10
defltExp = 4

(sqrt 0.5) -- wrong result, probably due to undetected overflows
-}

{- * types -}

type T = (Exponent, Mantissa)
type FixedPoint = (Integer, Mantissa)
type Mantissa = [Digit]
type Digit    = Int
type Exponent = Int
type Basis    = Int


{- * basic helpers -}

moveToZero :: Basis -> Digit -> (Digit,Digit)
moveToZero :: Basis -> Basis -> (Basis, Basis)
moveToZero Basis
b Basis
n =
   let b2 :: Basis
b2 = Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div Basis
b Basis
2
       (Basis
q,Basis
r) = Basis -> Basis -> (Basis, Basis)
forall a. C a => a -> a -> (a, a)
divMod (Basis
nBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
b2) Basis
b
   in  (Basis
q,Basis
rBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
b2)

checkPosDigit :: Basis -> Digit -> Digit
checkPosDigit :: Basis -> Basis -> Basis
checkPosDigit Basis
b Basis
d =
   if Basis
dBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
0 Bool -> Bool -> Bool
&& Basis
dBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<Basis
b
     then Basis
d
     else [Char] -> Basis
forall a. HasCallStack => [Char] -> a
error ([Char]
"digit " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Basis -> [Char]
forall a. Show a => a -> [Char]
show Basis
d [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" out of range [0," [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Basis -> [Char]
forall a. Show a => a -> [Char]
show Basis
b [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
")")

checkDigit :: Basis -> Digit -> Digit
checkDigit :: Basis -> Basis -> Basis
checkDigit Basis
b Basis
d =
   if Basis -> Basis
forall a. C a => a -> a
abs Basis
d Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
< Basis
b
     then Basis
d
     else [Char] -> Basis
forall a. HasCallStack => [Char] -> a
error ([Char]
"digit " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Basis -> [Char]
forall a. Show a => a -> [Char]
show Basis
d [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" out of range ("
                   [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Basis -> [Char]
forall a. Show a => a -> [Char]
show (-Basis
b) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"," [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Basis -> [Char]
forall a. Show a => a -> [Char]
show Basis
b [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
")")

{- |
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!)
-}
nonNegative :: Basis -> T -> T
nonNegative :: Basis -> T -> T
nonNegative Basis
b T
x =
   let (Basis
xe,Mantissa
xm) = Basis -> T -> T
compress Basis
b T
x
   in  (Basis
xe, Basis -> Mantissa -> Mantissa
nonNegativeMant Basis
b Mantissa
xm)

{- |
Requires, that no digit is @(basis-1)@ or @(1-basis)@.
The leading digit might be negative and might be @-basis@ or @basis@.
-}
nonNegativeMant :: Basis -> Mantissa -> Mantissa
nonNegativeMant :: Basis -> Mantissa -> Mantissa
nonNegativeMant Basis
b =
   let recurse :: Mantissa -> Mantissa
recurse (Basis
x0:Basis
x1:Mantissa
xs) =
          Mantissa -> [(Bool, Mantissa)] -> Mantissa
forall a. a -> [(Bool, a)] -> a
select (Basis -> Mantissa -> Mantissa
replaceZeroChain Basis
x0 (Basis
x1Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs))
             [(Basis
x1 Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=  Basis
1,  Basis
x0    Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Mantissa
recurse (Basis
x1Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs)),
              (Basis
x1 Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<= -Basis
1, (Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Mantissa
recurse ((Basis
x1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
b)Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs))]
       recurse Mantissa
xs = Mantissa
xs

       replaceZeroChain :: Basis -> Mantissa -> Mantissa
replaceZeroChain Basis
x Mantissa
xs =
          let (Mantissa
xZeros,Mantissa
xRem) = (Basis -> Bool) -> Mantissa -> (Mantissa, Mantissa)
forall a. (a -> Bool) -> [a] -> ([a], [a])
span (Basis
0Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
==) Mantissa
xs
          in  case Mantissa
xRem of
                [] -> (Basis
xBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs)  -- keep trailing zeros, because they show precision in 'show' functions
                (Basis
y:Mantissa
ys) ->
                  if Basis
yBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
0  -- equivalent to y>0
                    then Basis
x     Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Basis -> Mantissa
forall a b. [a] -> b -> [b]
Match.replicate Mantissa
xZeros Basis
0     Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Mantissa -> Mantissa
recurse Mantissa
xRem
                    else (Basis
xBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Basis -> Mantissa
forall a b. [a] -> b -> [b]
Match.replicate Mantissa
xZeros (Basis
bBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Mantissa -> Mantissa
recurse ((Basis
yBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
b) Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
ys)

   in  Mantissa -> Mantissa
recurse


{- |
May prepend a digit.
-}
compress :: Basis -> T -> T
compress :: Basis -> T -> T
compress Basis
_ x :: T
x@(Basis
_, []) = T
x
compress Basis
b (Basis
xe, Mantissa
xm) =
   let (Basis
hi:Mantissa
his,Mantissa
los) = [(Basis, Basis)] -> (Mantissa, Mantissa)
forall a b. [(a, b)] -> ([a], [b])
unzip ((Basis -> (Basis, Basis)) -> Mantissa -> [(Basis, Basis)]
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Basis -> (Basis, Basis)
moveToZero Basis
b) Mantissa
xm)
   in  Basis -> T -> T
prependDigit Basis
hi (Basis
xe, Mantissa -> Mantissa -> Mantissa
forall a. C a => [a] -> [a] -> [a]
Poly.add Mantissa
his Mantissa
los)

{- |
Compress first digit.
May prepend a digit.
-}
compressFirst :: Basis -> T -> T
compressFirst :: Basis -> T -> T
compressFirst Basis
_ x :: T
x@(Basis
_, []) = T
x
compressFirst Basis
b (Basis
xe, Basis
x:Mantissa
xs) =
   let (Basis
hi,Basis
lo) = Basis -> Basis -> (Basis, Basis)
moveToZero Basis
b Basis
x
   in  Basis -> T -> T
prependDigit Basis
hi (Basis
xe, Basis
loBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs)

{- |
Does not prepend a digit.
-}
compressMant :: Basis -> Mantissa -> Mantissa
compressMant :: Basis -> Mantissa -> Mantissa
compressMant Basis
_ [] = []
compressMant Basis
b (Basis
x:Mantissa
xs) =
   let (Mantissa
his,Mantissa
los) = [(Basis, Basis)] -> (Mantissa, Mantissa)
forall a b. [(a, b)] -> ([a], [b])
unzip ((Basis -> (Basis, Basis)) -> Mantissa -> [(Basis, Basis)]
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Basis -> (Basis, Basis)
moveToZero Basis
b) Mantissa
xs)
   in  Mantissa -> Mantissa -> Mantissa
forall a. C a => [a] -> [a] -> [a]
Poly.add Mantissa
his (Basis
xBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
los)

{- |
Compress second digit.
Sometimes this is enough to keep the digits in the admissible range.
Does not prepend a digit.
-}
compressSecondMant :: Basis -> Mantissa -> Mantissa
compressSecondMant :: Basis -> Mantissa -> Mantissa
compressSecondMant Basis
b (Basis
x0:Basis
x1:Mantissa
xs) =
   let (Basis
hi,Basis
lo) = Basis -> Basis -> (Basis, Basis)
moveToZero Basis
b Basis
x1
   in  Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
hi Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Basis
lo Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
xs
compressSecondMant Basis
_ Mantissa
xm = Mantissa
xm

prependDigit :: Basis -> T -> T
prependDigit :: Basis -> T -> T
prependDigit Basis
0 T
x = T
x
prependDigit Basis
x (Basis
xe, Mantissa
xs) = (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1, Basis
xBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs)

{- |
Eliminate leading zero digits.
This will fail for zero.
-}
trim :: T -> T
trim :: T -> T
trim (Basis
xe,Mantissa
xm) =
   let (Mantissa
xZero, Mantissa
xNonZero) = (Basis -> Bool) -> Mantissa -> (Mantissa, Mantissa)
forall a. (a -> Bool) -> [a] -> ([a], [a])
span (Basis
0 Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
==) Mantissa
xm
   in  (Basis
xe Basis -> Basis -> Basis
forall a. C a => a -> a -> a
- Mantissa -> Basis
forall (t :: * -> *) a. Foldable t => t a -> Basis
length Mantissa
xZero, Mantissa
xNonZero)

{- |
Trim until a minimum exponent is reached.
Safe for zeros.
-}
trimUntil :: Exponent -> T -> T
trimUntil :: Basis -> T -> T
trimUntil Basis
e =
   (T -> Bool) -> (T -> T) -> T -> T
forall a. (a -> Bool) -> (a -> a) -> a -> a
until (\(Basis
xe,Mantissa
xm) -> Basis
xeBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<=Basis
e Bool -> Bool -> Bool
||
              Bool -> Bool
not (Mantissa -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null Mantissa
xm Bool -> Bool -> Bool
|| Mantissa -> Basis
forall a. [a] -> a
head Mantissa
xm Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
== Basis
0)) T -> T
trimOnce

trimOnce :: T -> T
trimOnce :: T -> T
trimOnce (Basis
xe,[])   = (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1,[])
trimOnce (Basis
xe,Basis
0:Mantissa
xm) = (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1,Mantissa
xm)
trimOnce T
x = T
x

{- |
Accept a high leading digit for the sake of a reduced exponent.
This eliminates one leading digit.
Like 'pumpFirst' but with exponent management.
-}
decreaseExp :: Basis -> T -> T
decreaseExp :: Basis -> T -> T
decreaseExp Basis
b (Basis
xe,Mantissa
xm) =
   (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1, Basis -> Mantissa -> Mantissa
pumpFirst Basis
b Mantissa
xm)

{- |
Merge leading and second digit.
This is somehow an inverse of 'compressMant'.
-}
pumpFirst :: Basis -> Mantissa -> Mantissa
pumpFirst :: Basis -> Mantissa -> Mantissa
pumpFirst Basis
b Mantissa
xm =
   case Mantissa
xm of
     (Basis
x0:Basis
x1:Mantissa
xs) -> Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
bBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
x1Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs
     (Basis
x0:[])    -> Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
bBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:[]
     []         -> []

decreaseExpFP :: Basis -> (Exponent, FixedPoint) ->
                          (Exponent, FixedPoint)
decreaseExpFP :: Basis -> (Basis, FixedPoint) -> (Basis, FixedPoint)
decreaseExpFP Basis
b (Basis
xe,FixedPoint
xm) =
   (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1, Basis -> FixedPoint -> FixedPoint
pumpFirstFP Basis
b FixedPoint
xm)

pumpFirstFP :: Basis -> FixedPoint -> FixedPoint
pumpFirstFP :: Basis -> FixedPoint -> FixedPoint
pumpFirstFP Basis
b (Integer
x,Mantissa
xm) =
   let xb :: Integer
xb = Integer
x Integer -> Integer -> Integer
forall a. C a => a -> a -> a
* Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
b
   in  case Mantissa
xm of
         (Basis
x0:Mantissa
xs) -> (Integer
xb Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
x0, Mantissa
xs)
         []      -> (Integer
xb, [])

{- |
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.
-}
negativeExp :: Basis -> T -> T
negativeExp :: Basis -> T -> T
negativeExp Basis
b T
x =
   let tx :: T
tx  = Basis -> T -> T
trimUntil (-Basis
10) T
x
   in  if T -> Basis
forall a b. (a, b) -> a
fst T
tx Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
0 then Basis -> T -> T
decreaseExp Basis
b T
tx else T
tx


{- * conversions -}

{- ** integer -}

fromBaseCardinal :: Basis -> Integer -> T
fromBaseCardinal :: Basis -> Integer -> T
fromBaseCardinal Basis
b Integer
n =
   let mant :: Mantissa
mant = Basis -> Integer -> Mantissa
forall a. C a => Basis -> a -> Mantissa
mantissaFromCard Basis
b Integer
n
   in  (Mantissa -> Basis
forall (t :: * -> *) a. Foldable t => t a -> Basis
length Mantissa
mant Basis -> Basis -> Basis
forall a. C a => a -> a -> a
- Basis
1, Mantissa
mant)

fromBaseInteger :: Basis -> Integer -> T
fromBaseInteger :: Basis -> Integer -> T
fromBaseInteger Basis
b Integer
n =
   if Integer
nInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<Integer
0
     then Basis -> T -> T
neg Basis
b (Basis -> Integer -> T
fromBaseCardinal Basis
b (Integer -> Integer
forall a. C a => a -> a
negate Integer
n))
     else Basis -> Integer -> T
fromBaseCardinal Basis
b Integer
n

mantissaToNum :: Ring.C a => Basis -> Mantissa -> a
mantissaToNum :: Basis -> Mantissa -> a
mantissaToNum Basis
bInt =
   let b :: a
b = Basis -> a
forall a b. (C a, C b) => a -> b
fromIntegral Basis
bInt
   in  (a -> Basis -> a) -> a -> Mantissa -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\a
x Basis
d -> a
xa -> a -> a
forall a. C a => a -> a -> a
*a
b a -> a -> a
forall a. C a => a -> a -> a
+ Basis -> a
forall a b. (C a, C b) => a -> b
fromIntegral Basis
d) a
0

mantissaFromCard :: (ToInteger.C a) => Basis -> a -> Mantissa
mantissaFromCard :: Basis -> a -> Mantissa
mantissaFromCard Basis
bInt a
n =
   let b :: a
b = Basis -> a
forall a b. (C a, C b) => a -> b
NP.fromIntegral Basis
bInt
   in  Mantissa -> Mantissa
forall a. [a] -> [a]
reverse ((a -> Basis) -> [a] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map a -> Basis
forall a b. (C a, C b) => a -> b
NP.fromIntegral
          ([a] -> a -> [a]
forall a. (C a, C a) => [a] -> a -> [a]
Integral.decomposeVarPositional (a -> [a]
forall a. a -> [a]
repeat a
b) a
n))

mantissaFromInt :: (ToInteger.C a) => Basis -> a -> Mantissa
mantissaFromInt :: Basis -> a -> Mantissa
mantissaFromInt Basis
b a
n =
   if a
na -> a -> Bool
forall a. Ord a => a -> a -> Bool
<a
0
     then Mantissa -> Mantissa
forall a. C a => a -> a
negate (Basis -> a -> Mantissa
forall a. C a => Basis -> a -> Mantissa
mantissaFromCard Basis
b (a -> a
forall a. C a => a -> a
negate a
n))
     else Basis -> a -> Mantissa
forall a. C a => Basis -> a -> Mantissa
mantissaFromCard Basis
b a
n

mantissaFromFixedInt :: Basis -> Exponent -> Integer -> Mantissa
mantissaFromFixedInt :: Basis -> Basis -> Integer -> Mantissa
mantissaFromFixedInt Basis
bInt Basis
e Integer
n =
   let b :: Integer
b = Basis -> Integer
forall a b. (C a, C b) => a -> b
NP.fromIntegral Basis
bInt
   in  (Integer -> Basis) -> [Integer] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Basis
forall a b. (C a, C b) => a -> b
NP.fromIntegral ((Integer -> [Integer] -> [Integer])
-> (Integer, [Integer]) -> [Integer]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (:) ((Integer -> () -> (Integer, Integer))
-> Integer -> [()] -> (Integer, [Integer])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumR
          (\Integer
x () -> Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod Integer
x Integer
b)
          Integer
n (Basis -> () -> [()]
forall a. Basis -> a -> [a]
replicate (Basis -> Basis
forall a. Enum a => a -> a
pred Basis
e) ())))


{- ** rational -}

fromBaseRational :: Basis -> Rational -> T
fromBaseRational :: Basis -> Rational -> T
fromBaseRational Basis
bInt Rational
x =
   let b :: Integer
b = Basis -> Integer
forall a b. (C a, C b) => a -> b
NP.fromIntegral Basis
bInt
       xDen :: Integer
xDen = Rational -> Integer
forall a. T a -> a
denominator Rational
x
       (Integer
xInt,Integer
xNum) = Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod (Rational -> Integer
forall a. T a -> a
numerator Rational
x) Integer
xDen
       (Basis
xe,Mantissa
xm) = Basis -> Integer -> T
fromBaseInteger Basis
bInt Integer
xInt
       xFrac :: [Integer]
xFrac = (Integer -> Maybe (Integer, Integer)) -> Integer -> [Integer]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
List.unfoldr
                 (\Integer
num -> Bool -> (Integer, Integer) -> Maybe (Integer, Integer)
forall a. Bool -> a -> Maybe a
toMaybe (Integer
numInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=Integer
0) (Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod (Integer
numInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
b) Integer
xDen)) Integer
xNum
   in  (Basis
xe, Mantissa
xm Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ (Integer -> Basis) -> [Integer] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Basis
forall a. C a => Integer -> a
NP.fromInteger [Integer]
xFrac)

{- ** fixed point -}

{- |
Split into integer and fractional part.
-}
toFixedPoint :: Basis -> T -> FixedPoint
toFixedPoint :: Basis -> T -> FixedPoint
toFixedPoint Basis
b (Basis
xe,Mantissa
xm) =
   if Basis
xeBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
0
     then let (Mantissa
x0,Mantissa
x1) = Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1) Mantissa
xm
          in  (Basis -> Mantissa -> Integer
forall a. C a => Basis -> Mantissa -> a
mantissaToNum Basis
b Mantissa
x0, Mantissa
x1)
     else (Integer
0, Basis -> Basis -> Mantissa
forall a. Basis -> a -> [a]
replicate (- Basis -> Basis
forall a. Enum a => a -> a
succ Basis
xe) Basis
0 Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Mantissa
xm)

fromFixedPoint :: Basis -> FixedPoint -> T
fromFixedPoint :: Basis -> FixedPoint -> T
fromFixedPoint Basis
b (Integer
xInt,Mantissa
xFrac) =
   let (Basis
xe,Mantissa
xm) = Basis -> Integer -> T
fromBaseInteger Basis
b Integer
xInt
   in  (Basis
xe, Mantissa
xmMantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++Mantissa
xFrac)


{- ** floating point -}

toDouble :: Basis -> T -> Double
toDouble :: Basis -> T -> Double
toDouble Basis
b (Basis
xe,Mantissa
xm) =
   let txm :: Mantissa
txm = Basis -> Mantissa -> Mantissa
forall a. Basis -> [a] -> [a]
take (Basis -> Basis
mantLengthDouble Basis
b) Mantissa
xm
       bf :: Double
bf  = Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
b
       br :: Double
br  = Double -> Double
forall a. C a => a -> a
recip Double
bf
   in  Basis -> Double -> Double
forall a b. (C a, C b) => b -> a -> a
fieldPower Basis
xe Double
bf Double -> Double -> Double
forall a. C a => a -> a -> a
* (Basis -> Double -> Double) -> Double -> Mantissa -> Double
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Basis
xi Double
y -> Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
xi Double -> Double -> Double
forall a. C a => a -> a -> a
+ Double
yDouble -> Double -> Double
forall a. C a => a -> a -> a
*Double
br) Double
0 Mantissa
txm

{- |
cf. 'Numeric.floatToDigits'
-}
fromDouble :: Basis -> Double -> T
fromDouble :: Basis -> Double -> T
fromDouble Basis
b Double
x =
   let (Integer
n,Double
frac) = Double -> (Integer, Double)
forall a b. (C a, C b) => a -> (b, a)
splitFraction Double
x
       (Integer
mant,Basis
e) = Double -> (Integer, Basis)
forall a. RealFloat a => a -> (Integer, Basis)
P98.decodeFloat Double
frac
       fracR :: Mantissa
fracR    = Basis -> Basis -> T -> Mantissa
alignMant Basis
b (-Basis
1)
                     (Basis -> Rational -> T
fromBaseRational Basis
b (Integer
mant Integer -> Integer -> Rational
forall a. C a => a -> a -> T a
% Basis -> Integer -> Integer
forall a b. (C a, C b) => b -> a -> a
ringPower (-Basis
e) Integer
2))
   in  Basis -> FixedPoint -> T
fromFixedPoint Basis
b (Integer
n, Mantissa
fracR)

{- |
Only return as much digits as are contained in Double.
This will speedup further computations.
-}
fromDoubleApprox :: Basis -> Double -> T
fromDoubleApprox :: Basis -> Double -> T
fromDoubleApprox Basis
b Double
x =
   let (Basis
xe,Mantissa
xm) = Basis -> Double -> T
fromDouble Basis
b Double
x
   in  (Basis
xe, Basis -> Mantissa -> Mantissa
forall a. Basis -> [a] -> [a]
take (Basis -> Basis
mantLengthDouble Basis
b) Mantissa
xm)

fromDoubleRough :: Basis -> Double -> T
fromDoubleRough :: Basis -> Double -> T
fromDoubleRough Basis
b Double
x =
   let (Basis
xe,Mantissa
xm) = Basis -> Double -> T
fromDouble Basis
b Double
x
   in  (Basis
xe, Basis -> Mantissa -> Mantissa
forall a. Basis -> [a] -> [a]
take Basis
2 Mantissa
xm)

mantLengthDouble :: Basis -> Exponent
mantLengthDouble :: Basis -> Basis
mantLengthDouble Basis
b =
   let fi :: Basis -> Double
fi = Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral :: Int -> Double
       x :: Double
x  = Double
forall a. HasCallStack => a
undefined :: Double
   in  Double -> Basis
forall a b. (C a, C b) => a -> b
ceiling
          (Double -> Double -> Double
forall a. C a => a -> a -> a
logBase (Basis -> Double
fi Basis
b) (Integer -> Double
forall a. C a => Integer -> a
fromInteger (Double -> Integer
forall a. RealFloat a => a -> Integer
P98.floatRadix Double
x)) Double -> Double -> Double
forall a. C a => a -> a -> a
*
             Basis -> Double
fi (Double -> Basis
forall a. RealFloat a => a -> Basis
P98.floatDigits Double
x))

liftDoubleApprox :: Basis -> (Double -> Double) -> T -> T
liftDoubleApprox :: Basis -> (Double -> Double) -> T -> T
liftDoubleApprox Basis
b Double -> Double
f = Basis -> Double -> T
fromDoubleApprox Basis
b (Double -> T) -> (T -> Double) -> T -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
f (Double -> Double) -> (T -> Double) -> T -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Basis -> T -> Double
toDouble Basis
b

liftDoubleRough :: Basis -> (Double -> Double) -> T -> T
liftDoubleRough :: Basis -> (Double -> Double) -> T -> T
liftDoubleRough Basis
b Double -> Double
f = Basis -> Double -> T
fromDoubleRough Basis
b (Double -> T) -> (T -> Double) -> T -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
f (Double -> Double) -> (T -> Double) -> T -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Basis -> T -> Double
toDouble Basis
b


{- ** text -}

{- |
Show a number with respect to basis @10^e@.
-}
showDec :: Exponent -> T -> String
showDec :: Basis -> T -> [Char]
showDec = Basis -> Basis -> T -> [Char]
showBasis Basis
10

showHex :: Exponent -> T -> String
showHex :: Basis -> T -> [Char]
showHex = Basis -> Basis -> T -> [Char]
showBasis Basis
16

showBin :: Exponent -> T -> String
showBin :: Basis -> T -> [Char]
showBin = Basis -> Basis -> T -> [Char]
showBasis Basis
2

showBasis :: Basis -> Exponent -> T -> String
showBasis :: Basis -> Basis -> T -> [Char]
showBasis Basis
b Basis
e T
xBig =
   let x :: T
x = Basis -> Basis -> T -> T
rootBasis Basis
b Basis
e T
xBig
       ([Char]
sign,T
absX) =
          case Basis -> T -> T -> Ordering
cmp Basis
b T
x (T -> Basis
forall a b. (a, b) -> a
fst T
x,[]) of
             Ordering
LT -> ([Char]
"-", Basis -> T -> T
neg Basis
b T
x)
             Ordering
_  -> ([Char]
"", T
x)
       (Integer
int, Mantissa
frac) = Basis -> T -> FixedPoint
toFixedPoint Basis
b (Basis -> T -> T
nonNegative Basis
b T
absX)
       checkedFrac :: Mantissa
checkedFrac = (Basis -> Basis) -> Mantissa -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Basis -> Basis
checkPosDigit Basis
b) Mantissa
frac
       intStr :: [Char]
intStr =
          if Basis
bBasis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
==Basis
10
            then Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
int
            else (Basis -> Char) -> Mantissa -> [Char]
forall a b. (a -> b) -> [a] -> [b]
map Basis -> Char
intToDigit (Basis -> Integer -> Mantissa
forall a. C a => Basis -> a -> Mantissa
mantissaFromInt Basis
b Integer
int)
   in  [Char]
sign [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
intStr [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Char
'.' Char -> [Char] -> [Char]
forall a. a -> [a] -> [a]
: (Basis -> Char) -> Mantissa -> [Char]
forall a b. (a -> b) -> [a] -> [b]
map Basis -> Char
intToDigit Mantissa
checkedFrac


{- ** basis -}

{- |
Convert from a @b@ basis representation to a @b^e@ basis.

Works well with every exponent.
-}
powerBasis :: Basis -> Exponent -> T -> T
powerBasis :: Basis -> Basis -> T -> T
powerBasis Basis
b Basis
e (Basis
xe,Mantissa
xm) =
   let (Basis
ye,Basis
r)  = Basis -> Basis -> (Basis, Basis)
forall a. C a => a -> a -> (a, a)
divMod Basis
xe Basis
e
       (Mantissa
y0,Mantissa
y1) = Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero (Basis
rBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1) Mantissa
xm
       y1pad :: [Mantissa]
y1pad   = (Mantissa -> Mantissa) -> [Mantissa] -> [Mantissa]
forall a. (a -> a) -> [a] -> [a]
mapLast (Basis -> Basis -> Mantissa -> Mantissa
forall a. a -> Basis -> [a] -> [a]
padRight Basis
0 Basis
e) (Basis -> Mantissa -> [Mantissa]
forall a. Basis -> [a] -> [[a]]
sliceVertical Basis
e Mantissa
y1)
   in  (Basis
ye, (Mantissa -> Basis) -> [Mantissa] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Mantissa -> Basis
forall a. C a => Basis -> Mantissa -> a
mantissaToNum Basis
b) (Mantissa
y0 Mantissa -> [Mantissa] -> [Mantissa]
forall a. a -> [a] -> [a]
: [Mantissa]
y1pad))

{- |
Convert from a @b^e@ basis representation to a @b@ basis.

Works well with every exponent.
-}
rootBasis :: Basis -> Exponent -> T -> T
rootBasis :: Basis -> Basis -> T -> T
rootBasis Basis
b Basis
e (Basis
xe,Mantissa
xm) =
   let splitDigit :: a -> Mantissa
splitDigit a
d = Basis -> Basis -> Mantissa -> Mantissa
forall a. a -> Basis -> [a] -> [a]
padLeft Basis
0 Basis
e (Basis -> a -> Mantissa
forall a. C a => Basis -> a -> Mantissa
mantissaFromInt Basis
b a
d)
   in  Basis -> (T -> T) -> T -> T
forall a. Basis -> (a -> a) -> a -> a
nest (Basis
eBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) T -> T
trimOnce
            ((Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1)Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
eBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1, (Basis -> Mantissa) -> Mantissa -> Mantissa
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Basis -> Mantissa
forall a. C a => a -> Mantissa
splitDigit ((Basis -> Basis) -> Mantissa -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Basis -> Basis
checkDigit (Basis -> Basis -> Basis
forall a b. (C a, C b) => b -> a -> a
ringPower Basis
e Basis
b)) Mantissa
xm))

{- |
Convert between arbitrary bases.
This conversion is expensive (quadratic time).
-}
fromBasis :: Basis -> Basis -> T -> T
fromBasis :: Basis -> Basis -> T -> T
fromBasis Basis
bDst Basis
bSrc T
x =
   let (Integer
int,Mantissa
frac) = Basis -> T -> FixedPoint
toFixedPoint Basis
bSrc T
x
   in  Basis -> FixedPoint -> T
fromFixedPoint Basis
bDst (Integer
int, Basis -> Basis -> Mantissa -> Mantissa
fromBasisMant Basis
bDst Basis
bSrc Mantissa
frac)

fromBasisMant :: Basis -> Basis -> Mantissa -> Mantissa
fromBasisMant :: Basis -> Basis -> Mantissa -> Mantissa
fromBasisMant Basis
_    Basis
_    [] = []
fromBasisMant Basis
bDst Basis
bSrc Mantissa
xm =
   let {- We use a counter that alerts us,
          when the digits are grown too much by Poly.scale.
          Then it is time to do some carry/compression.
          'inc' is essentially the fractional number digits
          needed to represent the destination base in the source base.
          It is multiplied by 'unit' in order to allow integer computation. -}
       inc :: Basis
inc = Double -> Basis
forall a b. (C a, C b) => a -> b
ceiling
                (Double -> Double -> Double
forall a. C a => a -> a -> a
logBase (Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
bSrc) (Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
bDst)
                     Double -> Double -> Double
forall a. C a => a -> a -> a
* Double
forall a. C a => a
unit Double -> Double -> Double
forall a. C a => a -> a -> a
* Double
1.1 :: Double)
          -- Without the correction factor, invalid digits are emitted - why?
       unit :: Ring.C a => a
       unit :: a
unit = a
10000
       {-
       This would create finite representations
       in some cases (input is finite, and the result is finite)
       but will cause infinite loop otherwise.
       Rev.dropWhile (0==) . compressMant bDst
       -}
       cmpr :: (a, Mantissa) -> (a, Mantissa)
cmpr (a
mag,Mantissa
xs) = (a
mag a -> a -> a
forall a. C a => a -> a -> a
- a
forall a. C a => a
unit, Basis -> Mantissa -> Mantissa
compressMant Basis
bSrc Mantissa
xs)

       scl :: T -> Maybe (Basis, T)
scl (Basis
_,[]) = Maybe (Basis, T)
forall a. Maybe a
Nothing
       scl (Basis
mag,Mantissa
xs) =
          let (Basis
newMag,Basis
y:Mantissa
ys) =
                 (T -> Bool) -> (T -> T) -> T -> T
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<Basis
forall a. C a => a
unit) (Basis -> Bool) -> (T -> Basis) -> T -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T -> Basis
forall a b. (a, b) -> a
fst) T -> T
forall a. C a => (a, Mantissa) -> (a, Mantissa)
cmpr
                       (Basis
mag Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+ Basis
inc, Basis -> Mantissa -> Mantissa
forall a. C a => a -> [a] -> [a]
Poly.scale Basis
bDst Mantissa
xs)
              (Basis
d,Basis
y0) = Basis -> Basis -> (Basis, Basis)
moveToZero Basis
bSrc Basis
y
          in  (Basis, T) -> Maybe (Basis, T)
forall a. a -> Maybe a
Just (Basis
d, (Basis
newMag, Basis
y0Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
ys))

   in  (T -> Maybe (Basis, T)) -> T -> Mantissa
forall b a. (b -> Maybe (a, b)) -> b -> [a]
List.unfoldr T -> Maybe (Basis, T)
scl (Basis
0::Int,Mantissa
xm)


{- * comparison -}

{- |
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.
-}
cmp :: Basis -> T -> T -> Ordering
cmp :: Basis -> T -> T -> Ordering
cmp Basis
b T
x T
y =
   let (Basis
_,Mantissa
dm) = Basis -> T -> T -> T
sub Basis
b T
x T
y
       {- Only differences above 2 allow a safe decision,
          because 1(-9)(-9)(-9)(-9)... and (-1)9999...
          represent the same number, namely zero. -}
       recurse :: Mantissa -> Ordering
recurse [] = Ordering
EQ
       recurse (Basis
d:[]) = Basis -> Basis -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Basis
d Basis
0
       recurse (Basis
d0:Basis
d1:Mantissa
ds) =
          Ordering -> [(Bool, Ordering)] -> Ordering
forall a. a -> [(Bool, a)] -> a
select (Mantissa -> Ordering
recurse (Basis
d0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
bBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
d1 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
ds))
             [(Basis
d0 Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
< -Basis
2, Ordering
LT),
              (Basis
d0 Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>  Basis
2, Ordering
GT)]
   in  Mantissa -> Ordering
recurse Mantissa
dm

{-
Compare two numbers approximately.
This circumvents the infinite loop if both numbers are equal.
If @lessApprox bnd b x y@
then @x@ is definitely smaller than @y@.
The converse is not true.
You should use this one instead of 'cmp' for checking for bounds.
-}
lessApprox :: Basis -> Exponent -> T -> T -> Bool
lessApprox :: Basis -> Basis -> T -> T -> Bool
lessApprox Basis
b Basis
bnd T
x T
y =
   let tx :: T
tx = Basis -> T -> T
trunc Basis
bnd T
x
       ty :: T
ty = Basis -> T -> T
trunc Basis
bnd T
y
   in  Ordering
LT Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Basis -> T -> T -> Ordering
cmp Basis
b ((T Basis -> T Basis -> T Basis) -> T -> T -> T
liftLaurent2 T Basis -> T Basis -> T Basis
forall a. C a => T a -> T a -> T a
LPoly.add (Basis
bnd,[Basis
2]) T
tx) T
ty

trunc :: Exponent -> T -> T
trunc :: Basis -> T -> T
trunc Basis
bnd (Basis
xe, Mantissa
xm) =
   if Basis
bnd Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
> Basis
xe
     then (Basis
bnd, [])
     else (Basis
xe, Basis -> Mantissa -> Mantissa
forall a. Basis -> [a] -> [a]
take (Basis
1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
bnd) Mantissa
xm)

equalApprox :: Basis -> Exponent -> T -> T -> Bool
equalApprox :: Basis -> Basis -> T -> T -> Bool
equalApprox Basis
b Basis
bnd T
x T
y =
   T -> Basis
forall a b. (a, b) -> a
fst (Basis -> T -> T
trimUntil Basis
bnd (Basis -> T -> T -> T
sub Basis
b T
x T
y)) Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
== Basis
bnd


{- |
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@.
-}
ifLazy :: Basis -> Bool -> T -> T -> T
ifLazy :: Basis -> Bool -> T -> T -> T
ifLazy Basis
b Bool
c x :: T
x@(Basis
xe, Mantissa
_) y :: T
y@(Basis
ye, Mantissa
_) =
   let ze :: Basis
ze = Basis -> Basis -> Basis
forall a. Ord a => a -> a -> a
max Basis
xe Basis
ye
       xm :: Mantissa
xm = Basis -> Basis -> T -> Mantissa
alignMant Basis
b Basis
ze T
x
       ym :: Mantissa
ym = Basis -> Basis -> T -> Mantissa
alignMant Basis
b Basis
ze T
y
       recurse :: Mantissa -> Mantissa -> Mantissa
       recurse :: Mantissa -> Mantissa -> Mantissa
recurse Mantissa
xs0 Mantissa
ys0 =
          Mantissa
-> Mantissa -> Mantissa -> (T -> T -> Mantissa) -> Mantissa
forall a. Mantissa -> Mantissa -> a -> (T -> T -> a) -> a
withTwoMantissas Mantissa
xs0 Mantissa
ys0 [] ((T -> T -> Mantissa) -> Mantissa)
-> (T -> T -> Mantissa) -> Mantissa
forall a b. (a -> b) -> a -> b
$ \(Basis
x0,Mantissa
xs1) (Basis
y0,Mantissa
ys1) ->
          if Basis -> Basis
forall a. C a => a -> a
abs (Basis
y0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
x0) Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
> Basis
2
            then if Bool
c then Mantissa
xs0 else Mantissa
ys0
            else
              {-
              @x0==y0 || c@ means that in case of @x0==y0@
              we do not have to check @c@.
              -}
              Mantissa
-> Mantissa -> Mantissa -> (T -> T -> Mantissa) -> Mantissa
forall a. Mantissa -> Mantissa -> a -> (T -> T -> a) -> a
withTwoMantissas Mantissa
xs1 Mantissa
ys1 ((if Basis
x0Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
==Basis
y0 Bool -> Bool -> Bool
|| Bool
c then Basis
x0 else Basis
y0) Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: []) ((T -> T -> Mantissa) -> Mantissa)
-> (T -> T -> Mantissa) -> Mantissa
forall a b. (a -> b) -> a -> b
$
                  \(Basis
x1,Mantissa
xs2) (Basis
y1,Mantissa
ys2) ->
                {-
                We can choose @z0@ only when knowing also x1 and y1.
                Because of x0x1 = 09 and y0y1 = 19
                we may always choose the larger one of x0 and y0
                in order to get admissible digit z1.
                But this would be wrong for x0x1 = 0(-9) and y0y1 = 1(-9).
                -}
                let z0 :: Basis
z0  = Basis -> (Basis, Basis) -> (Basis, Basis) -> Basis
mean2 Basis
b (Basis
x0,Basis
x1) (Basis
y0,Basis
y1)
                    x1' :: Basis
x1' = Basis
x1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+(Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
z0)Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
b
                    y1' :: Basis
y1' = Basis
y1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+(Basis
y0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
z0)Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
b
                in  if Basis -> Basis
forall a. C a => a -> a
abs Basis
x1' Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
< Basis
b  Bool -> Bool -> Bool
&&  Basis -> Basis
forall a. C a => a -> a
abs Basis
y1' Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
< Basis
b
                      then Basis
z0 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Mantissa -> Mantissa
recurse (Basis
x1'Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
xs2) (Basis
y1'Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
ys2)
                      else if Bool
c then Mantissa
xs0 else Mantissa
ys0
   in  (Basis
ze, Mantissa -> Mantissa -> Mantissa
recurse Mantissa
xm Mantissa
ym)

{- |
> 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@
-}
{-# INLINE mean2 #-}
mean2 :: Basis -> (Digit,Digit) -> (Digit,Digit) -> Digit
mean2 :: Basis -> (Basis, Basis) -> (Basis, Basis) -> Basis
mean2 Basis
b (Basis
x0,Basis
x1) (Basis
y0,Basis
y1) =
   ((Basis
x0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
y0Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1)Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
b Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+ (Basis
x1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
y1)) Basis -> Basis -> Basis
forall a. C a => a -> a -> a
`div` (Basis
2Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
b)

{-
In a first trial I used

> zipMantissas :: Mantissa -> Mantissa -> [(Digit, Digit)]

for implementation of ifLazy.
However, this required to extract digits from the pairs
after the decision for an argument.
With withTwoMantissas we can just return a pointer to the original list.
-}
withTwoMantissas ::
   Mantissa -> Mantissa ->
   a ->
   ((Digit,Mantissa) -> (Digit,Mantissa) -> a) ->
   a
withTwoMantissas :: Mantissa -> Mantissa -> a -> (T -> T -> a) -> a
withTwoMantissas [] [] a
r T -> T -> a
_ = a
r
withTwoMantissas [] (Basis
y:Mantissa
ys) a
_ T -> T -> a
f = T -> T -> a
f (Basis
0,[]) (Basis
y,Mantissa
ys)
withTwoMantissas (Basis
x:Mantissa
xs) [] a
_ T -> T -> a
f = T -> T -> a
f (Basis
x,Mantissa
xs) (Basis
0,[])
withTwoMantissas (Basis
x:Mantissa
xs) (Basis
y:Mantissa
ys) a
_ T -> T -> a
f = T -> T -> a
f (Basis
x,Mantissa
xs) (Basis
y,Mantissa
ys)


align :: Basis -> Exponent -> T -> T
align :: Basis -> Basis -> T -> T
align Basis
b Basis
ye T
x = (Basis
ye, Basis -> Basis -> T -> Mantissa
alignMant Basis
b Basis
ye T
x)

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

@x@ and @(e, alignMant b e x)@ represent the same number.
-}
alignMant :: Basis -> Exponent -> T -> Mantissa
alignMant :: Basis -> Basis -> T -> Mantissa
alignMant Basis
b Basis
e (Basis
xe,Mantissa
xm) =
   if Basis
eBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
xe
     then Basis -> Basis -> Mantissa
forall a. Basis -> a -> [a]
replicate (Basis
eBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
xe) Basis
0 Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Mantissa
xm
     else let (Mantissa
xm0,Mantissa
xm1) = Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
eBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1) Mantissa
xm
          in  Basis -> Mantissa -> Basis
forall a. C a => Basis -> Mantissa -> a
mantissaToNum Basis
b Mantissa
xm0 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
xm1

absolute :: T -> T
absolute :: T -> T
absolute (Basis
xe,Mantissa
xm) = (Basis
xe, Mantissa -> Mantissa
absMant Mantissa
xm)

absMant :: Mantissa -> Mantissa
absMant :: Mantissa -> Mantissa
absMant xa :: Mantissa
xa@(Basis
x:Mantissa
xs) =
   case Basis -> Basis -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Basis
x Basis
0 of
      Ordering
EQ -> Basis
x Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa -> Mantissa
absMant Mantissa
xs
      Ordering
LT -> Mantissa -> Mantissa
forall a. C a => [a] -> [a]
Poly.negate Mantissa
xa
      Ordering
GT -> Mantissa
xa
absMant [] = []


{- * arithmetic -}

fromLaurent :: LPoly.T Digit -> T
fromLaurent :: T Basis -> T
fromLaurent (LPoly.Cons Basis
nxe Mantissa
xm) = (Basis -> Basis
forall a. C a => a -> a
NP.negate Basis
nxe, Mantissa
xm)

toLaurent :: T -> LPoly.T Digit
toLaurent :: T -> T Basis
toLaurent (Basis
xe, Mantissa
xm) = Basis -> Mantissa -> T Basis
forall a. Basis -> [a] -> T a
LPoly.Cons (Basis -> Basis
forall a. C a => a -> a
NP.negate Basis
xe) Mantissa
xm

liftLaurent2 ::
   (LPoly.T Digit -> LPoly.T Digit -> LPoly.T Digit) ->
      (T -> T -> T)
liftLaurent2 :: (T Basis -> T Basis -> T Basis) -> T -> T -> T
liftLaurent2 T Basis -> T Basis -> T Basis
f T
x T
y =
   T Basis -> T
fromLaurent (T Basis -> T Basis -> T Basis
f (T -> T Basis
toLaurent T
x) (T -> T Basis
toLaurent T
y))

liftLaurentMany ::
   ([LPoly.T Digit] -> LPoly.T Digit) ->
      ([T] -> T)
liftLaurentMany :: ([T Basis] -> T Basis) -> [T] -> T
liftLaurentMany [T Basis] -> T Basis
f =
   T Basis -> T
fromLaurent (T Basis -> T) -> ([T] -> T Basis) -> [T] -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [T Basis] -> T Basis
f ([T Basis] -> T Basis) -> ([T] -> [T Basis]) -> [T] -> T Basis
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (T -> T Basis) -> [T] -> [T Basis]
forall a b. (a -> b) -> [a] -> [b]
map T -> T Basis
toLaurent

{- |
Add two numbers but do not eliminate leading zeros.
-}
add :: Basis -> T -> T -> T
add :: Basis -> T -> T -> T
add Basis
b T
x T
y = Basis -> T -> T
compress Basis
b ((T Basis -> T Basis -> T Basis) -> T -> T -> T
liftLaurent2 T Basis -> T Basis -> T Basis
forall a. C a => T a -> T a -> T a
LPoly.add T
x T
y)

sub :: Basis -> T -> T -> T
sub :: Basis -> T -> T -> T
sub Basis
b T
x T
y = Basis -> T -> T
compress Basis
b ((T Basis -> T Basis -> T Basis) -> T -> T -> T
liftLaurent2 T Basis -> T Basis -> T Basis
forall a. C a => T a -> T a -> T a
LPoly.sub T
x T
y)

neg :: Basis -> T -> T
neg :: Basis -> T -> T
neg Basis
_ (Basis
xe, Mantissa
xm) = (Basis
xe, Mantissa -> Mantissa
forall a. C a => [a] -> [a]
Poly.negate Mantissa
xm)


{- |
Add at most @basis@ summands.
More summands will violate the allowed digit range.
-}
addSome :: Basis -> [T] -> T
addSome :: Basis -> [T] -> T
addSome Basis
b = Basis -> T -> T
compress Basis
b (T -> T) -> ([T] -> T) -> [T] -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([T Basis] -> T Basis) -> [T] -> T
liftLaurentMany [T Basis] -> T Basis
forall a. C a => [a] -> a
sum

{- |
Add many numbers efficiently by computing sums of sub lists
with only little carry propagation.
-}
addMany :: Basis -> [T] -> T
addMany :: Basis -> [T] -> T
addMany Basis
_ [] = T
zero
addMany Basis
b [T]
ys =
   let recurse :: [T] -> T
recurse [T]
xs =
          case ([T] -> T) -> [[T]] -> [T]
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> [T] -> T
addSome Basis
b) (Basis -> [T] -> [[T]]
forall a. Basis -> [a] -> [[a]]
sliceVertical Basis
b [T]
xs) of
            [T
s]  -> T
s
            [T]
sums -> [T] -> T
recurse [T]
sums
   in  [T] -> T
recurse [T]
ys


type Series = [(Exponent, T)]

{- |
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.
-}
series :: Basis -> Series -> T
series :: Basis -> Series -> T
series Basis
_ [] = [Char] -> T
forall a. HasCallStack => [Char] -> a
error [Char]
"empty series: don't know a good exponent"
-- series _ [] = (0,[]) -- unfortunate choice of exponent
series Basis
b Series
summands =
   {- Some pre-processing that asserts decreasing exponents.
      Increasing coefficients can appear legally
      due to non-unique number representation. -}
   let (Mantissa
es,[T]
xs)    = Series -> (Mantissa, [T])
forall a b. [(a, b)] -> ([a], [b])
unzip Series
summands
       safeSeries :: Series
safeSeries = Mantissa -> [T] -> Series
forall a b. [a] -> [b] -> [(a, b)]
zip ((Basis -> Basis -> Basis) -> Mantissa -> Mantissa
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Basis -> Basis -> Basis
forall a. Ord a => a -> a -> a
min Mantissa
es) [T]
xs
   in  Basis -> Series -> T
seriesPlain Basis
b Series
safeSeries

seriesPlain :: Basis -> Series -> T
seriesPlain :: Basis -> Series -> T
seriesPlain Basis
_ [] = [Char] -> T
forall a. HasCallStack => [Char] -> a
error [Char]
"empty series: don't know a good exponent"
seriesPlain Basis
b Series
summands =
   let (Mantissa
es,Mantissa
m:[Mantissa]
ms) = [T] -> (Mantissa, [Mantissa])
forall a b. [(a, b)] -> ([a], [b])
unzip (((Basis, T) -> T) -> Series -> [T]
forall a b. (a -> b) -> [a] -> [b]
map ((Basis -> T -> T) -> (Basis, T) -> T
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (Basis -> Basis -> T -> T
align Basis
b)) Series
summands)
       eDifs :: Mantissa
eDifs     = (Basis -> Basis -> Basis) -> Mantissa -> Mantissa
forall a b. (a -> a -> b) -> [a] -> [b]
mapAdjacent (-) Mantissa
es
       eDifLists :: [Mantissa]
eDifLists = Basis -> Mantissa -> [Mantissa]
forall a. Basis -> [a] -> [[a]]
sliceVertical (Basis -> Basis
forall a. Enum a => a -> a
pred Basis
b) Mantissa
eDifs
       mLists :: [[Mantissa]]
mLists    = Basis -> [Mantissa] -> [[Mantissa]]
forall a. Basis -> [a] -> [[a]]
sliceVertical (Basis -> Basis
forall a. Enum a => a -> a
pred Basis
b) [Mantissa]
ms
       accum :: Mantissa -> (Mantissa, [Mantissa]) -> (Mantissa, Mantissa)
accum Mantissa
sumM (Mantissa
eDifList,[Mantissa]
mList) =
          let subM :: Mantissa
subM = Mantissa -> [Mantissa] -> Mantissa
forall a. C a => Mantissa -> [[a]] -> [a]
LPoly.addShiftedMany Mantissa
eDifList (Mantissa
sumMMantissa -> [Mantissa] -> [Mantissa]
forall a. a -> [a] -> [a]
:[Mantissa]
mList)
              -- lazy unary sum
              len :: [()]
len = (Basis -> [()]) -> Mantissa -> [()]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((Basis -> () -> [()]) -> () -> Basis -> [()]
forall a b c. (a -> b -> c) -> b -> a -> c
flip Basis -> () -> [()]
forall a. Basis -> a -> [a]
replicate ()) Mantissa
eDifList
              (Mantissa
high,Mantissa
low)  = [()] -> Mantissa -> (Mantissa, Mantissa)
splitAtMatchPadZero [()]
len Mantissa
subM
          {-
          'compressMant' looks unsafe
          when the residue doesn't decrease for many summands.
          Then there is a leading digit of a chunk
          which is not compressed for long time.
          However, if the residue is estimated correctly
          there can be no overflow.
          -}
          in  (Basis -> Mantissa -> Mantissa
compressMant Basis
b Mantissa
low, Mantissa
high)
       (Mantissa
trailingDigits, [Mantissa]
chunks) =
          (Mantissa -> (Mantissa, [Mantissa]) -> (Mantissa, Mantissa))
-> Mantissa -> [(Mantissa, [Mantissa])] -> (Mantissa, [Mantissa])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL Mantissa -> (Mantissa, [Mantissa]) -> (Mantissa, Mantissa)
accum Mantissa
m ([Mantissa] -> [[Mantissa]] -> [(Mantissa, [Mantissa])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Mantissa]
eDifLists [[Mantissa]]
mLists)
   in  Basis -> T -> T
compress Basis
b (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
es, [Mantissa] -> Mantissa
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [Mantissa]
chunks Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Mantissa
trailingDigits)

{-
An alternative series implementation
could reduce carries by do the following cycle
(split, add sub-lists).
This would reduce carries to the minimum
but we must work hard in order to find out lazily
how many digits can be emitted.
-}


{- |
Like 'splitAt',
but it pads with zeros if the list is too short.
This way it preserves
 @ length (fst (splitAtPadZero n xs)) == n @
-}
splitAtPadZero :: Int -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero :: Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero Basis
n [] = (Basis -> Basis -> Mantissa
forall a. Basis -> a -> [a]
replicate Basis
n Basis
0, [])
splitAtPadZero Basis
0 Mantissa
xs = ([], Mantissa
xs)
splitAtPadZero Basis
n (Basis
x:Mantissa
xs) =
   let (Mantissa
ys, Mantissa
zs) = Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero (Basis
nBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) Mantissa
xs
   in  (Basis
xBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
ys, Mantissa
zs)
-- must get a case for negative index

splitAtMatchPadZero :: [()] -> Mantissa -> (Mantissa, Mantissa)
splitAtMatchPadZero :: [()] -> Mantissa -> (Mantissa, Mantissa)
splitAtMatchPadZero [()]
n  [] = ([()] -> Basis -> Mantissa
forall a b. [a] -> b -> [b]
Match.replicate [()]
n Basis
0, [])
splitAtMatchPadZero [] Mantissa
xs = ([], Mantissa
xs)
splitAtMatchPadZero [()]
n (Basis
x:Mantissa
xs) =
   let (Mantissa
ys, Mantissa
zs) = [()] -> Mantissa -> (Mantissa, Mantissa)
splitAtMatchPadZero ([()] -> [()]
forall a. [a] -> [a]
tail [()]
n) Mantissa
xs
   in  (Basis
xBasis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
:Mantissa
ys, Mantissa
zs)


{- |
help showing series summands
-}
truncSeriesSummands :: Series -> Series
truncSeriesSummands :: Series -> Series
truncSeriesSummands = ((Basis, T) -> (Basis, T)) -> Series -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\(Basis
e,T
x) -> (Basis
e,Basis -> T -> T
trunc (-Basis
20) T
x))



scale :: Basis -> Digit -> T -> T
scale :: Basis -> Basis -> T -> T
scale Basis
b Basis
y T
x = Basis -> T -> T
compress Basis
b (Basis -> T -> T
scaleSimple Basis
y T
x)

{-
scaleSimple :: ToInteger.C a => a -> T -> T
scaleSimple y (xe, xm) = (xe, Poly.scale (fromIntegral y) xm)
-}

scaleSimple :: Basis -> T -> T
scaleSimple :: Basis -> T -> T
scaleSimple Basis
y (Basis
xe, Mantissa
xm) = (Basis
xe, Basis -> Mantissa -> Mantissa
forall a. C a => a -> [a] -> [a]
Poly.scale Basis
y Mantissa
xm)

scaleMant :: Basis -> Digit -> Mantissa -> Mantissa
scaleMant :: Basis -> Basis -> Mantissa -> Mantissa
scaleMant Basis
b Basis
y Mantissa
xm = Basis -> Mantissa -> Mantissa
compressMant Basis
b (Basis -> Mantissa -> Mantissa
forall a. C a => a -> [a] -> [a]
Poly.scale Basis
y Mantissa
xm)



mulSeries :: Basis -> T -> T -> Series
mulSeries :: Basis -> T -> T -> Series
mulSeries Basis
_ (Basis
xe,[]) (Basis
ye,Mantissa
_) = [(Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
ye, T
zero)]
mulSeries Basis
b (Basis
xe,Mantissa
xm) (Basis
ye,Mantissa
ym) =
   let zes :: Mantissa
zes = (Basis -> Basis) -> Basis -> Mantissa
forall a. (a -> a) -> a -> [a]
iterate Basis -> Basis
forall a. Enum a => a -> a
pred (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
yeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
1)
       zs :: [T]
zs  = (Basis -> Basis -> T) -> Mantissa -> Mantissa -> [T]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Basis
xd Basis
e -> Basis -> Basis -> T -> T
scale Basis
b Basis
xd (Basis
e,Mantissa
ym)) Mantissa
xm (Mantissa -> Mantissa
forall a. [a] -> [a]
tail Mantissa
zes)
   in  Mantissa -> [T] -> Series
forall a b. [a] -> [b] -> [(a, b)]
zip Mantissa
zes [T]
zs

{- |
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.
-}
mul :: Basis -> T -> T -> T
mul :: Basis -> T -> T -> T
mul Basis
b T
x T
y = T -> T
trimOnce (Basis -> Series -> T
seriesPlain Basis
b (Basis -> T -> T -> Series
mulSeries Basis
b T
x T
y))



{- |
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.
-}
divide :: Basis -> T -> T -> T
divide :: Basis -> T -> T -> T
divide Basis
b (Basis
xe,Mantissa
xm) (Basis
ye',Mantissa
ym') =
   let (Basis
ye,Mantissa
ym) = (T -> Bool) -> (T -> T) -> T -> T
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
>=Basis
b) (Basis -> Bool) -> (T -> Basis) -> T -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Basis -> Basis
forall a. C a => a -> a
abs (Basis -> Basis) -> (T -> Basis) -> T -> Basis
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Mantissa -> Basis
forall a. [a] -> a
head (Mantissa -> Basis) -> (T -> Mantissa) -> T -> Basis
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T -> Mantissa
forall a b. (a, b) -> b
snd)
                       (Basis -> T -> T
decreaseExp Basis
b)
                       (Basis
ye',Mantissa
ym')
   in  if Mantissa -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null Mantissa
xm
         then (Basis
xe,Mantissa
xm)
         else Basis -> (T -> T) -> T -> T
forall a. Basis -> (a -> a) -> a -> a
nest Basis
3 T -> T
trimOnce (Basis -> T -> T
compress Basis
b (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
ye, Basis -> Mantissa -> Mantissa -> Mantissa
divMant Basis
b Mantissa
ym Mantissa
xm))

divMant :: Basis -> Mantissa -> Mantissa -> Mantissa
divMant :: Basis -> Mantissa -> Mantissa -> Mantissa
divMant Basis
_ [] Mantissa
_   = [Char] -> Mantissa
forall a. HasCallStack => [Char] -> a
error [Char]
"Number.Positional: division by zero"
divMant Basis
b Mantissa
ym Mantissa
xm0 =
   (Mantissa, Mantissa) -> Mantissa
forall a b. (a, b) -> b
snd ((Mantissa, Mantissa) -> Mantissa)
-> (Mantissa, Mantissa) -> Mantissa
forall a b. (a -> b) -> a -> b
$
   (Mantissa -> Bool -> (Mantissa, Basis))
-> Mantissa -> [Bool] -> (Mantissa, Mantissa)
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL
      (\Mantissa
xm Bool
fullCompress ->
       let z :: Basis
z = Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
xm) (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
ym)
           {- 'scaleMant' contains compression,
              which is not much of a problem,
              because it is always applied to @ym@.
              That is, there is no nested call. -}
           dif :: Mantissa
dif = Basis -> Mantissa -> Mantissa
pumpFirst Basis
b (Mantissa -> Mantissa -> Mantissa
forall a. C a => [a] -> [a] -> [a]
Poly.sub Mantissa
xm (Basis -> Basis -> Mantissa -> Mantissa
scaleMant Basis
b Basis
z Mantissa
ym))
           cDif :: Mantissa
cDif = if Bool
fullCompress
                    then Basis -> Mantissa -> Mantissa
compressMant       Basis
b Mantissa
dif
                    else Basis -> Mantissa -> Mantissa
compressSecondMant Basis
b Mantissa
dif
       in  (Mantissa
cDif, Basis
z))
   Mantissa
xm0 ([Bool] -> [Bool]
forall a. [a] -> [a]
cycle (Basis -> Bool -> [Bool]
forall a. Basis -> a -> [a]
replicate (Basis
bBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
1) Bool
False [Bool] -> [Bool] -> [Bool]
forall a. [a] -> [a] -> [a]
++ [Bool
True]))

divMantSlow :: Basis -> Mantissa -> Mantissa -> Mantissa
divMantSlow :: Basis -> Mantissa -> Mantissa -> Mantissa
divMantSlow Basis
_ [] = [Char] -> Mantissa -> Mantissa
forall a. HasCallStack => [Char] -> a
error [Char]
"Number.Positional: division by zero"
divMantSlow Basis
b Mantissa
ym =
   (Mantissa -> Maybe T) -> Mantissa -> Mantissa
forall b a. (b -> Maybe (a, b)) -> b -> [a]
List.unfoldr
      (\Mantissa
xm ->
       let z :: Basis
z = Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
xm) (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
ym)
           d :: Mantissa
d = Basis -> Mantissa -> Mantissa
compressMant Basis
b (Basis -> Mantissa -> Mantissa
pumpFirst Basis
b
                  (Mantissa -> Mantissa -> Mantissa
forall a. C a => [a] -> [a] -> [a]
Poly.sub Mantissa
xm (Basis -> Mantissa -> Mantissa
forall a. C a => a -> [a] -> [a]
Poly.scale Basis
z Mantissa
ym)))
       in  T -> Maybe T
forall a. a -> Maybe a
Just (Basis
z, Mantissa
d))

reciprocal :: Basis -> T -> T
reciprocal :: Basis -> T -> T
reciprocal Basis
b = Basis -> T -> T -> T
divide Basis
b T
one


{- |
Fast division for small integral divisors,
which occur for instance in summands of power series.
-}
divIntMant :: Basis -> Digit -> Mantissa -> Mantissa
divIntMant :: Basis -> Basis -> Mantissa -> Mantissa
divIntMant Basis
b Basis
y Mantissa
xInit =
   (T -> Maybe (Basis, T)) -> T -> Mantissa
forall b a. (b -> Maybe (a, b)) -> b -> [a]
List.unfoldr (\(Basis
r,Mantissa
rxs) ->
             let rb :: Basis
rb = Basis
rBasis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
b
                 (Basis
rbx, Mantissa
xs', Bool
run) =
                    case Mantissa
rxs of
                       []     -> (Basis
rb,   [], Basis
rBasis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
/=Basis
0)
                       (x:xs) -> (Basis
rbBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
x, Mantissa
xs, Bool
True)
                 (Basis
d,Basis
m) = Basis -> Basis -> (Basis, Basis)
forall a. C a => a -> a -> (a, a)
divMod Basis
rbx Basis
y
             in  Bool -> (Basis, T) -> Maybe (Basis, T)
forall a. Bool -> a -> Maybe a
toMaybe Bool
run (Basis
d, (Basis
m, Mantissa
xs')))
           (Basis
0,Mantissa
xInit)

-- this version is simple but ignores the possibility of a terminating result
divIntMantInf :: Basis -> Digit -> Mantissa -> Mantissa
divIntMantInf :: Basis -> Basis -> Mantissa -> Mantissa
divIntMantInf Basis
b Basis
y =
   ((Basis, Basis) -> Basis) -> [(Basis, Basis)] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map (Basis, Basis) -> Basis
forall a b. (a, b) -> a
fst ([(Basis, Basis)] -> Mantissa)
-> (Mantissa -> [(Basis, Basis)]) -> Mantissa -> Mantissa
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Basis, Basis)] -> [(Basis, Basis)]
forall a. [a] -> [a]
tail ([(Basis, Basis)] -> [(Basis, Basis)])
-> (Mantissa -> [(Basis, Basis)]) -> Mantissa -> [(Basis, Basis)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
      ((Basis, Basis) -> Basis -> (Basis, Basis))
-> (Basis, Basis) -> Mantissa -> [(Basis, Basis)]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\(Basis
_,Basis
r) Basis
x -> Basis -> Basis -> (Basis, Basis)
forall a. C a => a -> a -> (a, a)
divMod (Basis
rBasis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
bBasis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
x) Basis
y) (Basis
forall a. HasCallStack => a
undefined,Basis
0) (Mantissa -> [(Basis, Basis)])
-> (Mantissa -> Mantissa) -> Mantissa -> [(Basis, Basis)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
         (Mantissa -> Mantissa -> Mantissa
forall a. [a] -> [a] -> [a]
++ Basis -> Mantissa
forall a. a -> [a]
repeat Basis
0)

divInt :: Basis -> Digit -> T -> T
divInt :: Basis -> Basis -> T -> T
divInt Basis
b Basis
y (Basis
xe,Mantissa
xm) =
   -- (xe, divIntMant b y xm)
   let z :: T
z  = (Basis
xe, Basis -> Basis -> Mantissa -> Mantissa
divIntMant Basis
b Basis
y Mantissa
xm)
       {- Division by big integers may cause leading zeros.
          Eliminate as many as we can expect from the division.
          If the input number has leading zeros (it might be equal to zero),
          then the result may have, too. -}
       tz :: (Basis, T)
tz = ((Basis, T) -> Bool)
-> ((Basis, T) -> (Basis, T)) -> (Basis, T) -> (Basis, T)
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((Basis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<=Basis
1) (Basis -> Bool) -> ((Basis, T) -> Basis) -> (Basis, T) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Basis, T) -> Basis
forall a b. (a, b) -> a
fst) (\(Basis
yi,T
zi) -> (Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div Basis
yi Basis
b, T -> T
trimOnce T
zi)) (Basis
y,T
z)
   in  (Basis, T) -> T
forall a b. (a, b) -> b
snd (Basis, T)
tz


{- * algebraic functions -}


sqrt :: Basis -> T -> T
sqrt :: Basis -> T -> T
sqrt Basis
b = Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T
sqrtDriver Basis
b Basis -> FixedPoint -> Mantissa
sqrtFP

sqrtDriver :: Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T
sqrtDriver :: Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T
sqrtDriver Basis
_ Basis -> FixedPoint -> Mantissa
_ (Basis
xe,[]) = (Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div Basis
xe Basis
2, [])
sqrtDriver Basis
b Basis -> FixedPoint -> Mantissa
sqrtFPworker T
x =
   let (Basis
exe,Basis
ex0:Mantissa
exm) = if Basis -> Bool
forall a. (C a, C a) => a -> Bool
odd (T -> Basis
forall a b. (a, b) -> a
fst T
x) then Basis -> T -> T
decreaseExp Basis
b T
x else T
x
       (Basis
nxe,(Integer
nx0,Mantissa
nxm)) =
           ((Basis, FixedPoint) -> Bool)
-> ((Basis, FixedPoint) -> (Basis, FixedPoint))
-> (Basis, FixedPoint)
-> (Basis, FixedPoint)
forall a. (a -> Bool) -> (a -> a) -> a -> a
until (\(Basis, FixedPoint)
xi -> FixedPoint -> Integer
forall a b. (a, b) -> a
fst ((Basis, FixedPoint) -> FixedPoint
forall a b. (a, b) -> b
snd (Basis, FixedPoint)
xi) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
b Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^ Integer
2)
                 (Basis
-> ((Basis, FixedPoint) -> (Basis, FixedPoint))
-> (Basis, FixedPoint)
-> (Basis, FixedPoint)
forall a. Basis -> (a -> a) -> a -> a
nest Basis
2 (Basis -> (Basis, FixedPoint) -> (Basis, FixedPoint)
decreaseExpFP Basis
b))
                 (Basis
exe, (Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
ex0, Mantissa
exm))
   in  Basis -> T -> T
compress Basis
b (Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div Basis
nxe Basis
2, Basis -> FixedPoint -> Mantissa
sqrtFPworker Basis
b (Integer
nx0,Mantissa
nxm))


sqrtMant :: Basis -> Mantissa -> Mantissa
sqrtMant :: Basis -> Mantissa -> Mantissa
sqrtMant Basis
_ [] = []
sqrtMant Basis
b (Basis
x:Mantissa
xs) =
   Basis -> FixedPoint -> Mantissa
sqrtFP Basis
b (Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
x, Mantissa
xs)

{- |
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) `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 @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)@ ?
-}
sqrtFP :: Basis -> FixedPoint -> Mantissa
sqrtFP :: Basis -> FixedPoint -> Mantissa
sqrtFP Basis
b (Integer
x0,Mantissa
xs) =
   let y0 :: Basis
y0   = Double -> Basis
forall a b. (C a, C b) => a -> b
round (Double -> Double
forall a. C a => a -> a
NP.sqrt (Integer -> Double
forall a. C a => Integer -> a
fromInteger Integer
x0 :: Double))
       dyx0 :: Basis
dyx0 = Integer -> Basis
forall a. C a => Integer -> a
fromInteger (Integer
x0 Integer -> Integer -> Integer
forall a. C a => a -> a -> a
- Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
y0 Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^ Integer
2)

       accum :: Mantissa -> T -> (Mantissa, Basis)
accum Mantissa
dif (Basis
e,Mantissa
ty) =
          -- (e,ty) == xm - (trunc j y)^2
          let yj :: Basis
yj = Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div (Mantissa -> Basis
forall a. [a] -> a
head Mantissa
dif Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+ Basis
y0) (Basis
2Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
y0)
              newDif :: Mantissa
newDif = Basis -> Mantissa -> Mantissa
pumpFirst Basis
b (Mantissa -> Mantissa) -> Mantissa -> Mantissa
forall a b. (a -> b) -> a -> b
$
                 Basis -> Mantissa -> Mantissa -> Mantissa
forall a. C a => Basis -> [a] -> [a] -> [a]
LPoly.addShifted Basis
e
                    (Mantissa -> Mantissa -> Mantissa
forall a. C a => [a] -> [a] -> [a]
Poly.sub Mantissa
dif (Basis -> Basis -> Mantissa -> Mantissa
scaleMant Basis
b (Basis
2Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
yj) Mantissa
ty))
                    [Basis
yjBasis -> Basis -> Basis
forall a. C a => a -> a -> a
*Basis
yj]
              {- We could always compress the full difference number,
                 but it is not necessary,
                 and we save dependencies on less significant digits. -}
              cNewDif :: Mantissa
cNewDif =
                 if Basis -> Basis -> Basis
forall a. C a => a -> a -> a
mod Basis
e Basis
b Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
== Basis
0
                   then Basis -> Mantissa -> Mantissa
compressMant       Basis
b Mantissa
newDif
                   else Basis -> Mantissa -> Mantissa
compressSecondMant Basis
b Mantissa
newDif
          in  (Mantissa
cNewDif, Basis
yj)

       truncs :: [Mantissa]
truncs = Mantissa -> [Mantissa]
forall a. [a] -> [[a]]
lazyInits Mantissa
y
       y :: Mantissa
y = Basis
y0 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: (Mantissa, Mantissa) -> Mantissa
forall a b. (a, b) -> b
snd ((Mantissa -> T -> (Mantissa, Basis))
-> Mantissa -> [T] -> (Mantissa, Mantissa)
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL
                        Mantissa -> T -> (Mantissa, Basis)
accum
                        (Basis -> Mantissa -> Mantissa
pumpFirst Basis
b (Basis
dyx0 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
xs))
                        (Mantissa -> [Mantissa] -> [T]
forall a b. [a] -> [b] -> [(a, b)]
zip [Basis
1..] (Basis -> [Mantissa] -> [Mantissa]
forall a. Basis -> [a] -> [a]
drop Basis
2 [Mantissa]
truncs)))
   in  Mantissa
y


sqrtNewton :: Basis -> T -> T
sqrtNewton :: Basis -> T -> T
sqrtNewton Basis
b = Basis -> (Basis -> FixedPoint -> Mantissa) -> T -> T
sqrtDriver Basis
b Basis -> FixedPoint -> Mantissa
sqrtFPNewton

{- |
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.@
-}
sqrtFPNewton :: Basis -> FixedPoint -> Mantissa
sqrtFPNewton :: Basis -> FixedPoint -> Mantissa
sqrtFPNewton Basis
bInt (Integer
x0,Mantissa
xs) =
   let b :: Integer
b = Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
bInt
       chunkLengths :: Mantissa
chunkLengths = (Basis -> Basis) -> Basis -> Mantissa
forall a. (a -> a) -> a -> [a]
iterate (Basis
2Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*) Basis
1
       xChunks :: [Integer]
xChunks = (Mantissa -> Integer) -> [Mantissa] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Basis -> Mantissa -> Integer
forall a. C a => Basis -> Mantissa -> a
mantissaToNum Basis
bInt) ([Mantissa] -> [Integer]) -> [Mantissa] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Mantissa, [Mantissa]) -> [Mantissa]
forall a b. (a, b) -> b
snd ((Mantissa, [Mantissa]) -> [Mantissa])
-> (Mantissa, [Mantissa]) -> [Mantissa]
forall a b. (a -> b) -> a -> b
$
            (Mantissa -> Basis -> (Mantissa, Mantissa))
-> Mantissa -> Mantissa -> (Mantissa, [Mantissa])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL (\Mantissa
x Basis
cl -> (Mantissa, Mantissa) -> (Mantissa, Mantissa)
forall a b. (a, b) -> (b, a)
swap (Basis -> Mantissa -> (Mantissa, Mantissa)
splitAtPadZero Basis
cl Mantissa
x))
                           Mantissa
xs Mantissa
chunkLengths
       basisPowers :: [Integer]
basisPowers = (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2) Integer
b
       truncXs :: [Integer]
truncXs = (Integer -> (Integer, Integer) -> Integer)
-> Integer -> [(Integer, Integer)] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Integer
acc (Integer
bp,Integer
frac) -> Integer
accInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
bpInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
frac) Integer
x0
                       ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer]
basisPowers [Integer]
xChunks)
       accum :: Integer -> (Integer, Integer) -> (Integer, Integer)
accum Integer
y (Integer
bp, Integer
x) =
          let ybp :: Integer
ybp  = Integer
y Integer -> Integer -> Integer
forall a. C a => a -> a -> a
* Integer
bp
              newY :: Integer
newY = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
ybp Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
x Integer -> Integer -> Integer
forall a. C a => a -> a -> a
* Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
bp Integer
b) Integer
y) Integer
2
          in  (Integer
newY, Integer
newYInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
ybp)
       y0 :: Integer
y0 = Double -> Integer
forall a b. (C a, C b) => a -> b
round (Double -> Double
forall a. C a => a -> a
NP.sqrt (Integer -> Double
forall a. C a => Integer -> a
fromInteger Integer
x0 :: Double))
       yChunks :: [Integer]
yChunks = (Integer, [Integer]) -> [Integer]
forall a b. (a, b) -> b
snd ((Integer, [Integer]) -> [Integer])
-> (Integer, [Integer]) -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> (Integer, Integer) -> (Integer, Integer))
-> Integer -> [(Integer, Integer)] -> (Integer, [Integer])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL Integer -> (Integer, Integer) -> (Integer, Integer)
accum
                         Integer
y0 ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer]
basisPowers ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
truncXs))
       yFrac :: Mantissa
yFrac = [Mantissa] -> Mantissa
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([Mantissa] -> Mantissa) -> [Mantissa] -> Mantissa
forall a b. (a -> b) -> a -> b
$ (Basis -> Integer -> Mantissa)
-> Mantissa -> [Integer] -> [Mantissa]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Basis -> Basis -> Integer -> Mantissa
mantissaFromFixedInt Basis
bInt) Mantissa
chunkLengths [Integer]
yChunks
   in  Integer -> Basis
forall a. C a => Integer -> a
fromInteger Integer
y0 Basis -> Mantissa -> Mantissa
forall a. a -> [a] -> [a]
: Mantissa
yFrac


{- |
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 'List.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],[0,1],[0,1,2],[0,1,2,*** Exception: Prelude.undefined
-}
lazyInits :: [a] -> [[a]]
lazyInits :: [a] -> [[a]]
lazyInits ~(a
x:[a]
xs)  =  [] [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [[a]]
forall a. [a] -> [[a]]
lazyInits [a]
xs)
{-
The lazy match above is irrefutable,
so the pattern @[]@ would never be reached.
-}



{- * transcendent functions -}

{- ** exponential functions -}

expSeries :: Basis -> T -> Series
expSeries :: Basis -> T -> Series
expSeries Basis
b T
xOrig =
   let x :: T
x   = Basis -> T -> T
negativeExp Basis
b T
xOrig
       xps :: [T]
xps = (T -> Basis -> T) -> T -> Mantissa -> [T]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\T
p Basis
n -> Basis -> Basis -> T -> T
divInt Basis
b Basis
n (Basis -> T -> T -> T
mul Basis
b T
x T
p)) T
one [Basis
1..]
   in  (T -> (Basis, T)) -> [T] -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\T
xp -> (T -> Basis
forall a b. (a, b) -> a
fst T
xp, T
xp)) [T]
xps

{- |
Absolute value of argument should be below 1.
-}
expSmall :: Basis -> T -> T
expSmall :: Basis -> T -> T
expSmall Basis
b T
x = Basis -> Series -> T
series Basis
b (Basis -> T -> Series
expSeries Basis
b T
x)


expSeriesLazy :: Basis -> T -> Series
expSeriesLazy :: Basis -> T -> Series
expSeriesLazy Basis
b x :: T
x@(Basis
xe,Mantissa
_) =
   let xps :: [T]
xps = (T -> Basis -> T) -> T -> Mantissa -> [T]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\T
p Basis
n -> Basis -> Basis -> T -> T
divInt Basis
b Basis
n (Basis -> T -> T -> T
mul Basis
b T
x T
p)) T
one [Basis
1..]
       {- much effort for computing the residue exponents
          without touching the arguments mantissa -}
       es :: [Double]
       es :: [Double]
es = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-)
               ((Basis -> Double) -> Mantissa -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral ((Basis -> Basis) -> Basis -> Mantissa
forall a. (a -> a) -> a -> [a]
iterate ((Basis
1Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
xe)Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+) Basis
0))
               ((Double -> Double -> Double) -> Double -> [Double] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Double -> Double -> Double
forall a. C a => a -> a -> a
(+) Double
0
                  ((Integer -> Double) -> [Integer] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. C a => a -> a -> a
logBase (Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
b)
                          (Double -> Double) -> (Integer -> Double) -> Integer -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Double
forall a. C a => Integer -> a
fromInteger) [Integer
1..]))
   in  Mantissa -> [T] -> Series
forall a b. [a] -> [b] -> [(a, b)]
zip ((Double -> Basis) -> [Double] -> Mantissa
forall a b. (a -> b) -> [a] -> [b]
map Double -> Basis
forall a b. (C a, C b) => a -> b
ceiling [Double]
es) [T]
xps

expSmallLazy :: Basis -> T -> T
expSmallLazy :: Basis -> T -> T
expSmallLazy Basis
b T
x = Basis -> Series -> T
series Basis
b (Basis -> T -> Series
expSeriesLazy Basis
b T
x)


exp :: Basis -> T -> T
exp :: Basis -> T -> T
exp Basis
b T
x =
   let (Integer
xInt,Mantissa
xFrac) = Basis -> T -> FixedPoint
toFixedPoint Basis
b (Basis -> T -> T
compress Basis
b T
x)
       yFrac :: T
yFrac = Basis -> T -> T
expSmall Basis
b (-Basis
1,Mantissa
xFrac)
       {-
       (xFrac0,xFrac1) = splitAt 2 xFrac
       yFrac = mul b
                 -- slow convergence but simple argument
                 (expSmall b (-1, xFrac0))
                 -- fast convergence but big argument
                 (expSmall b (-3, xFrac1))
       -}
   in  Basis -> Integer -> T -> T -> T -> T
intPower Basis
b Integer
xInt T
yFrac (Basis -> T
recipEConst Basis
b) (Basis -> T
eConst Basis
b)

intPower :: Basis -> Integer -> T -> T -> T -> T
intPower :: Basis -> Integer -> T -> T -> T -> T
intPower Basis
b Integer
expon T
neutral T
recipX T
x =
   if Integer
expon Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
0
     then Basis -> Integer -> T -> T -> T
cardPower Basis
b   Integer
expon  T
neutral T
x
     else Basis -> Integer -> T -> T -> T
cardPower Basis
b (-Integer
expon) T
neutral T
recipX

cardPower :: Basis -> Integer -> T -> T -> T
cardPower :: Basis -> Integer -> T -> T -> T
cardPower Basis
b Integer
expon T
neutral T
x =
   if Integer
expon Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
0
     then (T -> T -> T) -> T -> T -> Integer -> T
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative (Basis -> T -> T -> T
mul Basis
b) T
neutral T
x Integer
expon
     else [Char] -> T
forall a. HasCallStack => [Char] -> a
error [Char]
"negative exponent - use intPower"


{- |
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.
-}
powerSeries :: Basis -> Rational -> T -> Series
powerSeries :: Basis -> Rational -> T -> Series
powerSeries Basis
b Rational
expon T
xOrig =
   let scaleRat :: Basis -> Rational -> T -> T
scaleRat Basis
ni Rational
yi =
          Basis -> Basis -> T -> T
divInt Basis
b (Integer -> Basis
forall a. C a => Integer -> a
fromInteger (Rational -> Integer
forall a. T a -> a
denominator Rational
yi) Basis -> Basis -> Basis
forall a. C a => a -> a -> a
* Basis
ni) (T -> T) -> (T -> T) -> T -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          Basis -> T -> T
scaleSimple (Integer -> Basis
forall a. C a => Integer -> a
fromInteger (Rational -> Integer
forall a. T a -> a
numerator Rational
yi))
       x :: T
x   = Basis -> T -> T
negativeExp Basis
b (Basis -> T -> T -> T
sub Basis
b T
xOrig T
one)
       xps :: [T]
xps = (T -> (Basis, Rational) -> T) -> T -> [(Basis, Rational)] -> [T]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\T
p (Basis, Rational)
fac -> (Basis -> Rational -> T -> T) -> (Basis, Rational) -> T -> T
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Basis -> Rational -> T -> T
scaleRat (Basis, Rational)
fac (Basis -> T -> T -> T
mul Basis
b T
x T
p))
                   T
one (Mantissa -> [Rational] -> [(Basis, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Basis
1..] ((Rational -> Rational) -> Rational -> [Rational]
forall a. (a -> a) -> a -> [a]
iterate (Rational -> Rational -> Rational
forall a. C a => a -> a -> a
subtract Rational
1) Rational
expon))
   in  (T -> (Basis, T)) -> [T] -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\T
xp -> (T -> Basis
forall a b. (a, b) -> a
fst T
xp, T
xp)) [T]
xps

powerSmall :: Basis -> Rational -> T -> T
powerSmall :: Basis -> Rational -> T -> T
powerSmall Basis
b Rational
y T
x = Basis -> Series -> T
series Basis
b (Basis -> Rational -> T -> Series
powerSeries Basis
b Rational
y T
x)

power :: Basis -> Rational -> T -> T
power :: Basis -> Rational -> T -> T
power Basis
b Rational
expon T
x =
   let num :: Integer
num   = Rational -> Integer
forall a. T a -> a
numerator   Rational
expon
       den :: Integer
den   = Rational -> Integer
forall a. T a -> a
denominator Rational
expon
       rootX :: T
rootX = Basis -> Integer -> T -> T
root Basis
b Integer
den T
x
   in  Basis -> Integer -> T -> T -> T -> T
intPower Basis
b Integer
num T
one (Basis -> T -> T
reciprocal Basis
b T
rootX) T
rootX

root :: Basis -> Integer -> T -> T
root :: Basis -> Integer -> T -> T
root Basis
b Integer
expon T
x =
   let estimate :: T
estimate = Basis -> (Double -> Double) -> T -> T
liftDoubleApprox Basis
b (Double -> Double -> Double
forall a. C a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. C a => a -> a -> a
/ Integer -> Double
forall a. C a => Integer -> a
fromInteger Integer
expon)) T
x
       estPower :: T
estPower = Basis -> Integer -> T -> T -> T
cardPower Basis
b Integer
expon T
one T
estimate
       residue :: T
residue  = Basis -> T -> T -> T
divide Basis
b T
x T
estPower
   in  Basis -> T -> T -> T
mul Basis
b T
estimate (Basis -> Rational -> T -> T
powerSmall Basis
b (Integer
1 Integer -> Integer -> Rational
forall a. C a => a -> a -> T a
% Integer -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Integer
expon) T
residue)



{- |
Absolute value of argument should be below 1.
-}
cosSinhSmall :: Basis -> T -> (T, T)
cosSinhSmall :: Basis -> T -> (T, T)
cosSinhSmall Basis
b T
x =
   let (Series
coshXps, Series
sinhXps) = [((Basis, T), (Basis, T))] -> (Series, Series)
forall a b. [(a, b)] -> ([a], [b])
unzip (Series -> [((Basis, T), (Basis, T))]
forall a. [a] -> [(a, a)]
sliceVertPair (Basis -> T -> Series
expSeries Basis
b T
x))
   in  (Basis -> Series -> T
series Basis
b Series
coshXps,
        Basis -> Series -> T
series Basis
b Series
sinhXps)

{- |
Absolute value of argument should be below 1.
-}
cosSinSmall :: Basis -> T -> (T, T)
cosSinSmall :: Basis -> T -> (T, T)
cosSinSmall Basis
b T
x =
   let (Series
coshXps, Series
sinhXps) = [((Basis, T), (Basis, T))] -> (Series, Series)
forall a b. [(a, b)] -> ([a], [b])
unzip (Series -> [((Basis, T), (Basis, T))]
forall a. [a] -> [(a, a)]
sliceVertPair (Basis -> T -> Series
expSeries Basis
b T
x))
       alternate :: [(a, T)] -> [(a, T)]
alternate [(a, T)]
s =
          (Bool -> (a, T) -> (a, T) -> (a, T))
-> [Bool] -> [(a, T)] -> [(a, T)] -> [(a, T)]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 Bool -> (a, T) -> (a, T) -> (a, T)
forall a. Bool -> a -> a -> a
if' ([Bool] -> [Bool]
forall a. [a] -> [a]
cycle [Bool
True,Bool
False])
             [(a, T)]
s (((a, T) -> (a, T)) -> [(a, T)] -> [(a, T)]
forall a b. (a -> b) -> [a] -> [b]
map (\(a
e,T
y) -> (a
e, Basis -> T -> T
neg Basis
b T
y)) [(a, T)]
s)
   in  (Basis -> Series -> T
series Basis
b (Series -> Series
forall a. [(a, T)] -> [(a, T)]
alternate Series
coshXps),
        Basis -> Series -> T
series Basis
b (Series -> Series
forall a. [(a, T)] -> [(a, T)]
alternate Series
sinhXps))


{- |
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.
-}
cosSinFourth :: Basis -> T -> (T, T)
cosSinFourth :: Basis -> T -> (T, T)
cosSinFourth Basis
b T
x =
   let (T
cosx, T
sinx) = Basis -> T -> (T, T)
cosSinSmall Basis
b (Basis -> Basis -> T -> T
divInt Basis
b Basis
4 T
x)
       sinx2 :: T
sinx2   = Basis -> T -> T -> T
mul Basis
b T
sinx T
sinx
       cosx2 :: T
cosx2   = Basis -> T -> T -> T
mul Basis
b T
cosx T
cosx
       sincosx :: T
sincosx = Basis -> T -> T -> T
mul Basis
b T
sinx T
cosx
   in  (Basis -> T -> T -> T
add Basis
b T
one (Basis -> Basis -> T -> T
scale Basis
b Basis
8 (Basis -> T -> T -> T
mul Basis
b T
cosx2 (Basis -> T -> T -> T
sub Basis
b T
cosx2 T
one))),
        Basis -> Basis -> T -> T
scale Basis
b Basis
4 (Basis -> T -> T -> T
mul Basis
b T
sincosx (Basis -> T -> T -> T
sub Basis
b T
one (Basis -> Basis -> T -> T
scale Basis
b Basis
2 T
sinx2))))


cosSin :: Basis -> T -> (T, T)
cosSin :: Basis -> T -> (T, T)
cosSin Basis
b T
x =
   let pi2 :: T
pi2 = Basis -> Basis -> T -> T
divInt Basis
b Basis
2 (Basis -> T
piConst Basis
b)
       {- @compress@ ensures that the leading digit of the fractional part
          is close to zero -}
       (Integer
quadrant, Mantissa
frac) = Basis -> T -> FixedPoint
toFixedPoint Basis
b (Basis -> T -> T
compress Basis
b (Basis -> T -> T -> T
divide Basis
b T
x T
pi2))
       -- it's possibly faster if we subtract quadrant*pi/4
       wrapped :: T
wrapped = if Integer
quadrantInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0 then T
x else Basis -> T -> T -> T
mul Basis
b T
pi2 (-Basis
1, Mantissa
frac)
       (T
cosW,T
sinW) = Basis -> T -> (T, T)
cosSinSmall Basis
b T
wrapped
   in  case Integer -> Integer -> Integer
forall a. C a => a -> a -> a
mod Integer
quadrant Integer
4 of
          Integer
0 -> (      T
cosW,       T
sinW)
          Integer
1 -> (Basis -> T -> T
neg Basis
b T
sinW,       T
cosW)
          Integer
2 -> (Basis -> T -> T
neg Basis
b T
cosW, Basis -> T -> T
neg Basis
b T
sinW)
          Integer
3 -> (      T
sinW, Basis -> T -> T
neg Basis
b T
cosW)
          Integer
_ -> [Char] -> (T, T)
forall a. HasCallStack => [Char] -> a
error [Char]
"error in implementation of 'mod'"

tan :: Basis -> T -> T
tan :: Basis -> T -> T
tan Basis
b T
x = (T -> T -> T) -> (T, T) -> T
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((T -> T -> T) -> T -> T -> T
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Basis -> T -> T -> T
divide Basis
b)) (Basis -> T -> (T, T)
cosSin Basis
b T
x)

cot :: Basis -> T -> T
cot :: Basis -> T -> T
cot Basis
b T
x = (T -> T -> T) -> (T, T) -> T
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (Basis -> T -> T -> T
divide Basis
b) (Basis -> T -> (T, T)
cosSin Basis
b T
x)


{- ** logarithmic functions -}

lnSeries :: Basis -> T -> Series
lnSeries :: Basis -> T -> Series
lnSeries Basis
b T
xOrig =
   let x :: T
x   = Basis -> T -> T
negativeExp Basis
b (Basis -> T -> T -> T
sub Basis
b T
xOrig T
one)
       mx :: T
mx  = Basis -> T -> T
neg Basis
b T
x
       xps :: [T]
xps = (Basis -> T -> T) -> Mantissa -> [T] -> [T]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Basis -> Basis -> T -> T
divInt Basis
b) [Basis
1..] ((T -> T) -> T -> [T]
forall a. (a -> a) -> a -> [a]
iterate (Basis -> T -> T -> T
mul Basis
b T
mx) T
x)
   in  (T -> (Basis, T)) -> [T] -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\T
xp -> (T -> Basis
forall a b. (a, b) -> a
fst T
xp, T
xp)) [T]
xps

lnSmall :: Basis -> T -> T
lnSmall :: Basis -> T -> T
lnSmall Basis
b T
x = Basis -> Series -> T
series Basis
b (Basis -> T -> Series
lnSeries Basis
b T
x)

{- |
@
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/1.4140) / 2 == 1.414213578500707
(1.4149 + 2/1.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
lnNewton :: Basis -> T -> T
lnNewton Basis
b T
y =
   let estimate :: T
estimate = Basis -> (Double -> Double) -> T -> T
liftDoubleApprox Basis
b Double -> Double
forall a. C a => a -> a
log T
y
       expRes :: T
expRes   = Basis -> T -> T -> T
mul Basis
b T
y (Basis -> T -> T
expSmall Basis
b (Basis -> T -> T
neg Basis
b T
estimate))
       -- try to reduce dependencies by feeding expSmall with a small argument
       residue :: T
residue =
          Basis -> T -> T -> T
sub Basis
b (Basis -> T -> T -> T
mul Basis
b T
expRes (Basis -> T -> T
expSmallLazy Basis
b (Basis -> T -> T
neg Basis
b T
resTrim))) T
one
       resTrim :: T
resTrim =
          -- (-3, replicate 4 0 ++ alignMant b (-7) residue)
          Basis -> Basis -> T -> T
align Basis
b (- Basis -> Basis
mantLengthDouble Basis
b) T
residue
       lazyAdd :: (Basis, [a]) -> (Basis, [a]) -> (Basis, [a])
lazyAdd (Basis
xe,[a]
xm) (Basis
ye,[a]
ym) =
          (Basis
xe, Basis -> [a] -> [a] -> [a]
forall a. C a => Basis -> [a] -> [a] -> [a]
LPoly.addShifted (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
ye) [a]
xm [a]
ym)
       x :: T
x = T -> T -> T
forall a. C a => (Basis, [a]) -> (Basis, [a]) -> (Basis, [a])
lazyAdd T
estimate T
resTrim
   in  T
x

lnNewton' :: Basis -> T -> T
lnNewton' :: Basis -> T -> T
lnNewton' Basis
b T
y =
   let estimate :: T
estimate = Basis -> (Double -> Double) -> T -> T
liftDoubleApprox Basis
b Double -> Double
forall a. C a => a -> a
log T
y
       residue :: T
residue  =
          Basis -> T -> T -> T
sub Basis
b (Basis -> T -> T -> T
mul Basis
b T
y (Basis -> T -> T
expSmall Basis
b (Basis -> T -> T
neg Basis
b T
x))) T
one
          -- sub b (mul b y (expSmall b (neg b estimate))) one
          -- sub b (mul b y (expSmall b (neg b
          --     (fst estimate, snd estimate ++ [undefined])))) one
       resTrim :: T
resTrim =
          -- align b (-6) residue
          Basis -> Basis -> T -> T
align Basis
b (- Basis -> Basis
mantLengthDouble Basis
b) T
residue
             -- align returns the new exponent immediately
          -- nest (mantLengthDouble b) trimOnce residue
          -- negativeExp b residue
       lazyAdd :: (Basis, [a]) -> (Basis, [a]) -> (Basis, [a])
lazyAdd (Basis
xe,[a]
xm) (Basis
ye,[a]
ym) =
          (Basis
xe, Basis -> [a] -> [a] -> [a]
forall a. C a => Basis -> [a] -> [a] -> [a]
LPoly.addShifted (Basis
xeBasis -> Basis -> Basis
forall a. C a => a -> a -> a
-Basis
ye) [a]
xm [a]
ym)
       x :: T
x = T -> T -> T
forall a. C a => (Basis, [a]) -> (Basis, [a]) -> (Basis, [a])
lazyAdd T
estimate T
resTrim
          -- add b estimate resTrim
                -- LPoly.add checks for empty lists and is thus too strict
   in  T
x


ln :: Basis -> T -> T
ln :: Basis -> T -> T
ln Basis
b x :: T
x@(Basis
xe,Mantissa
_) =
   let e :: Basis
e  = Double -> Basis
forall a b. (C a, C b) => a -> b
round (Double -> Double
forall a. C a => a -> a
log (Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
b) Double -> Double -> Double
forall a. C a => a -> a -> a
* Basis -> Double
forall a b. (C a, C b) => a -> b
fromIntegral Basis
xe :: Double)
       ei :: Integer
ei = Basis -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Basis
e
       y :: T
y  = T -> T
trim (T -> T) -> T -> T
forall a b. (a -> b) -> a -> b
$
          if Basis
eBasis -> Basis -> Bool
forall a. Ord a => a -> a -> Bool
<Basis
0
            then (T -> T -> T) -> T -> T -> Integer -> T
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative (Basis -> T -> T -> T
mul Basis
b) T
x (Basis -> T
eConst Basis
b)    (-Integer
ei)
            else (T -> T -> T) -> T -> T -> Integer -> T
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative (Basis -> T -> T -> T
mul Basis
b) T
x (Basis -> T
recipEConst Basis
b) Integer
ei
       estimate :: T
estimate = Basis -> (Double -> Double) -> T -> T
liftDoubleApprox Basis
b Double -> Double
forall a. C a => a -> a
log T
y
       residue :: T
residue  = Basis -> T -> T -> T
mul Basis
b (Basis -> T -> T
expSmall Basis
b (Basis -> T -> T
neg Basis
b T
estimate)) T
y
   in  Basis -> [T] -> T
addSome Basis
b [(Basis
0,[Basis
e]), T
estimate, Basis -> T -> T
lnSmall Basis
b T
residue]


{- |
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.
-}
angle :: Basis -> (T,T) -> T
angle :: Basis -> (T, T) -> T
angle Basis
b (T
cosx, T
sinx) =
   let wd :: Double
wd      = Double -> Double -> Double
forall a. C a => a -> a -> a
atan2 (Basis -> T -> Double
toDouble Basis
b T
sinx) (Basis -> T -> Double
toDouble Basis
b T
cosx)
       wApprox :: T
wApprox = Basis -> Double -> T
fromDoubleApprox Basis
b Double
wd
       (T
cosApprox, T
sinApprox) = Basis -> T -> (T, T)
cosSin Basis
b T
wApprox
       (T
cosD,T
sinD) =
          (Basis -> T -> T -> T
add Basis
b (Basis -> T -> T -> T
mul Basis
b T
cosx T
cosApprox)
                 (Basis -> T -> T -> T
mul Basis
b T
sinx T
sinApprox),
           Basis -> T -> T -> T
sub Basis
b (Basis -> T -> T -> T
mul Basis
b T
sinx T
cosApprox)
                 (Basis -> T -> T -> T
mul Basis
b T
cosx T
sinApprox))
       sinDSmall :: T
sinDSmall = Basis -> T -> T
negativeExp Basis
b T
sinD
   in  Basis -> T -> T -> T
add Basis
b T
wApprox (Basis -> T -> T
arctanSmall Basis
b (Basis -> T -> T -> T
divide Basis
b T
sinDSmall T
cosD))


{- |
Arcus tangens of arguments with absolute value less than @1 \/ sqrt 3@.
-}
arctanSeries :: Basis -> T -> Series
arctanSeries :: Basis -> T -> Series
arctanSeries Basis
b T
xOrig =
   let x :: T
x   = Basis -> T -> T
negativeExp Basis
b T
xOrig
       mx2 :: T
mx2 = Basis -> T -> T
neg Basis
b (Basis -> T -> T -> T
mul Basis
b T
x T
x)
       xps :: [T]
xps = (Basis -> T -> T) -> Mantissa -> [T] -> [T]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Basis -> Basis -> T -> T
divInt Basis
b) [Basis
1,Basis
3..] ((T -> T) -> T -> [T]
forall a. (a -> a) -> a -> [a]
iterate (Basis -> T -> T -> T
mul Basis
b T
mx2) T
x)
   in  (T -> (Basis, T)) -> [T] -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\T
xp -> (T -> Basis
forall a b. (a, b) -> a
fst T
xp, T
xp)) [T]
xps

arctanSmall :: Basis -> T -> T
arctanSmall :: Basis -> T -> T
arctanSmall Basis
b T
x = Basis -> Series -> T
series Basis
b (Basis -> T -> Series
arctanSeries Basis
b T
x)

{- |
Efficient computation of Arcus tangens of an argument of the form @1\/n@.
-}
arctanStem :: Basis -> Digit -> T
arctanStem :: Basis -> Basis -> T
arctanStem Basis
b Basis
n =
   let x :: T
x = (Basis
0, Basis -> Basis -> Mantissa -> Mantissa
divIntMant Basis
b Basis
n [Basis
1])
       divN2 :: T -> T
divN2 = Basis -> Basis -> T -> T
divInt Basis
b Basis
n (T -> T) -> (T -> T) -> T -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Basis -> Basis -> T -> T
divInt Basis
b (-Basis
n)
       {- this one can cause overflows in piConst too easily
       mn2 = - n*n
       divN2 = divInt b mn2
       -}
       xps :: [T]
xps = (Basis -> T -> T) -> Mantissa -> [T] -> [T]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Basis -> Basis -> T -> T
divInt Basis
b) [Basis
1,Basis
3..] ((T -> T) -> T -> [T]
forall a. (a -> a) -> a -> [a]
iterate (T -> T
trim (T -> T) -> (T -> T) -> T -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T -> T
divN2) T
x)
   in  Basis -> Series -> T
series Basis
b ((T -> (Basis, T)) -> [T] -> Series
forall a b. (a -> b) -> [a] -> [b]
map (\T
xp -> (T -> Basis
forall a b. (a, b) -> a
fst T
xp, T
xp)) [T]
xps)


{- |
This implementation gets the first decimal place for free
by calling the arcus tangens implementation for 'Double's.
-}
arctan :: Basis -> T -> T
arctan :: Basis -> T -> T
arctan Basis
b T
x =
   let estimate :: T
estimate = Basis -> (Double -> Double) -> T -> T
liftDoubleRough Basis
b Double -> Double
forall a. C a => a -> a
atan T
x
       tanEst :: T
tanEst   = Basis -> T -> T
tan Basis
b T
estimate
       residue :: T
residue  = Basis -> T -> T -> T
divide Basis
b (Basis -> T -> T -> T
sub Basis
b T
x T
tanEst) (Basis -> T -> T -> T
add Basis
b T
one (Basis -> T -> T -> T
mul Basis
b T
x T
tanEst))
   in  Basis -> [T] -> T
addSome Basis
b [T
estimate, Basis -> T -> T
arctanSmall Basis
b T
residue]

{- |
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.
-}
arctanClassic :: Basis -> T -> T
arctanClassic :: Basis -> T -> T
arctanClassic Basis
b T
x =
   let absX :: T
absX = T -> T
absolute T
x
       pi2 :: T
pi2  = Basis -> Basis -> T -> T
divInt Basis
b Basis
2 (Basis -> T
piConst Basis
b)
   in  T -> [(Bool, T)] -> T
forall a. a -> [(Bool, a)] -> a
select
          (Basis -> Basis -> T -> T
divInt Basis
b Basis
2 (Basis -> T -> T -> T
sub Basis
b T
pi2
              (Basis -> T -> T
arctanSmall Basis
b
                  (Basis -> Basis -> T -> T
divInt Basis
b Basis
2 (Basis -> T -> T -> T
sub Basis
b (Basis -> T -> T
reciprocal Basis
b T
x) T
x)))))
          [(Basis -> Basis -> T -> T -> Bool
lessApprox Basis
b (-Basis
5) T
absX (Basis -> Rational -> T
fromBaseRational Basis
b (Integer
11Integer -> Integer -> Rational
forall a. C a => a -> a -> T a
%Integer
19)),
               Basis -> T -> T
arctanSmall Basis
b T
x),
           (Basis -> Basis -> T -> T -> Bool
lessApprox Basis
b (-Basis
5) (Basis -> Rational -> T
fromBaseRational Basis
b (Integer
19Integer -> Integer -> Rational
forall a. C a => a -> a -> T a
%Integer
11)) T
absX,
               Basis -> T -> T -> T
sub Basis
b T
pi2 (Basis -> T -> T
arctanSmall Basis
b (Basis -> T -> T
reciprocal Basis
b T
x)))]



{- * constants -}

{- ** elementary -}

zero :: T
zero :: T
zero = (Basis
0,[])

one :: T
one :: T
one = (Basis
0,[Basis
1])

minusOne :: T
minusOne :: T
minusOne = (Basis
0,[-Basis
1])


{- ** transcendental -}

eConst :: Basis -> T
eConst :: Basis -> T
eConst Basis
b = Basis -> T -> T
expSmall Basis
b T
one

recipEConst :: Basis -> T
recipEConst :: Basis -> T
recipEConst Basis
b = Basis -> T -> T
expSmall Basis
b T
minusOne

piConst :: Basis -> T
piConst :: Basis -> T
piConst Basis
b =
   let numCompress :: Mantissa
numCompress = (Basis -> Bool) -> Mantissa -> Mantissa
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Basis
0Basis -> Basis -> Bool
forall a. Eq a => a -> a -> Bool
/=)
          ((Basis -> Basis) -> Basis -> Mantissa
forall a. (a -> a) -> a -> [a]
iterate ((Basis -> Basis -> Basis) -> Basis -> Basis -> Basis
forall a b c. (a -> b -> c) -> b -> a -> c
flip Basis -> Basis -> Basis
forall a. C a => a -> a -> a
div Basis
b) (Basis
4Basis -> Basis -> Basis
forall a. C a => a -> a -> a
*(Basis
44Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
7Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
12Basis -> Basis -> Basis
forall a. C a => a -> a -> a
+Basis
24)))
       stArcTan :: Basis -> Basis -> T
stArcTan Basis
k Basis
den = Basis -> T -> T
scaleSimple Basis
k (Basis -> Basis -> T
arctanStem Basis
b Basis
den)
       sum' :: T
sum' = Basis -> [T] -> T
addSome Basis
b
                 [Basis -> Basis -> T
stArcTan   Basis
44     Basis
57,
                  Basis -> Basis -> T
stArcTan    Basis
7    Basis
239,
                  Basis -> Basis -> T
stArcTan (-Basis
12)   Basis
682,
                  Basis -> Basis -> T
stArcTan   Basis
24  Basis
12943]
   in  (T -> Basis -> T) -> T -> Mantissa -> T
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (T -> Basis -> T
forall a b. a -> b -> a
const (T -> Basis -> T) -> (T -> T) -> T -> Basis -> T
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Basis -> T -> T
compress Basis
b)
             (Basis -> T -> T
scaleSimple Basis
4 T
sum') Mantissa
numCompress



{- * auxilary functions -}

sliceVertPair :: [a] -> [(a,a)]
sliceVertPair :: [a] -> [(a, a)]
sliceVertPair (a
x0:a
x1:[a]
xs) = (a
x0,a
x1) (a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
: [a] -> [(a, a)]
forall a. [a] -> [(a, a)]
sliceVertPair [a]
xs
sliceVertPair [] = []
sliceVertPair [a]
_ = [Char] -> [(a, a)]
forall a. HasCallStack => [Char] -> a
error [Char]
"odd number of elements"



{-
Pi as a zero of trigonometric functions. -
  Is a corresponding computation that bad?
Newton converges quadratically,
  but the involved trigonometric series converge only slightly more than linearly.

-- lift cos to higher frequencies, in order to shift the zero to smaller values, which let trigonometric series converge faster

take 10 $ Numerics.Newton.zero 0.7 (\x -> (cos (2*x), -2 * sin (2*x)))

(\x -> (2 * cos x ^ 2 - 1, -4 * cos x * sin x))
(\x -> (cos x ^ 2 - sin x ^ 2, -4 * cos x * sin x))
(\x -> (tan x ^ 2 - 1, 4 * tan x))


-- compute arctan as inverse of tan by Newton

zero 0.7 (\x -> (tan x - 1, 1 + tan x ^ 2))
zero 0.7 (\x -> (tan x - 1, 1 / cos x ^ 2))
iterate (\x -> x + (cos x - sin x) * cos x) 0.7
iterate (\x -> x + (cos x - sin x) * sqrt 0.5) 0.7
iterate (\x -> x + cos x ^ 2 - sin x * cos x) 0.7
iterate (\x -> x + 0.5 - sin x * cos x) 0.7
iterate (\x -> x + cos x ^ 2 - 0.5) 0.7


-- compute section of tan and cot

zero 0.7 (\x -> (tan x - 1 / tan x, (1 + tan x ^ 2) * (1 + 1 / tan x ^ 2))
zero 0.7 (\x -> ((tan x ^ 2 - 1) * tan x, (1 + tan x ^ 2) ^ 2)
iterate (\x -> x - (sin x ^ 2 - cos x ^ 2) * sin x * cos x) 0.7
iterate (\x -> x - (sin x ^ 2 - cos x ^ 2) * 0.5) 0.7
iterate (\x -> x + 1/2 - sin x ^ 2) 0.7

For using the last formula,
the n-th digit of (sin x) must depend only on the n-th digit of x.
The same holds for (^2).
This means that no interim carry compensation is possible.
This will certainly force usage of Integer for digits,
otherwise the multiplication will overflow sooner or later.
-}