{-# LANGUAGE RebindableSyntax #-}
module MathObj.PowerSeries.Example where

import qualified MathObj.PowerSeries.Core as PS

import qualified Algebra.Field          as Field
import qualified Algebra.Ring           as Ring
import qualified Algebra.ZeroTestable   as ZeroTestable
import qualified Algebra.Transcendental as Transcendental

import Algebra.Additive (zero, subtract, negate)

import Data.List (intersperse, )
import Data.List.HT (sieve, )

import NumericPrelude.Numeric (one, (*), (/),
                       fromInteger, {-fromRational,-} pi)
import NumericPrelude.Base -- (Bool, const, map, zipWith, id, (&&), (==))


{- $setup
>>> import qualified MathObj.PowerSeries.Core as PS
>>> import qualified MathObj.PowerSeries.Example as PSE
>>> import Test.NumericPrelude.Utility (equalTrunc)
>>> import NumericPrelude.Numeric as NP
>>> import NumericPrelude.Base as P
>>> import Prelude ()
-}


{- * Default implementations. -}

recip :: (Ring.C a) => [a]
recip :: [a]
recip = [a]
forall a. C a => [a]
recipExpl

exp, sin, cos,
  log, asin, atan, sqrt :: (Field.C a) => [a]
acos :: (Transcendental.C a) => [a]
tan :: (ZeroTestable.C a, Field.C a) => [a]
exp :: [a]
exp = [a]
forall a. C a => [a]
expODE
sin :: [a]
sin = [a]
forall a. C a => [a]
sinODE
cos :: [a]
cos = [a]
forall a. C a => [a]
cosODE
tan :: [a]
tan = [a]
forall a. (C a, C a) => [a]
tanExplSieve
log :: [a]
log = [a]
forall a. C a => [a]
logODE
asin :: [a]
asin = [a]
forall a. C a => [a]
asinODE
acos :: [a]
acos = [a]
forall a. C a => [a]
acosODE
atan :: [a]
atan = [a]
forall a. C a => [a]
atanODE

sinh, cosh, atanh :: (Field.C a) => [a]
sinh :: [a]
sinh  = [a]
forall a. C a => [a]
sinhODE
cosh :: [a]
cosh  = [a]
forall a. C a => [a]
coshODE
atanh :: [a]
atanh = [a]
forall a. C a => [a]
atanhODE


-- | prop> \m n -> equalTrunc 30 (PS.mul (PSE.pow m) (PSE.pow n)) (PSE.pow (m+n))
pow :: (Field.C a) => a -> [a]
pow :: a -> [a]
pow = a -> [a]
forall a. C a => a -> [a]
powExpl
sqrt :: [a]
sqrt = [a]
forall a. C a => [a]
sqrtExpl


{- * Generate Taylor series explicitly. -}

recipExpl :: (Ring.C a) => [a]
recipExpl :: [a]
recipExpl = [a] -> [a]
forall a. [a] -> [a]
cycle [a
1,-a
1]

-- | prop> equalTrunc 500 PSE.expExpl PSE.expODE
expExpl :: (Field.C a) => [a]
expExpl :: [a]
expExpl = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. C a => a -> a -> a
(*) a
forall a. C a => a
one [a]
forall a. C a => [a]
PS.recipProgression
-- | prop> equalTrunc 500 PSE.sinExpl PSE.sinODE
sinExpl :: (Field.C a) => [a]
sinExpl :: [a]
sinExpl = a
forall a. C a => a
zero a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2alternate ([a] -> [a]
forall a. [a] -> [a]
tail [a]
forall a. C a => [a]
expExpl)
-- | prop> equalTrunc 500 PSE.cosExpl PSE.cosODE
cosExpl :: (Field.C a) => [a]
cosExpl :: [a]
cosExpl = [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2alternate [a]
forall a. C a => [a]
expExpl

-- | prop> equalTrunc 50 PSE.tanExpl PSE.tanODE
tanExpl :: (ZeroTestable.C a, Field.C a) => [a]
tanExpl :: [a]
tanExpl = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PS.divide [a]
forall a. C a => [a]
sinExpl [a]
forall a. C a => [a]
cosExpl
-- ignore zero values
-- | prop> equalTrunc 50 PSE.tanExpl PSE.tanExplSieve
tanExplSieve :: (ZeroTestable.C a, Field.C a) => [a]
tanExplSieve :: [a]
tanExplSieve =
   (a -> [a]) -> [a] -> [a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap
      (\a
x -> [a
forall a. C a => a
zero,a
x])
      ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PS.divide (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
sieve Int
2 ([a] -> [a]
forall a. [a] -> [a]
tail [a]
forall a. C a => [a]
sin)) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
sieve Int
2 [a]
forall a. C a => [a]
cos))

-- | prop> equalTrunc 500 PSE.logExpl PSE.logODE
logExpl :: (Field.C a) => [a]
logExpl :: [a]
logExpl  = a
forall a. C a => a
zero a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. C a => [a] -> [a]
PS.alternate       [a]
forall a. C a => [a]
PS.recipProgression
-- | prop> equalTrunc 500 PSE.atanExpl PSE.atanODE
atanExpl :: (Field.C a) => [a]
atanExpl :: [a]
atanExpl = a
forall a. C a => a
zero a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2alternate [a]
forall a. C a => [a]
PS.recipProgression

-- | prop> equalTrunc 500 PSE.sinhExpl PSE.sinhODE
sinhExpl :: (Field.C a) => [a]
sinhExpl :: [a]
sinhExpl  = a
forall a. C a => a
zero a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2 ([a] -> [a]
forall a. [a] -> [a]
tail [a]
forall a. C a => [a]
expExpl)
-- | prop> equalTrunc 500 PSE.coshExpl PSE.coshODE
coshExpl :: (Field.C a) => [a]
coshExpl :: [a]
coshExpl  =        [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2       [a]
forall a. C a => [a]
expExpl
-- | prop> equalTrunc 500 PSE.atanhExpl PSE.atanhODE
atanhExpl :: (Field.C a) => [a]
atanhExpl :: [a]
atanhExpl = a
forall a. C a => a
zero a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. C a => [a] -> [a]
PS.holes2 [a]
forall a. C a => [a]
PS.recipProgression

{- * Power series of (1+x)^expon using the binomial series. -}

-- | prop> \expon -> equalTrunc 50 (PSE.powODE expon) (PSE.powExpl expon)
powExpl :: (Field.C a) => a -> [a]
powExpl :: a -> [a]
powExpl a
expon =
   (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. C a => a -> a -> a
(*) a
1 ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(/)
      ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. C a => a -> a -> a
subtract a
1) a
expon) [a]
forall a. C a => [a]
PS.progression)

