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

import qualified MathObj.Polynomial.Core as Poly

import qualified Algebra.Module         as Module
import qualified Algebra.Transcendental as Transcendental
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 Data.List.Match as Match
import qualified NumericPrelude.Numeric as NP
import qualified NumericPrelude.Base as P

import NumericPrelude.Base    hiding (const)
import NumericPrelude.Numeric hiding (negate, stdUnit, divMod,
                              sqrt, exp, log,
                              sin, cos, tan, asin, acos, atan)


{- $setup
>>> import qualified MathObj.PowerSeries.Core as PS
>>> import qualified MathObj.PowerSeries.Example as PSE
>>> import Test.NumericPrelude.Utility (equalTrunc, (/\))
>>> import qualified Test.QuickCheck as QC
>>> import NumericPrelude.Numeric as NP
>>> import NumericPrelude.Base as P
>>> import Prelude ()
>>> import Control.Applicative (liftA3)
>>>
>>> checkHoles ::
>>>    Int -> ([Rational] -> [Rational]) ->
>>>    Rational -> [Rational] -> QC.Property
>>> checkHoles trunc f x xs =
>>>    QC.choose (1,10) /\ \expon ->
>>>    equalTrunc trunc
>>>       (f (PS.insertHoles expon (x:xs)) ++ repeat zero)
>>>       (PS.insertHoles expon (f (x:xs)) ++ repeat zero)
>>>
>>> genInvertible :: QC.Gen [Rational]
>>> genInvertible =
>>>    liftA3 (\x0 x1 xs -> x0:x1:xs)
>>>       QC.arbitrary (fmap QC.getNonZero QC.arbitrary) QC.arbitrary
-}


