{-# OPTIONS -fno-implicit-prelude #-}
module MathObj.PowerSeries.Example where

import qualified MathObj.PowerSeries as PS

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

import Algebra.Additive (zero, subtract, negate)

import Data.List (map, tail, cycle, zipWith, scanl, intersperse)
import NumericPrelude.List (sieve)

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


{- * Default implementations. -}

recip :: (Ring.C a) => [a]
recip = 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 = expODE
sin = sinODE
cos = cosODE
tan = tanExplSieve
log = logODE
asin = asinODE
acos = acosODE
atan = atanODE

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

pow :: (Field.C a) => a -> [a]
pow = powExpl
sqrt = sqrtExpl


{- * Generate Taylor series explicitly. -}

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

expExpl, sinExpl, cosExpl :: (Field.C a) => [a]
expExpl = scanl (*) one PS.recipProgression
sinExpl = zero : PS.holes2alternate (tail expExpl)
cosExpl =        PS.holes2alternate       expExpl

tanExpl, tanExplSieve :: (ZeroTestable.C a, Field.C a) => [a]
tanExpl = PS.divide sinExpl cosExpl
-- ignore zero values
tanExplSieve =
   concatMap
      (\x -> [zero,x])
      (PS.divide (sieve 2 (tail sin)) (sieve 2 cos))

logExpl, atanExpl, sqrtExpl :: (Field.C a) => [a]
logExpl  = zero : PS.alternate       PS.recipProgression
atanExpl = zero : PS.holes2alternate PS.recipProgression

sinhExpl, coshExpl, atanhExpl :: (Field.C a) => [a]
sinhExpl  = zero : PS.holes2 (tail expExpl)
coshExpl  =        PS.holes2       expExpl
atanhExpl = zero : PS.holes2 PS.recipProgression

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

powExpl :: (Field.C a) => a -> [a]
powExpl expon =
   scanl (*) 1 (zipWith (/)
      (iterate (subtract 1) expon) PS.progression)
sqrtExpl = powExpl (1/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 = PS.integrate 0 $ intersperse 0 $ PS.alternate 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, tanODESieve :: (Field.C a) => [a]
expODE = PS.integrate 1 expODE
sinODE = PS.integrate 0 cosODE
cosODE = PS.integrate 1 (PS.negate sinODE)
tanODE = PS.integrate 0 (PS.add [1] (PS.mul tanODE tanODE))
tanODESieve =
   -- sieve is too strict here because it wants to detect end of lists
   let tan2 = map head (iterate (drop 2) (tail tanODESieve))
   in  PS.integrate 0 (intersperse zero (1 : PS.mul tan2 tan2))

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

logODE, recipCircle, asinODE, atanODE, sqrtODE :: (Field.C a) => [a]
logODE  = PS.integrate zero recip
recipCircle = intersperse zero (PS.alternate (powODE (-1/2)))
asinODE = PS.integrate 0 recipCircle
atanODE = PS.integrate zero (cycle [1,0,-1,0])
sqrtODE = powODE (1/2)

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

sinhODE, coshODE, atanhODE :: (Field.C a) => [a]
sinhODE = PS.integrate 0 coshODE
coshODE = PS.integrate 1 sinhODE
atanhODE = PS.integrate zero (cycle [1,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 expon =
   let y  = PS.integrate 1 y'
       y' = PS.scale expon (scanl1 subtract y)
   in  y