{-# LANGUAGE NoImplicitPrelude #-}moduleMathObj.PowerSeries.CorewhereimportqualifiedMathObj.Polynomial.CoreasPolyimportqualifiedAlgebra.ModuleasModuleimportqualifiedAlgebra.TranscendentalasTranscendentalimportqualifiedAlgebra.FieldasFieldimportqualifiedAlgebra.RingasRingimportqualifiedAlgebra.AdditiveasAdditiveimportqualifiedAlgebra.ZeroTestableasZeroTestableimportqualifiedData.List.MatchasMatchimportqualifiedNumericPrelude.NumericasNPimportqualifiedNumericPrelude.BaseasPimportNumericPrelude.Basehiding(const)importNumericPrelude.Numerichiding(negate,stdUnit,divMod,sqrt,exp,log,sin,cos,tan,asin,acos,atan){-# INLINE evaluate #-}evaluate::Ring.Ca=>[a]->a->aevaluate=flipPoly.horner{-# INLINE evaluateCoeffVector #-}evaluateCoeffVector::Module.Cav=>[v]->a->vevaluateCoeffVector=flipPoly.hornerCoeffVector{-# INLINE evaluateArgVector #-}evaluateArgVector::(Module.Cav,Ring.Cv)=>[a]->v->vevaluateArgVector=flipPoly.hornerArgVector{-# INLINE approximate #-}approximate::Ring.Ca=>[a]->a->[a]approximateyx=scanl(+)zero(zipWith(*)(iterate(x*)1)y){-# INLINE approximateCoeffVector #-}approximateCoeffVector::Module.Cav=>[v]->a->[v]approximateCoeffVectoryx=scanl(+)zero(zipWith(*>)(iterate(x*)1)y){-# INLINE approximateArgVector #-}approximateArgVector::(Module.Cav,Ring.Cv)=>[a]->v->[v]approximateArgVectoryx=scanl(+)zero(zipWith(*>)y(iterate(x*)1)){- * Simple series manipulation -}{- |
For the series of a real function @f@
compute the series for @\x -> f (-x)@
-}alternate::Additive.Ca=>[a]->[a]alternate=zipWithid(cycle[id,NP.negate]){- |
For the series of a real function @f@
compute the series for @\x -> (f x + f (-x)) \/ 2@
-}holes2::Additive.Ca=>[a]->[a]holes2=zipWithid(cycle[id,P.constzero]){- |
For the series of a real function @f@
compute the real series for @\x -> (f (i*x) + f (-i*x)) \/ 2@
-}holes2alternate::Additive.Ca=>[a]->[a]holes2alternate=zipWithid(cycle[id,P.constzero,NP.negate,P.constzero]){- * Series arithmetic -}add,sub::(Additive.Ca)=>[a]->[a]->[a]add=Poly.addsub=Poly.subnegate::(Additive.Ca)=>[a]->[a]negate=Poly.negatescale::Ring.Ca=>a->[a]->[a]scale=Poly.scalemul::Ring.Ca=>[a]->[a]->[a]mul=Poly.mulstripLeadZero::(ZeroTestable.Ca)=>[a]->[a]->([a],[a])stripLeadZero(x:xs)(y:ys)=ifisZerox&&isZeroythenstripLeadZeroxsyselse(x:xs,y:ys)stripLeadZeroxsys=(xs,ys)divMod::(ZeroTestable.Ca,Field.Ca)=>[a]->[a]->([a],[a])divModxsys=let(yZero,yRem)=spanisZeroys(xMod,xRem)=Match.splitAtyZeroxsin(dividexRemyRem,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.Ca)=>[a]->[a]->[a]divide(x:xs)(y:ys)=letzs=map(/y)(x:subxs(mulzsys))inzsdivide[]_=[]divide_[]=error"PowerSeries.divide: division by empty series"{- |
Divide two series also if the divisor has leading zeros.
-}divideStripZero::(ZeroTestable.Ca,Field.Ca)=>[a]->[a]->[a]divideStripZerox'y'=let(x0,y0)=stripLeadZerox'y'inifnully0||isZero(heady0)thenerror"PowerSeries.divideStripZero: Division by zero."elsedividex0y0progression::Ring.Ca=>[a]progression=Poly.progressionrecipProgression::(Field.Ca)=>[a]recipProgression=maprecipprogressiondifferentiate::(Ring.Ca)=>[a]->[a]differentiate=Poly.differentiateintegrate::(Field.Ca)=>a->[a]->[a]integrate=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.
-}sqrt::Field.Ca=>(a->a)->[a]->[a]sqrt_[]=[]sqrtf0(x:xs)=lety=f0xys=map(/(y+y))(xs-(0:mulysys))iny:ys{-
pow alpha t = t^alpha
(pow alpha . x)' = alpha * (pow (alpha-1) . x) * x'
alpha * (pow alpha . x) = x * x' * (pow alpha . x)'
y = pow alpha . x
alpha * y = x * x' * y'
-}{- |
Input series must start with non-zero term.
-}pow::(Field.Ca)=>(a->a)->a->[a]->[a]powf0exponx=lety=integrate(f0(headx))y'y'=scaleexpon(dividey(mulx(differentiatex)))iny{- |
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'
-}exp::Field.Ca=>(a->a)->[a]->[a]expf0x=letx'=differentiatexy=integrate(f0(headx))(mulyx')inysinCos::Field.Ca=>(a->(a,a))->[a]->([a],[a])sinCosf0x=let(y0Sin,y0Cos)=f0(headx)x'=differentiatexySin=integratey0Sin(mulyCosx')yCos=integratey0Cos(negate(mulySinx'))in(ySin,yCos)sinCosScalar::Transcendental.Ca=>a->(a,a)sinCosScalarx=(Transcendental.sinx,Transcendental.cosx)sin,cos::Field.Ca=>(a->(a,a))->[a]->[a]sinf0=fst.sinCosf0cosf0=snd.sinCosf0tan::(Field.Ca)=>(a->(a,a))->[a]->[a]tanf0=uncurrydivide.sinCosf0{-
(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.
-}log::(Field.Ca)=>(a->a)->[a]->[a]logf0x=integrate(f0(headx))(derivedLogx){- |
Computes @(log x)'@, that is @x'\/x@
-}derivedLog::(Field.Ca)=>[a]->[a]derivedLogx=divide(differentiatex)xatan::(Field.Ca)=>(a->a)->[a]->[a]atanf0x=letx'=differentiatexinintegrate(f0(headx))(dividex'([1]+mulxx))asin,acos::(Field.Ca)=>(a->a)->(a->a)->[a]->[a]asinsqrt0f0x=letx'=differentiatexinintegrate(f0(headx))(dividex'(sqrtsqrt0([1]-mulxx)))acos=asin{- |
Since the inner series must start with a zero,
the first term is omitted in y.
-}compose::(Ring.Ca)=>[a]->[a]->[a]composexsy=foldr(\xacc->x:mulyacc)[]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.Ca=>(a->[a])->[a]->[a]composeTaylorx(y:ys)=compose(xy)yscomposeTaylorx[]=x0{-
(x . y) = id
(x' . y) * y' = 1
y' = 1 / (x' . y)
-}{- |
This function returns the series of the function in the form:
(point of the expansion, power series)
This is exceptionally slow and needs cubic run-time.
-}inv::(Field.Ca)=>[a]->(a,[a])invx=lety'=divide[1](compose(differentiatex)(taily))y=integrate0y'-- the first term is zero, which is required for compositionin(headx,y)