{-# INLINE evaluate #-}
evaluate :: Ring.C a => [a] -> a -> a
evaluate :: [a] -> a -> a
evaluate = (a -> [a] -> a) -> [a] -> a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> [a] -> a
forall a. C a => a -> [a] -> a
Poly.horner

{-# INLINE evaluateCoeffVector #-}
evaluateCoeffVector :: Module.C a v => [v] -> a -> v
evaluateCoeffVector :: [v] -> a -> v
evaluateCoeffVector = (a -> [v] -> v) -> [v] -> a -> v
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> [v] -> v
forall a v. C a v => a -> [v] -> v
Poly.hornerCoeffVector

{-# INLINE evaluateArgVector #-}
evaluateArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> v
evaluateArgVector :: [a] -> v -> v
evaluateArgVector = (v -> [a] -> v) -> [a] -> v -> v
forall a b c. (a -> b -> c) -> b -> a -> c
flip v -> [a] -> v
forall a v. (C a v, C v) => v -> [a] -> v
Poly.hornerArgVector


{-# INLINE approximate #-}
approximate :: Ring.C a => [a] -> a -> [a]
approximate :: [a] -> a -> [a]
approximate [a]
y a
x =
   (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
zero ((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
xa -> a -> a
forall a. C a => a -> a -> a
*) a
1) [a]
y)

{-# INLINE approximateCoeffVector #-}
approximateCoeffVector :: Module.C a v => [v] -> a -> [v]
approximateCoeffVector :: [v] -> a -> [v]
approximateCoeffVector [v]
y a
x =
   (v -> v -> v) -> v -> [v] -> [v]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl v -> v -> v
forall a. C a => a -> a -> a
(+) v
forall a. C a => a
zero ((a -> v -> v) -> [a] -> [v] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> v -> v
forall a v. C a v => a -> v -> v
(*>) ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
xa -> a -> a
forall a. C a => a -> a -> a
*) a
1) [v]
y)

{-# INLINE approximateArgVector #-}
approximateArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> [v]
approximateArgVector :: [a] -> v -> [v]
approximateArgVector [a]
y v
x =
   (v -> v -> v) -> v -> [v] -> [v]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl v -> v -> v
forall a. C a => a -> a -> a
(+) v
forall a. C a => a
zero ((a -> v -> v) -> [a] -> [v] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> v -> v
forall a v. C a v => a -> v -> v
(*>) [a]
y ((v -> v) -> v -> [v]
forall a. (a -> a) -> a -> [a]
iterate (v
xv -> v -> v
forall a. C a => a -> a -> a
*) v
1))


{- * Simple series manipulation -}

{- |
For the series of a real function @f@
compute the series for @\x -> f (-x)@
-}

alternate :: Additive.C a => [a] -> [a]
alternate :: [a] -> [a]
alternate = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a. a -> a
id ([a -> a] -> [a -> a]
forall a. [a] -> [a]
cycle [a -> a
forall a. a -> a
id, a -> a
forall a. C a => a -> a
NP.negate])

{- |
For the series of a real function @f@
compute the series for @\x -> (f x + f (-x)) \/ 2@
-}

holes2 :: Additive.C a => [a] -> [a]
holes2 :: [a] -> [a]
holes2 = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a. a -> a
id ([a -> a] -> [a -> a]
forall a. [a] -> [a]
cycle [a -> a
forall a. a -> a
id, a -> a -> a
forall a b. a -> b -> a
P.const a
forall a. C a => a
zero])

{- |
For the series of a real function @f@
compute the real series for @\x -> (f (i*x) + f (-i*x)) \/ 2@
-}
holes2alternate :: Additive.C a => [a] -> [a]
holes2alternate :: [a] -> [a]
holes2alternate =
   ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a. a -> a
id ([a -> a] -> [a -> a]
forall a. [a] -> [a]
cycle [a -> a
forall a. a -> a
id, a -> a -> a
forall a b. a -> b -> a
P.const a
forall a. C a => a
zero, a -> a
forall a. C a => a -> a
NP.negate, a -> a -> a
forall a b. a -> b -> a
P.const a
forall a. C a => a
zero])


{- |
For power series of @f x@, compute the power series of @f(x^n)@.

prop> QC.choose (1,10) /\ \m -> QC.choose (1,10) /\ \n xs -> equalTrunc 100 (PS.insertHoles m $ PS.insertHoles n xs) (PS.insertHoles (m*n) xs)
-}
insertHoles :: Additive.C a => Int -> [a] -> [a]
insertHoles :: Int -> [a] -> [a]
insertHoles Int
n =
   if Int
nInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<=Int
0
     then [Char] -> [a] -> [a]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [a] -> [a]) -> [Char] -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [Char]
"insertHoles requires positive exponent, but got " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
n
     else (a -> [a]) -> [a] -> [a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\a
x -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) a
forall a. C a => a
zero)


{- * Series arithmetic -}

add, sub :: (Additive.C a) => [a] -> [a] -> [a]
add :: [a] -> [a] -> [a]
add = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
Poly.add
sub :: [a] -> [a] -> [a]
sub = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
Poly.sub

negate :: (Additive.C a) => [a] -> [a]
negate :: [a] -> [a]
negate = [a] -> [a]
forall a. C a => [a] -> [a]
Poly.negate

scale :: Ring.C a => a -> [a] -> [a]
scale :: a -> [a] -> [a]
scale = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
Poly.scale

mul :: Ring.C a => [a] -> [a] -> [a]
mul :: [a] -> [a] -> [a]
mul = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
Poly.mul


stripLeadZero :: (ZeroTestable.C a) => [a] -> [a] -> ([a],[a])
stripLeadZero :: [a] -> [a] -> ([a], [a])
stripLeadZero (a
x:[a]
xs) (a
y:[a]
ys) =
  if a -> Bool
forall a. C a => a -> Bool
isZero a
x Bool -> Bool -> Bool
&& a -> Bool
forall a. C a => a -> Bool
isZero a
y
    then [a] -> [a] -> ([a], [a])
forall a. C a => [a] -> [a] -> ([a], [a])
stripLeadZero [a]
xs [a]
ys
    else (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs,a
ya -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ys)
stripLeadZero [a]
xs [a]
ys = ([a]
xs,[a]
ys)


divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a],[a])
divMod :: [a] -> [a] -> ([a], [a])
divMod [a]
xs [a]
ys =
   let ([a]
yZero,[a]
yRem) = (a -> Bool) -> [a] -> ([a], [a])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span a -> Bool
forall a. C a => a -> Bool
isZero [a]
ys
       ([a]
xMod, [a]
xRem) = [a] -> [a] -> ([a], [a])
forall b a. [b] -> [a] -> ([a], [a])
Match.splitAt [a]
yZero [a]
xs
   in  ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a]
