{- |
Compute Taylor expansions of some standard functions at given points.
The outputs are essentially the derivatives divided by factorials.
-}
module PowerSeries.Taylor where

import qualified PowerSeries as PowerSeries
import qualified Numerics.Differentiation as Diff
import Combinatorics (binomialSeq, binomialSeqGen)

import qualified Prelude
import Prelude hiding (id, recip, exp, log, sin, cos, atan)


type T a = a -> PowerSeries.T a



compose :: Num a => T a -> PowerSeries.T a -> PowerSeries.T a
compose x (y:ys) = PowerSeries.compose (x y) ys
compose x []     = x 0


id :: Fractional a => T a
id x = [x,1]

recip :: Fractional a => T a
recip x = iterate (/ negate x) $ Prelude.recip x

{- |
The exponent must be non-negative.
-}
cardinalPower :: Num a => Integer -> T a
cardinalPower n x =
   reverse $
   zipWith (*) (map fromInteger $ binomialSeq n) $
   iterate (x*) 1

realPower :: Floating a => a -> T a
realPower expon x =
   zipWith (*) (binomialSeqGen expon) $
   iterate (/x) (x**expon)



exp :: Floating a => T a
exp = PowerSeries.fromDerivatives . Diff.exp

log :: Floating a => T a
log x = PowerSeries.integrate (Prelude.log x) (recip x)




cos :: Floating a => T a
cos = PowerSeries.fromDerivatives . Diff.cos

sin :: Floating a => T a
sin = PowerSeries.fromDerivatives . Diff.sin

tan :: Floating a => T a
tan x = PowerSeries.divide (sin x) (cos x)



recipCircle :: Floating a => T a
recipCircle x =
   PowerSeries.divide [1]
      (PowerSeries.sqRt [1-x^(2::Int),-2*x,-1])

acos :: Floating a => T a
acos x =
   PowerSeries.integrate
      (Prelude.acos x)
      (PowerSeries.neg $ recipCircle x)

asin :: Floating a => T a
asin x = PowerSeries.integrate (Prelude.asin x) (recipCircle x)

atan :: Floating a => T a
atan x =
   PowerSeries.integrate
      (Prelude.atan x)
      (PowerSeries.divide [1] [1+x^(2::Int),2*x,1])