-- | prop> equalTrunc 100 PSE.sqrtExpl PSE.sqrtODE
sqrtExpl :: (Field.C a) => [a]
sqrtExpl :: [a]
sqrtExpl = a -> [a]
forall a. C a => a -> [a]
powExpl (a
1a -> a -> a
forall a. C a => a -> a -> a
/a
2)

{- |
Power series of error function (almost).
More precisely @ erf = 2 \/ sqrt pi * integrate (\x -> exp (-x^2)) @,
with @erf 0 = 0@.
-}

erf :: (Field.C a) => [a]
erf :: [a]
erf = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ a -> [a] -> [a]
forall a. a -> [a] -> [a]
intersperse a
0 ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. C a => [a] -> [a]
PS.alternate [a]
forall a. C a => [a]
exp

{-
integrate (\x -> exp (-x^2/2)) :

erf = PS.integrate 0 $ intersperse 0 $
    snd $ mapAccumL (\twoPow c -> (twoPow/(-2), twoPow*c)) 1 exp
-}


{- * Generate Taylor series from differential equations. -}

{-
exp' x == exp x
sin' x == cos x
cos' x == - sin x

tan' x == 1 + tan x ^ 2
       == cos x ^ (-2)
-}

expODE, sinODE, cosODE, tanODE :: (Field.C a) => [a]
expODE :: [a]
expODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
1 [a]
forall a. C a => [a]
expODE
sinODE :: [a]
sinODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 [a]
forall a. C a => [a]
cosODE
cosODE :: [a]
cosODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
1 ([a] -> [a]
forall a. C a => [a] -> [a]
PS.negate [a]
forall a. C a => [a]
sinODE)
tanODE :: [a]
tanODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PS.add [a
1] ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PS.mul [a]
forall a. C a => [a]
tanODE [a]
forall a. C a => [a]
tanODE))
-- | prop> equalTrunc 50 PSE.tanODE PSE.tanODESieve
tanODESieve :: (Field.C a) => [a]
tanODESieve :: [a]
tanODESieve =
   -- sieve is too strict here because it wants to detect end of lists
   let tan2 :: [a]