xRem [a]
yRem, [a]
xMod)

{- |
Divide two series where the absolute term of the divisor is non-zero.
That is, power series with leading non-zero terms are the units
in the ring of power series.

Knuth: Seminumerical algorithms
-}
divide :: (Field.C a) => [a] -> [a] -> [a]
divide :: [a] -> [a] -> [a]
divide (a
x:[a]
xs) (a
y:[a]
ys) =
   let zs :: [a]
zs = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. C a => a -> a -> a
/a
y) (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
sub [a]
xs ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
zs [a]
ys))
   in  [a]
zs
divide [] [a]
_ = []
divide [a]
_ [] = [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"PowerSeries.divide: division by empty series"

{- |
Divide two series also if the divisor has leading zeros.
-}
divideStripZero :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> [a]
divideStripZero :: [a] -> [a] -> [a]
divideStripZero [a]
x' [a]
y' =
   let ([a]
x0,[a]
y0) = [a] -> [a] -> ([a], [a])
forall a. C a => [a] -> [a] -> ([a], [a])
stripLeadZero [a]
x' [a]
y'
   in  if [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
y0 Bool -> Bool -> Bool
|| a -> Bool
forall a. C a => a -> Bool
isZero ([a] -> a
forall a. [a] -> a
head [a]
y0)
         then [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"PowerSeries.divideStripZero: Division by zero."
         else [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a]
x0 [a]
y0


progression :: Ring.C a => [a]
progression :: [a]
progression = [a]
forall a. C a => [a]
Poly.progression

recipProgression :: (Field.C a) => [a]
recipProgression :: [a]
recipProgression = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. C a => a -> a
recip [a]
forall a. C a => [a]
progression

differentiate :: (Ring.C a) => [a] -> [a]
differentiate :: [a] -> [a]
differentiate = [a] -> [a]
forall a. C a => [a] -> [a]
Poly.differentiate

integrate :: (Field.C a) => a -> [a] -> [a]
integrate :: a -> [a] -> [a]
integrate = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
Poly.integrate


{- |
We need to compute the square root only of the first term.
That is, if the first term is rational,
then all terms of the series are rational.

prop> equalTrunc 50 PSE.sqrtExpl (PS.sqrt (\1 -> 1) [1,1])
prop> equalTrunc 500 (1:1:repeat 0) (PS.sqrt (\1 -> 1) (PS.mul [1,1] [1,1]))
prop> checkHoles 50 (PS.sqrt (\1 -> 1)) 1
-}
sqrt :: Field.C a => (a -> a) -> [a] -> [a]
sqrt :: (a -> a) -> [a] -> [a]
sqrt a -> a
_ [] = []
sqrt a -> a
f0 (a
x:[a]
xs) =
   let y :: a
y  = a -> a
f0 a
x
       ys :: [a]
ys = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. C a => a -> a -> a
/(a
ya -> a -> a
forall a. C a => a -> a -> a
+a
y)) ([a]
xs [a] -> [a] -> [a]
forall a. C a => a -> a -> a
- (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
ys [a]
ys))
   in  a
ya -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ys

{-
pow alpha t = t^alpha
(pow alpha . x)' = alpha * (pow (alpha-1) . x) * x'
(pow alpha . x)' * x = alpha * (pow alpha . x) * x'

y = pow alpha . x
y' * x = alpha * y * x'

This yields an implementation that is a fused
exp (alpha * log x)
-}

{- |
Input series must start with a non-zero term,
even better with a positive one.

prop> equalTrunc 100 (PSE.powExpl (-1/3)) (PS.pow (\1 -> 1) (-1/3) [1,1])
prop> equalTrunc 50 (PSE.powExpl (-1/3)) (PS.exp (\0 -> 1) (PS.scale (-1/3) PSE.log))
prop> checkHoles 30 (PS.pow (\1 -> 1) (1/3)) 1
prop> checkHoles 30 (PS.pow (\1 -> 1) (2/5)) 1
-}
pow :: (Field.C a) => (a -> a) -> a -> [a] -> [a]
pow :: (a -> a) -> a -> [a] -> [a]
pow a -> a
f0 a
expon [a]
x =
   let y :: [a]
y  = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate (a -> a
f0 ([a] -> a
forall a. [a] -> a
head [a]
x)) [a]
y'
       y' :: [a]
y' = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
scale a
expon ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
y ([a] -> [a]
forall a. C a => [a] -> [a]
derivedLog [a]
x))
   in  [a]
y


{- |
The first term needs a transcendent computation but the others do not.
That's why we accept a function which computes the first term.

> (exp . x)' =   (exp . x) * x'
> (sin . x)' =   (cos . x) * x'
> (cos . x)' = - (sin . x) * x'

prop> equalTrunc 500 PSE.expExpl (PS.exp (\0 -> 1) [0,1])
prop> equalTrunc 100 (1:1:repeat 0) (PS.exp (\0 -> 1) PSE.log)
prop> checkHoles 30 (PS.exp (\0 -> 1)) 0
-}
exp :: Field.C a => (a -> a) -> [a] -> [a]
exp :: (a -> a) -> [a] -> [a]
exp a -> a
f0 [a]
x =
   let x' :: [a]
x' = [a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x
       y :: [a]
y  = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate (a -> a
f0 ([a] -> a
forall a. [a] -> a
head [a]
x)) ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
y [a]
x')
   in  [a]
y

sinCos :: Field.C a => (a -> (a,a)) -> [a] -> ([a],[a])
sinCos :: (a -> (a, a)) -> [a] -> ([a], [a])
sinCos a -> (a, a)
f0 [a]
x =
   let (a
y0Sin, a
y0Cos) = a -> (a, a)
f0 ([a] -> a
forall a. [a] -> a
head [a]
x)
       x' :: [a]
x'   = [a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x
       ySin :: [a]
ySin = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate a
y0Sin         ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
yCos [a]
x')
       yCos :: [a]
yCos = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate a
y0Cos ([a] -> [a]
forall a. C a => [a] -> [a]
negate ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
ySin [a]
x'))
   in  ([a]
ySin, [a]
yCos)

sinCosScalar :: Transcendental.C a => a -> (a,a)
sinCosScalar :: a -> (a, a)
sinCosScalar a
x = (a -> a
forall a. C a => a -> a
Transcendental.sin a
x, a -> a
forall a. C a => a -> a
Transcendental.cos a
x)

{- |
prop> equalTrunc 500 PSE.sinExpl (PS.sin (\0 -> (0,1)) [0,1])
prop> equalTrunc 50 (0:1:repeat 0) (PS.sin (\0 -> (0,1)) PSE.asin)
prop> checkHoles 20 (PS.sin (\0 -> (0,1))) 0
-}
sin :: Field.C a => (a -> (a,a)) -> [a] -> [a]
sin :: (a -> (a, a)) -> [a] -> [a]
sin a -> (a, a)
f0 = ([a], [a]) -> [a]
forall a b. (a, b) -> a
fst (([a], [a]) -> [a]) -> ([a] -> ([a], [a])) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> (a, a)) -> [a] -> ([a], [a])
forall a. C a => (a -> (a, a)) -> [a] -> ([a], [a])
sinCos a -> (a, a)
f0
{- |
prop> equalTrunc 500 PSE.cosExpl (PS.cos (\0 -> (0,1)) [0,1])
prop> checkHoles 20 (PS.cos (\0 -> (0,1))) 0
-}
cos :: Field.C a => (a -> (a,a)) -> [a] -> [a]
cos :: (a -> (a, a)) -> [a] -> [a]
cos a -> (a, a)
f0 = ([a], [a]) -> [a]
forall a b. (a, b) -> b
snd (([a], [a]) -> [a]) -> ([a] -> ([a], [a])) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> (a, a)) -> [a] -> ([a], [a])
forall a. C a => (a -> (a, a)) -> [a] -> ([a], [a])
sinCos a -> (a, a)
f0

