{-# LANGUAGE NoImplicitPrelude #-}{-# LANGUAGE MultiParamTypeClasses #-}{-# LANGUAGE FlexibleInstances #-}{- |
Polynomials and rational functions in a single indeterminate.
Polynomials are represented by a list of coefficients.
All non-zero coefficients are listed, but there may be extra '0's at the end.
Usage:
Say you have the ring of 'Integer' numbers
and you want to add a transcendental element @x@,
that is an element, which does not allow for simplifications.
More precisely, for all positive integer exponents @n@
the power @x^n@ cannot be rewritten as a sum of powers with smaller exponents.
The element @x@ must be represented by the polynomial @[0,1]@.
In principle, you can have more than one transcendental element
by using polynomials whose coefficients are polynomials as well.
However, most algorithms on multi-variate polynomials
prefer a different (sparse) representation,
where the ordering of elements is not so fixed.
If you want division, you need "Number.Ratio"s
of polynomials with coefficients from a "Algebra.Field".
You can also compute with an algebraic element,
that is an element which satisfies an algebraic equation like
@x^3-x-1==0@.
Actually, powers of @x@ with exponents above @3@ can be simplified,
since it holds @x^3==x+1@.
You can perform these computations with "Number.ResidueClass" of polynomials,
where the divisor is the polynomial equation that determines @x@.
If the polynomial is irreducible
(in our case @x^3-x-1@ cannot be written as a non-trivial product)
then the residue classes also allow unrestricted division
(except by zero, of course).
That is, using residue classes of polynomials
you can work with roots of polynomial equations
without representing them by radicals
(powers with fractional exponents).
It is well-known, that roots of polynomials of degree above 4
may not be representable by radicals.
-}moduleMathObj.Polynomial(T,fromCoeffs,coeffs,showsExpressionPrec,const,evaluate,evaluateCoeffVector,evaluateArgVector,compose,equal,add,sub,negate,horner,hornerCoeffVector,hornerArgVector,shift,unShift,mul,scale,divMod,tensorProduct,tensorProductAlt,mulShear,mulShearTranspose,progression,differentiate,integrate,integrateInt,fromRoots,alternate,reverse,)whereimportqualifiedAlgebra.DifferentialasDifferentialimportqualifiedAlgebra.VectorSpaceasVectorSpaceimportqualifiedAlgebra.ModuleasModuleimportqualifiedAlgebra.VectorasVectorimportqualifiedAlgebra.FieldasFieldimportqualifiedAlgebra.PrincipalIdealDomainasPIDimportqualifiedAlgebra.UnitsasUnitsimportqualifiedAlgebra.IntegralDomainasIntegralimportqualifiedAlgebra.RingasRingimportqualifiedAlgebra.AdditiveasAdditiveimportqualifiedAlgebra.ZeroTestableasZeroTestableimportqualifiedAlgebra.IndexableasIndexableimportAlgebra.Module((*>))importAlgebra.ZeroTestable(isZero)importControl.Monad(liftM,)importqualifiedData.ListasListimportNumericPrelude.List(zipWithOverlap,)importData.List.HT(dropWhileRev,shear,shearTranspose,outerProduct,)importTest.QuickCheck(Arbitrary(arbitrary,coarbitrary))importqualifiedPreludeasP98importqualifiedPreludeBaseasPimportqualifiedNumericPreludeasNPimportPreludeBasehiding(const,reverse,)importNumericPreludehiding(divMod,negate,stdUnit,)newtypeTa=Cons{coeffs::[a]}{-# INLINE fromCoeffs #-}fromCoeffs::[a]->TafromCoeffs=lift0{-# INLINE lift0 #-}lift0::[a]->Talift0=Cons{-# INLINE lift1 #-}lift1::([a]->[a])->(Ta->Ta)lift1f(Consx0)=Cons(fx0){-# INLINE lift2 #-}lift2::([a]->[a]->[a])->(Ta->Ta->Ta)lift2f(Consx0)(Consx1)=Cons(fx0x1){-
Functor instance is e.g. useful for showing polynomials in residue rings.
@fmap (ResidueClass.concrete 7) (polynomial [1,4,4::ResidueClass.T Integer] * polynomial [1,5,6])@
-}instanceFunctorTwherefmapf(Consxs)=Cons(mapfxs){-# INLINE plusPrec #-}{-# INLINE appPrec #-}plusPrec,appPrec::IntplusPrec=6appPrec=10instance(Showa)=>Show(Ta)whereshowsPrecp(Consxs)=showParen(p>=appPrec)(showString"Polynomial.fromCoeffs ".showsxs){-# INLINE showsExpressionPrec #-}showsExpressionPrec::(Showa,ZeroTestable.Ca,Additive.Ca)=>Int->String->Ta->String->StringshowsExpressionPrecpvarpoly=ifisZeropolythenshowString"0"elseletterms=filter(not.isZero.fst)(zip(coeffspoly)monomials)monomials=id:showString"*".showStringvar:map(\k->showString"*".showStringvar.showString"^".showsk)[(2::Int)..]showsTermxshowsMon=showsPrec(plusPrec+1)x.showsMoninshowParen(p>plusPrec)(foldl(.)id$List.intersperse(showString" + ")$map(uncurryshowsTerm)terms){- |
Horner's scheme for evaluating a polynomial in a ring.
-}{-# INLINE horner #-}horner::Ring.Ca=>a->[a]->ahornerx=foldr(\cval->c+x*val)zero{- |
Horner's scheme for evaluating a polynomial in a module.
-}{-# INLINE hornerCoeffVector #-}hornerCoeffVector::Module.Cav=>a->[v]->vhornerCoeffVectorx=foldr(\cval->c+x*>val)zero{-# INLINE hornerArgVector #-}hornerArgVector::(Module.Cav,Ring.Cv)=>v->[a]->vhornerArgVectorx=foldr(\cval->c*>one+val*x)zero{-# INLINE evaluate #-}evaluate::Ring.Ca=>Ta->a->aevaluate(Consy)x=hornerxy{- |
Here the coefficients are vectors,
for example the coefficients are real and the coefficents are real vectors.
-}{-# INLINE evaluateCoeffVector #-}evaluateCoeffVector::Module.Cav=>Tv->a->vevaluateCoeffVector(Consy)x=hornerCoeffVectorxy{- |
Here the argument is a vector,
for example the coefficients are complex numbers or square matrices
and the coefficents are reals.
-}{-# INLINE evaluateArgVector #-}evaluateArgVector::(Module.Cav,Ring.Cv)=>Ta->v->vevaluateArgVector(Consy)x=hornerArgVectorxy{- |
'compose' is the functional composition of polynomials.
It fulfills
@ eval x . eval y == eval (compose x y) @
-}-- compose :: Module.C a b => T b -> T a -> T a-- compose (Cons x) y = horner y (map const x){-# INLINE compose #-}compose::(Ring.Ca)=>Ta->Ta->Tacompose(Consx)y=hornery(mapconstx){- |
It's also helpful to put a polynomial in canonical form.
'normalize' strips leading coefficients that are zero.
-}{-# INLINE normalize #-}normalize::(ZeroTestable.Ca)=>[a]->[a]normalize=dropWhileRevisZero{- |
Multiply by the variable, used internally.
-}{-# INLINE shift #-}shift::(Additive.Ca)=>[a]->[a]shift[]=[]shiftl=zero:l{-# INLINE unShift #-}unShift::[a]->[a]unShift[]=[]unShift(_:xs)=xs{-# INLINE const #-}const::a->Taconstx=lift0[x]{-# INLINE equal #-}equal::(Eqa,ZeroTestable.Ca)=>[a]->[a]->Boolequalxy=and(zipWithOverlapisZeroisZero(==)xy)instance(Eqa,ZeroTestable.Ca)=>Eq(Ta)where(Consx)==(Consy)=equalxyinstance(Indexable.Ca,ZeroTestable.Ca)=>Indexable.C(Ta)wherecompare=Indexable.liftComparecoeffsinstance(ZeroTestable.Ca)=>ZeroTestable.C(Ta)whereisZero(Consx)=isZeroxadd,sub::(Additive.Ca)=>[a]->[a]->[a]add=(+)sub=(-){-# INLINE negate #-}negate::(Additive.Ca)=>[a]->[a]negate=mapNP.negateinstance(Additive.Ca)=>Additive.C(Ta)where(+)=lift2add(-)=lift2subzero=lift0[]negate=lift1negate{-# INLINE scale #-}scale::Ring.Ca=>a->[a]->[a]scales=map(s*)instanceVector.CTwherezero=zero(<+>)=(+)(*>)=Vector.functorScaleinstance(Module.Cab)=>Module.Ca(Tb)where(*>)x=lift1(x*>)instance(Field.Ca,Module.Cab)=>VectorSpace.Ca(Tb){-# INLINE tensorProduct #-}tensorProduct::Ring.Ca=>[a]->[a]->[[a]]tensorProduct=outerProduct(*)tensorProductAlt::Ring.Ca=>[a]->[a]->[[a]]tensorProductAltxsys=map(flipscaleys)xs{- |
'mul' is fast if the second argument is a short polynomial,
'MathObj.PowerSeries.**' relies on that fact.
-}{-# INLINE mul #-}mul::Ring.Ca=>[a]->[a]->[a]{- prevent from generation of many zeros
if the first operand is the empty list -}mul[]=P.const[]mulxs=foldr(\yzs->let(v:vs)=scaleyxsinv:addvszs)[]-- this one fails on infinite lists-- mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) []{-# INLINE mulShear #-}mulShear::Ring.Ca=>[a]->[a]->[a]mulShearxsys=mapsum(shear(tensorProductxsys)){-# INLINE mulShearTranspose #-}mulShearTranspose::Ring.Ca=>[a]->[a]->[a]mulShearTransposexsys=mapsum(shearTranspose(tensorProductxsys))instance(Ring.Ca)=>Ring.C(Ta)whereone=constonefromInteger=const.fromInteger(*)=lift2muldivMod::(ZeroTestable.Ca,Field.Ca)=>[a]->[a]->([a],[a])divModxy=let(y0:ys)=dropWhileisZero(List.reversey)auxlxs'=ifl<0then([],xs')elselet(x0:xs)=xs'q0=x0/y0(d',m')=aux(l-1)(subxs(scaleq0ys))in(q0:d',m')(d,m)=aux(lengthx-lengthy)(List.reversex)inifisZeroythenerror"MathObj.Polynomial: division by zero"else(List.reversed,List.reversem)instance(ZeroTestable.Ca,Field.Ca)=>Integral.C(Ta)wheredivMod(Consx)(Consy)=let(d,m)=divModxyin(Consd,Consm){-# INLINE stdUnit #-}stdUnit::(ZeroTestable.Ca,Ring.Ca)=>[a]->astdUnitx=casenormalizexof[]->onel->lastlinstance(ZeroTestable.Ca,Field.Ca)=>Units.C(Ta)whereisUnit(Cons[])=FalseisUnit(Cons(x0:xs))=not(isZerox0)&&allisZeroxsstdUnit(Consx)=const(stdUnitx)stdUnitInv(Consx)=const(recip(stdUnitx)){-
Polynomials are a Euclidean domain, so no instance is necessary
(although it might be faster).
-}instance(ZeroTestable.Ca,Field.Ca)=>PID.C(Ta){-# INLINE progression #-}progression::Ring.Ca=>[a]progression=iterate(one+)one{-# INLINE differentiate #-}differentiate::(Ring.Ca)=>[a]->[a]differentiate=zipWith(*)progression.tail{-# INLINE integrate #-}integrate::(Field.Ca)=>a->[a]->[a]integratecx=c:zipWith(/)xprogression{- |
Integrates if it is possible to represent the integrated polynomial
in the given ring.
Otherwise undefined coefficients occur.
-}{-# INLINE integrateInt #-}integrateInt::(ZeroTestable.Ca,Integral.Ca)=>a->[a]->[a]integrateIntcx=c:zipWithIntegral.safeDivxprogressioninstance(Ring.Ca)=>Differential.C(Ta)wheredifferentiate=lift1differentiate{-# INLINE fromRoots #-}fromRoots::(Ring.Ca)=>[a]->TafromRoots=Cons.foldl(flipmulLinearFactor)[1]{-# INLINE mulLinearFactor #-}mulLinearFactor::Ring.Ca=>a->[a]->[a]mulLinearFactorxyt@(y:ys)=Additive.negate(x*y):yt-scalexysmulLinearFactor_[]=[]{-# INLINE alternate #-}alternate::Additive.Ca=>[a]->[a]alternate=zipWith($)(cycle[id,Additive.negate]){-# INLINE reverse #-}reverse::Additive.Ca=>Ta->Tareverse=lift1alternate{-
see htam: Wavelet/DyadicResultant
resultant :: Ring.C a => [a] -> [a] -> [a]
resultant xs ys =
discriminant :: Ring.C a => [a] -> a
discriminant xs =
let degree = genericLength xs
in parityFlip (safeDiv (degree*(degree-1)) 2)
(resultant xs (differentiate xs))
`safeDiv` last xs
-}instance(Arbitrarya,ZeroTestable.Ca)=>Arbitrary(Ta)wherearbitrary=liftM(fromCoeffs.normalize)arbitrarycoarbitrary=undefined{- * legacy instances -}{- |
It is disputable whether polynomials shall be represented by number literals or not.
An advantage is, that one can write
let x = polynomial [0,1]
in (x^2+x+1)*(x-1)
However the output looks much different.
-}{-# INLINE legacyInstance #-}legacyInstance::alegacyInstance=error"legacy Ring.C instance for simple input of numeric literals"instance(Ring.Ca,Eqa,Showa,ZeroTestable.Ca)=>P98.Num(Ta)wherefromInteger=const.fromIntegernegate=Additive.negate-- for unary minus(+)=legacyInstance(*)=legacyInstanceabs=legacyInstancesignum=legacyInstanceinstance(Field.Ca,Eqa,Showa,ZeroTestable.Ca)=>P98.Fractional(Ta)wherefromRational=const.fromRational(/)=legacyInstance