tan2 = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> a
forall a. [a] -> a
head (([a] -> [a]) -> [a] -> [[a]]
forall a. (a -> a) -> a -> [a]
iterate (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
2) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
forall a. C a => [a]
tanODESieve))
   in  a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 (a -> [a] -> [a]
forall a. a -> [a] -> [a]
intersperse a
forall a. C a => a
zero (a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PS.mul [a]
tan2 [a]
tan2))

{-
log' (1+x) == 1/(1+x)
asin' x == acos' x == 1/sqrt(1-x^2)
atan' x == 1/(1+x^2)
-}

logODE, recipCircle, atanODE, sqrtODE :: (Field.C a) => [a]
logODE :: [a]
logODE  = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
forall a. C a => a
zero [a]
forall a. C a => [a]
recip
recipCircle :: [a]
recipCircle = a -> [a] -> [a]
forall a. a -> [a] -> [a]
intersperse a
forall a. C a => a
zero ([a] -> [a]
forall a. C a => [a] -> [a]
PS.alternate (a -> [a]
forall a. C a => a -> [a]
powODE (-a
1a -> a -> a
forall a. C a => a -> a -> a
/a
2)))
-- | prop> equalTrunc 50 PSE.asinODE (snd $ PS.inv PSE.sinODE)
asinODE :: (Field.C a) => [a]
asinODE :: [a]
asinODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 [a]
forall a. C a => [a]
recipCircle
atanODE :: [a]
atanODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
forall a. C a => a
zero ([a] -> [a]
forall a. [a] -> [a]
cycle [a
1,a
0,-a
1,a
0])
sqrtODE :: [a]
sqrtODE = a -> [a]
forall a. C a => a -> [a]
powODE (a
1a -> a -> a
forall a. C a => a -> a -> a
/a
2)

acosODE :: (Transcendental.C a) => [a]
acosODE :: [a]
acosODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate (a
forall a. C a => a
pia -> a -> a
forall a. C a => a -> a -> a
/a
2) [a]
forall a. C a => [a]
recipCircle

sinhODE, coshODE, atanhODE :: (Field.C a) => [a]
sinhODE :: [a]
sinhODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
0 [a]
forall a. C a => [a]
coshODE
coshODE :: [a]
coshODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
1 [a]
forall a. C a => [a]
sinhODE
atanhODE :: [a]
atanhODE = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
forall a. C a => a
zero ([a] -> [a]
forall a. [a] -> [a]
cycle [a
1,a
0])


{-
Power series for y with
   y x = (1+x) ** alpha
by solving the differential equation
   alpha * y x = (1+x) * y' x
-}

powODE :: (Field.C a) => a -> [a]
powODE :: a -> [a]
powODE a
expon =
   let y :: [a]
y  = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.integrate a
1 [a]
y'
       y' :: [a]
y' = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PS.scale a
expon ((a -> a -> a) -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 a -> a -> a
forall a. C a => a -> a -> a
subtract [a]
y)
   in  [a]
y