{- |
prop> equalTrunc 50 PSE.tanExpl (PS.tan (\0 -> (0,1)) [0,1])
prop> equalTrunc 50 (0:1:repeat 0) (PS.tan (\0 -> (0,1)) PSE.atan)
prop> checkHoles 20 (PS.tan (\0 -> (0,1))) 0
-}
tan :: (Field.C a) => (a -> (a,a)) -> [a] -> [a]
tan :: (a -> (a, a)) -> [a] -> [a]
tan a -> (a, a)
f0 = ([a] -> [a] -> [a]) -> ([a], [a]) -> [a]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide (([a], [a]) -> [a]) -> ([a] -> ([a], [a])) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> (a, a)) -> [a] -> ([a], [a])
forall a. C a => (a -> (a, a)) -> [a] -> ([a], [a])
sinCos a -> (a, a)
f0

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

{- |
Input series must start with non-zero term.

prop> equalTrunc 500 PSE.logExpl (PS.log (\1 -> 0) [1,1])
prop> equalTrunc 100 (0:1:repeat 0) (PS.log (\1 -> 0) PSE.exp)
prop> checkHoles 30 (PS.log (\1 -> 0)) 1
-}
log :: (Field.C a) => (a -> a) -> [a] -> [a]
log :: (a -> a) -> [a] -> [a]
log a -> a
f0 [a]
x = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate (a -> a
f0 ([a] -> a
forall a. [a] -> a
head [a]
x)) ([a] -> [a]
forall a. C a => [a] -> [a]
derivedLog [a]
x)

