-- Copyright (c) David Amos, 2008. All rights reserved.{-# OPTIONS_GHC -fglasgow-exts #-}moduleMath.Algebra.Commutative.MPolywhereimportqualifiedData.MapasMimportData.ListasLimportControl.Arrow(first,second)importData.Ratio(denominator)importMath.Algebra.Field.BaseimportMath.Algebra.Commutative.Monomial-- MULTIVARIATE POLYNOMIALS-- |Type for multivariate polynomials.-- ord is a phantom type defining how terms are ordered, r is the type of the ring we are working over.-- For example, a common choice will be MPoly Grevlex Q, meaning polynomials over Q with the grevlex term orderingnewtypeMPolyordr=MP[(Monomialord,r)]deriving(Eq)-- deriving instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r)-- standalone deriving supported from GHC 6.8instance(Ord(Monomialord),Ordr)=>Ord(MPolyordr)wherecompare(MPts)(MPus)=comparetsusinstance(Showr,Numr)=>Show(MPolyordr)whereshow(MP[])="0"show(MPts)=let(c:cs)=concatMapshowTermtsinifc=='+'thencselsec:cswhereshowTerm(m,c)=caseshowcof"1"->"+"++showm"-1"->"-"++showmcs@(x:_)->(ifx=='-'thencselse'+':cs)++(ifm==1then""elseshowm)instance(Ord(Monomialord),Numr)=>Num(MPolyordr)whereMPts+MPus=MP(mergeTermstsus)negate(MPts)=MP$map(secondnegate)tsMPts*MPus=MP$collect$sortBycmpTerm$[(g*h,c*d)|(g,c)<-ts,(h,d)<-us]{-
-- The following appears to be slightly slower, perhaps because sortBy is compiled
MP (t@(g,c):ts) * MP (u@(h,d):us) =
let MP vs = MP ts * MP us
in MP $ mergeTerms ((g*h,c*d):vs) $ mergeTerms [(g*h,c*d) | (h,d) <- us] [(g*h,c*d) | (g,c) <- ts]
_ * _ = MP []
-}fromInteger0=MP[]fromIntegern=MP[(fromInteger1,fromIntegern)]cmpTerm(a,c)(b,d)=casecompareabofEQ->EQ;GT->LT;LT->GT-- in mpolys we put "larger" terms first-- inputs in descending ordermergeTerms(t@(g,c):ts)(u@(h,d):us)=casecompareghofGT->t:mergeTermsts(u:us)LT->u:mergeTerms(t:ts)usEQ->ife==0thenmergeTermstsuselse(g,e):mergeTermstsuswheree=c+dmergeTermstsus=ts++us-- one of them is nullcollect(t1@(g,c):t2@(h,d):ts)|g==h=collect$(g,c+d):ts|c==0=collect$t2:ts|otherwise=t1:collect(t2:ts)collectts=ts-- Fractional instance so that we can enter fractional coefficients-- Only lets us divide by field elements (with unit monomial), not any other polynomialsinstance(Ord(Monomialord),Fractionalr)=>Fractional(MPolyordr)whererecip(MP[(m,c)])=MP[(recipm,recipc)]-- recip (MP [(m,c)]) | m == fromInteger 1 = MP [(m, recip c)]recip_=error"MPoly.recip: only supported for (non-zero) constants or monomials"-- |Create a variable with the supplied name.-- By convention, variable names should usually be a single letter followed by none, one or two digits.var::String->MPolyGrevlexQvarv=MP[(Monomial$M.singletonv1,1)]::MPolyGrevlexQa,b,c,d,s,t,u,v,w,x,y,z::MPolyGrevlexQa=var"a"b=var"b"c=var"c"d=var"d"s=var"s"t=var"t"u=var"u"v=var"v"w=var"w"x=var"x"y=var"y"z=var"z"x_i=var("x"++showi)x0,x1,x2,x3::MPolyGrevlexQx0=x_0x1=x_1x2=x_2x3=x_3-- convertMP :: Ord (Monomial ord') => MPoly ord k -> MPoly ord' kconvertMP(MPts)=MP$sortBycmpTerm$map(firstconvertM)ts-- |Convert a polynomial to lex term orderingtoLex::MPolyordk->MPolyLexktoLex=convertMP-- |Convert a polynomial to glex term orderingtoGlex::MPolyordk->MPolyGlexktoGlex=convertMP-- |Convert a polynomial to grevlex term orderingtoGrevlex::MPolyordk->MPolyGrevlexktoGrevlex=convertMPtoElim::MPolyordk->MPolyElimktoElim=convertMPvarLexv=toLex$varvvarElimv=toElim$varv-- DIVISION ALGORITHMlt(MP(t:ts))=tlm=fst.ltdeg0=-1deg(MPts)=maximum[degMm|(m,c)<-ts]-- the true degree of the polynomial, not the degree of the leading term-- required for sugar strategy when computing Groebner basismulT(m,c)(m',c')=(m*m',c*c')divT(m,c)(m',c')=(m/m',c/c')dividesT(m,_)(m',_)=dividesMmm'properlyDividesT(m,_)(m',_)=dividesMmm'&&m/=m'lcmT(m,c)(m',c')=(lcmMmm',1)infixl7.*t.*MPts=MP$map(mulTt)ts-- preserves term order-- given f, gs, find as, r such that f = sum (zipWith (*) as gs) + r, with r not divisible by any gquotRemMPfgs=quotRemMP'f(replicaten0,0)wheren=lengthgsquotRemMP'0(us,r)=(us,r)quotRemMP'h(us,r)=divisionSteph(gs,[],us,r)divisionSteph(g:gs,us',u:us,r)=ifltg`dividesT`lththenlett=MP[lth`divT`ltg]h'=h-t*gu'=u+tinquotRemMP'h'(reverseus'++u':us,r)elsedivisionSteph(gs,u:us',us,r)divisionSteph([],us',[],r)=let(lth,h')=splitlthinquotRemMP'h'(reverseus',r+lth)splitlt(MP(t:ts))=(MP[t],MPts)infixl7%%f%%gs=rwhere(_,r)=quotRemMPfgs-- div and mod by single mpolydivModMPfg=(q,r)where([q],r)=quotRemMPf[g]divMPfg=qwhere([q],_)=quotRemMPf[g]modMPfg=rwhere(_,r)=quotRemMPf[g]-- OTHER STUFF-- injection of field elements into polynomial ringinject0=MP[]injectc=MP[(fromInteger1,c)]toMonic0=0toMonic(MPts@((_,c):_))|c==1=MPts|otherwise=MP$map(second(/c))ts-- multiply out all denominatorstoZ(MPts)=MP$map(second(*c))tswherec=fromInteger$foldllcm1$[denominatorc|(m,Qc)<-ts]-- substitute terms for variables in an MPoly-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2substvts(MPus)=sum[injectc*substMm|(m,c)<-us]wheresubstM(Monomialm)=product[substVv^i|(v,i)<-M.toListm]substVv=letv'=MP[(Monomial$M.singletonv1,1)]incaseL.lookupv'vtsofJustt->tNothing->v'-- no substitute, so keep as issupport(MPts)=[MP[(m,1)]|m<-reverse$L.sort$support'ts]wheresupport'ts=foldlL.union[][supportMm|(m,c)<-ts]