-- Copyright (c) David Amos, 2008. All rights reserved.{-# OPTIONS_GHC -fglasgow-exts #-}moduleMath.Algebra.Commutative.GBasiswhereimportData.ListimportqualifiedData.MapasMimportMath.Algebra.Commutative.MonomialimportMath.Algebra.Commutative.MPoly-- Sources:-- Cox, Little, O'Shea: Ideals, Varieties and Algorithms-- Giovini, Mora, Niesi, Robbiano, Traverso, "One sugar cube please, or Selection strategies in the Buchberger algorithm"sPolyfg=leth=lcmT(ltf)(ltg)inh`divT`ltf.*f-h`divT`ltg.*g-- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second termsisGBfs=all(\h->h%%fs==0)(pairWithsPolyfs)-- Cox p87gb1fs=gb'fs(pairWithsPolyfs)wheregb'gs(h:hs)=leth'=h%%gsinifh'==0thengb'gshselsegb'(h':gs)(hs++map(sPolyh')gs)gb'gs[]=gs-- [f xi xj | xi <- xs, xj <- xs, i < j]pairWithf(x:xs)=map(fx)xs++pairWithfxspairWith_[]=[]-- Cox p89-90reducegs=reduce'[]gswherereduce'gs'(g:gs)|g'==0=reduce'gs'gs|otherwise=reduce'(g':gs')gswhereg'=g%%(gs'++gs)reduce'gs'[]=reverse$sort$maptoMonicgs'-- the reverse means that when using an elimination order, the elimination ideal will be at the end-- Cox et al p106-7-- No need to calculate an spoly fi fj if-- 1. the lm fi and lm fj are coprime, or-- 2. there exists some fk, with (i,k) (j,k) already considered, and lm fk divides lcm (lm fi) (lm fj) -- some slight inefficiencies from looking up fi, fj repeatedlygb2fs=reduce$gb'fs(pairs[1..s])swheres=lengthfsgb'gs((i,j):ps)t=letfi=gs!i;fj=gs!jinifcoprimeM(lmfi)(lmfj)||criterionfifj-- if lcmM (lm fi) (lm fj) == lm fi * lm fj || criterion fi fjthengb'gspstelseleth=sPolyfifj%%gsinifh==0thengb'gspstelsegb'(gs++[h])(ps++[(i,t+1)|i<-[1..t]])(t+1)wherecriterionfifj=letl=lcmM(lmfi)(lmfj)inany(testl)[1..t]testlk=k`notElem`[i,j]&&ordpairik`notElem`ps&&ordpairjk`notElem`ps&&lm(gs!k)`dividesM`lgb'gs[]_=gspairs(x:xs)=map(\y->(x,y))xs++pairsxspairs[]=[]xs!i=xs!!(i-1)-- in other words, index the list from 1 not 0ordpairxy|x<y=(x,y)|otherwise=(y,x)-- version of gb2 where we eliminate pairs as they're created, rather than as they're processedgb2bfs=reduce$gb1'[]fs[]0wheregb1'gs(f:fs)pst=gb1'(gs++[f])fsps'(t+1)whereps'=updatePairsgspsftgb1'ls[]pst=gb2'lspstgb2'gs((i,j):ps)t=leth=sPoly(gs!i)(gs!j)%%gsinifh==0thengb2'gspstelseletps'=updatePairsgs((i,j):ps)htingb2'(gs++[h])ps'(t+1)gb2'gs[]_=gsupdatePairsgspsft=[p|p@(i,j)<-ps,not(lmf`dividesM`lcmM(lm(gs!i))(lm(gs!j)))]++[(i,t+1)|(gi,i)<-zipgs[1..t],not(coprimeM(lmgi)(lmf)),not(criterion(lcmM(lmgi)(lmf))i)]wherecriterionli=any(`dividesM`l)[lmgk|(gk,k)<-zipgs[1..t],k/=i,ordpairik`notElem`ps]-- Cox et al 108-- 1. list smallest fs first, as more likely to reduce-- 2. order the pairs with smallest lcm fi fj first ("normal selection strategy")gb3bfs=letfs'=sort$filter(/=0)fsinreduce$gb1'[]fs'[]0wheregb1'gs(f:fs)pst=gb1'(gs++[f])fsps'(t+1)whereps'=updatePairsgspsftgb1'ls[]pst=gb2'lspstgb2'gs((l,(i,j)):ps)t=leth=sPoly(gs!i)(gs!j)%%gsinifh==0thengb2'gspstelseletps'=updatePairsgs((l,(i,j)):ps)htingb2'(gs++[h])ps'(t+1)gb2'gs[]_=gsupdatePairs::(Ord(Monomialord),Fractionalr)=>[MPolyordr]->[(Monomialord,(Int,Int))]->(MPolyordr)->Int->[(Monomialord,(Int,Int))]updatePairsgspsft=mergeBycmpFst[p|p@(l,(i,j))<-ps,not(lmf`dividesM`l)]$sortBycmpFst[(l,(i,t+1))|(gi,i)<-zipgs[1..t],l<-[lcmM(lmgi)(lmf)],not(coprimeM(lmgi)(lmf)),not(criterionli)]wherecriterionli=any(`dividesM`l)[lmgk|(gk,k)<-zipgs[1..t],k/=i,ordpairik`notElem`mapsndps]cmpFst(a,_)(b,_)=compareabmergeBycmp(t:ts)(u:us)=casecmptuofLT->t:mergeBycmpts(u:us)EQ->t:mergeBycmpts(u:us)GT->u:mergeBycmp(t:ts)usmergeBy_tsus=ts++us-- one of them is null-- naive implementation of "sugar strategy"gb4bfs=letfs'=sort$filter(/=0)fsinreduce$gb1'[]fs'[]0wheregb1'gs(f:fs)pst=gb1'(gs++[f])fsps'(t+1)whereps'=updatePairsgspsftgb1'ls[]pst=gb2'lspstgb2'gs((sl,(i,j)):ps)t=leth=sPoly(gs!i)(gs!j)%%gsinifh==0thengb2'gspstelseletps'=updatePairsgs((sl,(i,j)):ps)htingb2'(gs++[h])ps'(t+1)gb2'gs[]_=gsupdatePairs::(Ord(Monomialord),Fractionalr)=>[MPolyordr]->[((Int,Monomialord),(Int,Int))]->(MPolyordr)->Int->[((Int,Monomialord),(Int,Int))]updatePairsgspsft=mergeBycmpFst[p|p@((s,l),(i,j))<-ps,not(lmf`dividesM`l)]$sortBycmpFst[((s,l),(i,t+1))|(gi,i)<-zipgs[1..t],l<-[lcmM(lmgi)(lmf)],s<-[sugargifl],not(coprimeM(lmgi)(lmf)),not(criterionli)]wherecriterionli=any(`dividesM`l)[lmgk|(gk,k)<-zipgs[1..t],k/=i,ordpairik`notElem`mapsndps]-- Giovini et al-- The point of sugar is, given fi, fj, to give an upper bound on the degree of sPoly fi fj without having to calculate it-- We can then select by preference pairs with lower sugar, expecting therefore that the s-polys will have lower degree-- |Given a list of polynomials over a field, return a Groebner basis for the ideal generated by the polynomialsgb::(Ord(Monomialord),Fractionalk,Ordk)=>[MPolyordk]->[MPolyordk]gbfs=-- let fs' = sort $ filter (/=0) fsletfs'=sort$maptoMonic$filter(/=0)fsinreduce$gb1'[]fs'[]0wheregb1'gs(f:fs)pst=gb1'(gs++[f])fsps'(t+1)whereps'=updatePairsgspsf(t+1)gb1'ls[]pst=gb2'lspstgb2'gs(p@(_,(i,j)):ps)t=ifh==0thengb2'gspstelsegb2'(gs++[h])ps'(t+1)whereh=toMonic$sPoly(gs!i)(gs!j)%%gsps'=updatePairsgs(p:ps)h(t+1)gb2'gs[]_=gsupdatePairs::(Ord(Monomialord),Fractionalr)=>[MPolyordr]->[((Int,Monomialord),(Int,Int))]->(MPolyordr)->Int->[((Int,Monomialord),(Int,Int))]updatePairsgspsgkk=letnewps=[letl=lcmM(lmgi)(lmgk)in((sugargigkl,l),(i,k))|(gi,i)<-zipgs[1..k-1]]ps'=[p|p@((sij,tij),(i,j))<-ps,((sik,tik),_)<-[newps!i],((sjk,tjk),_)<-[newps!j],not((tik`properlyDividesM`tij)&&(tjk`properlyDividesM`tij))]-- sloppy variantnewps'=discard1[]newpsnewps''=sortBycmpSug$discard2[]$sortBycmpNormalnewps'inmergeBycmpSugps'newps''wherediscard1ls(r@((_sik,tik),(i,_k)):rs)=iflm(gs!i)`coprimeM`lmgk-- then discard [l | l@((_,tjk),_) <- ls, tjk /= tik] [r | r@((_,tjk),_) <- ls, tjk /= tik]thendiscard1(filter(\((_,tjk),_)->tjk/=tik)ls)(filter(\((_,tjk),_)->tjk/=tik)rs)elsediscard1(r:ls)rsdiscard1ls[]=lsdiscard2ls(r@((_sik,tik),(i,k)):rs)=discard2(r:ls)$filter(\((_sjk,tjk),(j,k'))->not(k==k'&&tik`dividesM`tjk))rsdiscard2ls[]=ls-- The type annotation on updatePairs appears to be required-- The two calls to toMonic are designed to prevent coefficient explosion, but it is unproven that they are effective-- sugar of sPoly f g, where h = lcm (lt f) (lt g)-- this is an upper bound on deg (sPoly f g)sugarfgh=degMh+max(degf-degM(lmf))(degg-degM(lmg))cmpNormal((s1,t1),(i1,j1))((s2,t2),(i2,j2))=compare(t1,j1)(t2,j2)cmpSug((s1,t1),(i1,j1))((s2,t2),(i2,j2))=compare(s1,t1,j1)(s2,t2,j2)-- earlier version of gb3bgb3fs=letgs=sort$filter(/=0)fsps=sortBycmpFst$pairWith(\(i,fi)(j,fj)->(lcmM(lmfi)(lmfj),(i,j)))$zip[1..]gsinreduce$gb'gspsswheres=lengthfsgb'::(Ord(Monomialord),Fractionalr)=>[MPolyordr]->[(Monomialord,(Int,Int))]->Int->[MPolyordr]gb'gs((l,(i,j)):ps)t=letfi=gs!i;fj=gs!jinifcoprimeM(lmfi)(lmfj)||any(testl)[1..t]thengb'gspstelseleth=sPolyfifj%%gsinifh==0thengb'gspstelseletps'=mergeBycmpFstps$sortBycmpFst$zip[lcmM(lmh)(lmfi)|fi<-gs][(i,t+1)|i<-[1..t]]-- else let ps' = mergeBy cmpFst ps $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]]ingb'(gs++[h])ps'(t+1)wheretestlk=k`notElem`[i,j]&&ordpairik`notElem`mapsndps&&ordpairjk`notElem`mapsndps&&lm(gs!k)`dividesM`lgb'gs[]_=gs-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm-- earlier version of gb4bgb4fs=letgs=sort$filter(/=0)fsps=sortBycmpFst$pairWith(\(i,fi)(j,fj)->letl=lcmM(lmfi)(lmfj)in((sugarfifjl,l),(i,j)))$zip[1..]gsinreduce$gb'gspsswheres=lengthfsgb'::(Ord(Monomialord),Fractionalr)=>[MPolyordr]->[((Int,Monomialord),(Int,Int))]->Int->[MPolyordr]gb'gs(((s,l),(i,j)):ps)t=letfi=gs!i;fj=gs!jinifcoprimeM(lmfi)(lmfj)||any(testl)[1..t]thengb'gspstelseleth=sPolyfifj%%gsinifh==0thengb'gspstelseletps'=mergeBycmpFstps$sortBycmpFst$zip[letl=lcmM(lmfi)(lmh)in(sugarfihl,l)|fi<-gs][(i,t+1)|i<-[1..t]]ingb'(gs++[h])ps'(t+1)wheretestlk=k`notElem`[i,j]&&ordpairik`notElem`mapsndps&&ordpairjk`notElem`mapsndps&&lm(gs!k)`dividesM`lgb'gs[]_=gs-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm{-
merge (t:ts) (u:us) =
if t <= u
then t : merge ts (u:us)
else u : merge (t:ts) us
merge ts us = ts ++ us -- one of them is null
-}-- OPERATIONS ON IDEALS-- Cox et al, p181-- Geometric interpretation: V(I+J) = V(I) `intersect` V(J)sumIfsgs=gb$fs++gs-- Cox et al, p183-- Geometric interpretation: V(I.J) = V(I) `union` V(J)productIfsgs=gb[f*g|f<-fs,g<-gs]