{- |
Computes @(log x)'@, that is @x'\/x@
-}
derivedLog :: (Field.C a) => [a] -> [a]
derivedLog :: [a] -> [a]
derivedLog [a]
x = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide ([a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x) [a]
x

{- |
prop> equalTrunc 500 PSE.atan (PS.atan (\0 -> 0) [0,1])
prop> equalTrunc 50 (0:1:repeat 0) (PS.atan (\0 -> 0) PSE.tan)
prop> checkHoles 20 (PS.atan (\0 -> 0)) 0
-}
atan :: (Field.C a) => (a -> a) -> [a] -> [a]
atan :: (a -> a) -> [a] -> [a]
atan a -> a
f0 [a]
x =
   let x' :: [a]
x' = [a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x
   in  a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate (a -> a
f0 ([a] -> a
forall a. [a] -> a
head [a]
x)) ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a]
x' ([a
1] [a] -> [a] -> [a]
forall a. C a => a -> a -> a
+ [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
x [a]
x))

{- |
prop> equalTrunc 100 (0:1:repeat 0) (PS.asin (\1 -> 1) (\0 -> 0) PSE.sin)
prop> equalTrunc 50 PSE.asin (PS.asin (\1 -> 1) (\0 -> 0) [0,1])
prop> checkHoles 30 (PS.asin (\1 -> 1) (\0 -> 0)) 0
-}
asin :: (Field.C a) => (a -> a) -> (a -> a) -> [a] -> [a]
asin :: (a -> a) -> (a -> a) -> [a] -> [a]
asin a -> a
sqrt0 a -> a
f0 [a]
x =
   let x' :: [a]
x' = [a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x
   in  a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate (a -> a
f0 ([a] -> a
forall a. [a] -> a
head [a]
x))
                 ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a]
x' ((a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> [a] -> [a]
sqrt a -> a
sqrt0 ([a
1] [a] -> [a] -> [a]
forall a. C a => a -> a -> a
- [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
x [a]
x)))

{- |
Would be a nice test, but we cannot compute exactly with 'pi':

> equalTrunc 50 PSE.acos (PS.acos (\1 -> 1) (\0 -> pi/2) [0,1])
-}
acos :: (Field.C a) => (a -> a) -> (a -> a) -> [a] -> [a]
acos :: (a -> a) -> (a -> a) -> [a] -> [a]
acos = (a -> a) -> (a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> (a -> a) -> [a] -> [a]
asin

{- |
Since the inner series must start with a zero,
the first term is omitted in y.
-}
compose :: (Ring.C a) => [a] -> [a] -> [a]
compose :: [a] -> [a] -> [a]
compose [a]
xs [a]
y = (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
x [a]
acc -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
mul [a]
y [a]
acc) [] [a]
xs


{- |
Compose two power series where the outer series
can be developed for any expansion point.
To be more precise:
The outer series must be expanded with respect to the leading term
of the inner series.
-}
composeTaylor :: Ring.C a => (a -> [a]) -> [a] -> [a]
composeTaylor :: (a -> [a]) -> [a] -> [a]
composeTaylor a -> [a]
x (a
y:[a]
ys) = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
compose (a -> [a]
x a
y) [a]
ys
composeTaylor a -> [a]
x []     = a -> [a]
x a
0


{-
X(t) = t*x(t)
R(t) = t*r(t)

r(t) = 1 / (x(r(t)*t))
R(t)/t
   = 1 / (x(R(t)))
   = 1 / (X(R(t)) / R(t))
   = 1 / (t / R(t))
-}

{- |
This function returns the series of the inverse function in the form:
(point of the expansion, power series).

That is, say we have the equation:

> y = a + f(x)

where function f is given by a power series with f(0) = 0.
We want to solve for x:

> x = f^-1(y-a)

If you pass the power series of @a+f(x)@ to 'inv',
you get @(a, f^-1)@ as answer, where @f^-1@ is a power series.

The linear term of @f@ (the coefficient of @x@) must be non-zero.

This needs cubic run-time and thus is exceptionally slow.
Computing inverse series for special power series might be faster.

prop> genInvertible /\ \xs -> let (y,ys) = PS.inv xs; (z,zs) = PS.invDiff xs in y==z && equalTrunc 15 ys zs
-}
-- how about NonEmpty.T here?
inv :: (Eq a, Field.C a) => [a] -> (a, [a])
inv :: [a] -> (a, [a])
inv [] = [Char] -> (a, [a])
forall a. HasCallStack => [Char] -> a
error [Char]
"inv: power series must be non-zero"
inv (a
x:[a]
xs) =
   (a
x, let r :: [a]
r = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a
1] ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
compose [a]
xs [a]
r) in a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
r)


{-
(x . y) = id
(x' . y) * y' = 1
y' = 1 / (x' . y)
-}

{-
Like 'inv' but with a slightly cumbersome implementation.
-}
invDiff :: (Field.C a) => [a] -> (a, [a])
invDiff :: [a] -> (a, [a])
invDiff [a]
x =
   let y' :: [a]
y' = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
divide [a
1] ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
compose ([a] -> [a]
forall a. C a => [a] -> [a]
differentiate [a]
x) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
y))
       y :: [a]
y  = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
integrate a
0 [a]
y'
            -- the first term is zero, which is required for composition
   in  ([a] -> a
forall a. [a] -> a
head [a]
x, [a]
y)