Faug\303\250re's F5 algorithm Implementation based on Jean-Charles Faug\303\250re's 2000 paper, A new efficient algorithm for computing Gr\303\266bner bases without reduction to zero F5 and on Till Stegers' 2006 thesis, F5 Revisited. Copyright John Perry, 2007; developed at the University of Southern Mississippi. Be advised that this is SLOW SLOW SLOW. However, the implementation appears to work perfectly for regular sequences. The problem may lie in my choice of inappropriate data structures; I am working on this. Some documentation & explanation provided, but read it carefully. Caveat lector. Some of the commentary may be wrong: writing this helped me understand F5 better, but I don't quite understand how everything works. The implementation has more reductions to zero than Faug\303\250re's for nonregular sequences, but the results are similar to Stegers'. I'm not sure why. The "benchmark examples" below helped squash some very subtle bugs. What works correctly is due to Faug\303\250re and Stegers; anything that's broken is my fault. The reader may also be interested in a version of F5 that uses the Extended First Criterion, which will be available from http://www.math.usm.edu/perry/research.html Real Soon Now (TM). It works on the systems I've tested, but needs a lot more work. Comments are welcome and earnestly desired. Contact me at john dot perry at delete-this-hyphenated-phrase usm dot edu. Released under the GNU public license. See below. restart:

f5_computation := module() export # procedures Crit_Pair, SPols, Reduction, Top_Reduction, Find_Reductor, Top_Reducible, Find_Rewriting, Is_Rewritable, Basis, Partial_Basis, Report_Number_Of_SPolys, Report_Number_Of_Zero_Reductions, Report_Pairs_Giving_Zero_Reductions, Report_Set_Of_Critical_Pairs, Report_Rules, Report_Rule_Counts, Report_Labeled_Polys; local # procedures Add_Rule, Add_Labeled_Poly, Poly, Sig, Sig_Prod, Sig_Comp, First_Sig_Min, Rule_Element, ModuleLoad, # data lps, Rules, last_lp, last_Rules, ordering, characteristic, critical_pairs_used, checked_pairs, number_of_spolys, cps_of_zero_reductions, warn_on_zero_reduction: option package: # DESCRIPTION OF MODULE VARIABLES # lps is the "result" of labeled polynomials with record-keeping; # eventually it contains a Gr\303\266bner basis. # Rules keeps track of which original polynomial generates which # new polynomial. # last_lp tracks the number of labeled polynomials in lps # last_Rules tracks the number of Rules # ordering is the term ordering, made visible to the entire module # characteristic is the characteristic of the ground field, determined # from the ordering # critical_pairs_used is the set of critical pairs computed # (but not necessarily of S-polynomials reduced) # number_of_spolys is the number of S-polynomials generated # cps_of_zero_reductions is a set of critical pairs that reduces to # zero (if the polynomials are a regular sequence, this set # should be empty) # warn_on_zero_reduction is a boolean: do you want to be warned when # an S-polynomial reduces to zero? This should not happen for # regular sequences, but does happen for non-regular sequences (see # examples below). ModuleLoad:=proc() printf("F5, implemented in Maple.\134n"); printf("Based on Faugere's paper and Till Stegers' thesis `F5 Revisited'.\134n"); printf("Copyright 2007 by John Perry\134n"); printf("Developed at the University of Southern Mississippi\134n"); printf("Released under the GNU General Public License, available at\134n"); printf("http://www.gnu.org/licenses/gpl.html\134n"); printf("This program is free software; you can redistribute it and/or modify\134n"); printf("it under the terms of the GNU General Public License as published by\134n"); printf("the Free Software Foundation; either version 2 of the License, or\134n"); printf("(at your option) any later version.\134n\134n"); printf("This program is distributed in the hope that it will be useful,\134n"); printf("but WITHOUT ANY WARRANTY; without even the implied warranty of\134n"); printf("MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\134n"); printf("See the GNU General Public License for more details.\134n\134n"); printf("You should have received a copy of the GNU General Public License\134n"); printf("along with this program; if not, write to the Free Software\134n"); printf("Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA\134n"); end: # ModuleLoad Basis is the "outer loop" of the algorithm. It proceeds incrementally through the polynomials, computing first the Gr\303\266bner basis of , then the Gr\303\266bner basis of , and so on until it computes the Gr\303\266bner basis of F_in. # Main procedure: call this with a list of polynomials # and a term ordering. The output is a Gr\303\266bner basis of the inputs. Basis:=proc(F_in::list,tord) # DESCRIPTION OF LOCAL VARIABLES # m is the number of input polynomials # G is a list of indices of lps that corresponds to the # Gr\303\266bner bases generated incrementally # (so G[i] is a subset of G[i+1]) # i,j,k are counters # sorter is a local procedure that sorts F by decreasing lt # F is the sorted form of F_in local m,G,i,j,k,sorter,F: warn_on_zero_reduction := false: # obtain the characteristic characteristic := tord["algebra"]["commutation"]["ground_ring"] ["characteristic"]; if not type(characteristic,posint) then characteristic := 0: end if: # Sort F by decreasing lt, to start with the "smallest" ones. ordering:=tord: sorter := proc(f,g) local tf,tg: tf:=Groebner[LeadingMonomial](f,ordering): tg:=Groebner[LeadingMonomial](g,ordering): if degree(tf)<>degree(tg) then return evalb(degree(tf)>degree(tg)): else return not Groebner[TestOrder](tf,tg,ordering): end if: end: # sorter F:=sort (F_in,sorter): # Initialize some diagnostics number_of_spolys := 0: cps_of_zero_reductions := {}: critical_pairs_used := {}: checked_pairs := {}: # Initialize the rewrite rules and labeled polynomials, # and the lists of GBs m:=nops(F): Rules := [seq(table(),seqctr=1..m)]: last_Rules := [seq(0,seqctr=1..m)]: lps := table(): G:=[seq([],seqctr=1..m)]: # The first Groebner basis is {F[m]} if characteristic=0 then lps[m]:=[[1,m], 1/Groebner[LeadingCoefficient](F[m],tord)*F[m]]: else lps[m]:=[[1,m], (1/Groebner[LeadingCoefficient](F[m],tord) mod characteristic)*F[m] mod characteristic]: end if: last_lp:=m: G[m]:={m}: # Incrementally compute Groebner bases using Partial_Basis i := m-1: while i > 0 do G[i] := Partial_Basis(F[i],i,G[i+1]): for k in G[i] do if Poly(k)=1 then return [1]: end if: end do: # k i := i - 1: end do: # i return [seq(Poly(seqctr),seqctr in G[1])]: end: # Basis Partial_Basis organizes the essential structure of Faug\303\250re's algorithm. It computes a Gr\303\266bner basis of the input set. After creating a list of essential critical pairs, it loops through this set. On each loop, it removes the pairs of smallest degree, computes the corresponding S-polynomials, reduces them modulo the previous Gr\303\266bner basis, then adds new critical pairs if necessary. # f is a polynomial # i is the current iteration of the outer algorithm # G is a list of lists of indices which indicate previous bases Partial_Basis:=proc(f,i,G_prev) # DESCRIPTION OF LOCAL VARIABLES # G_curr a list of generators of the current ideal; # when the algorithm completes, G_curr is a Groebner basis # P is the list of critical pairs considered in this phase # P_d is the list of critical pairs of degree d, where # d is the smallest degree of a critical pair in P # p is a critical pair in P # S_d is the list of S-polynomials to be considered # R_d is the list of indices of new polynomials in the basis # poly_index is an index used to loop through R_d # new_cps is a list of new critical pairs to add to P # after Reduction # k is a counter # sort_by_increasing_lcm is a local procedure local G_curr,P,P_d,d,p,S_d,R_d,poly_index,new_cps,k, sort_by_increasing_lcm; # add ith polynomial to basis (should be from original set) if characteristic > 0 then lps[i]:=[[1,i], (1/Groebner[LeadingCoefficient](f,ordering) mod characteristic) * f mod characteristic]: else lps[i]:=[[1,i],1/Groebner[LeadingCoefficient](f,ordering)*f]: end if: # initialize record keeping G_curr:=G_prev union {i}: # generate critical pairs from ith polynomial P := `union`(seq(Crit_Pair(i,seqi,G_prev),seqi in G_prev)): # determine which S-polynomials reduce to zero; add the ones that # don't to the basis, generate new critical pairs for those new # polynomials, repeat until finished while P<>{} do # pick the critical pairs of lowest degree, put them in P_d P_d := {P[1]}: d := degree(P[1][1]): for p in P do if degree(p[1])<d then P_d := {p}: d := degree(p[1]): elif degree(p[1])=d then P_d := P_d union {p}: end if: end do: P := P minus P_d: # determine which S-polynomials corresponding to the critical # pairs in P_d need reduction, and reduce them S_d := SPols(P_d): R_d := Reduction(S_d,G_prev,G_curr): critical_pairs_used := critical_pairs_used union P_d: number_of_spolys := number_of_spolys+nops(S_d): # did we add new polynomials to the set? if so, we must consider # the new critical pairs generated for poly_index in R_d do new_cps:=`union`( seq(Crit_Pair(poly_index,g_ctr,G_prev),g_ctr in G_curr) ): P := P union new_cps: G_curr := G_curr union {poly_index}; end do: end do: return G_curr: end: # Partial_Basis Crit_Pair uses Buchberger's criteria and Faug\303\250re's record keeping to determine whether the critical pair for polynomials LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2I1EhRictSSVtc3ViR0YkNiUtRiw2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYsNiVRImtGJ0Y1RjgvRjlRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSIsRidGQC8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0ZLLyUqc3ltbWV0cmljR0ZLLyUobGFyZ2VvcEdGSy8lLm1vdmFibGVsaW1pdHNHRksvJSdhY2NlbnRHRksvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GMDYlRjItRiM2JC1GLDYlUSJsRidGNUY4RkBGQkYrRkA= needs to be computed. # k,l are the indices of the critical pair in question # G_prev is the previous Groebner basis Crit_Pair:=proc(k,l,G_prev) # DESCRIPTION OF LOCAL VARIABLES # tk,tl are the leading terms of the polynomials labeled by lps[k],lab_poly[l] # t is the lcm of tk,tl # u1,u2 are the multiples necessary to obtain the S-polynomial # of lps[k],lps[l] # kcp,lcp are the indices of lps[k],lps[l] (maybe not in that order) # Sigk1,Sigk2 are the signatures of lps[kcp],lps[lcp] # k1,k2 correspond to the generators of lps[kcp],lps[lcp], resp. # kprime is a temporary variable for swapping kcp,lcp if needed local tk,tl,t,u1,u2,Sigk1,Sigk2,k1,k2,kcp,lcp,d,g: # initialize and get polys in desired order tk:=Groebner[LeadingMonomial](Poly(k),ordering): tl:=Groebner[LeadingMonomial](Poly(l),ordering): t:=lcm(tk,tl); u1:=t/tk: u2:=t/tl: # If the signatures are equal, then the critical pairs # are not normalized. For some reason, Faug\303\250re and # Stegers do not include this in their algorithm for # computing critical pairs. The Rewrite Rules catch # this later, but it is better to checked here. if Sig_Prod(u1,k)=Sig_Prod(u2,l) then return {}: end if: if First_Sig_Min(Sig_Prod(u1,k),Sig_Prod(u2,l),ordering) then kcp := l: lcp := k: # SWAP u1,u2 := u2,u1: else kcp := k: lcp := l: end: Sigk1:=Sig(kcp): Sigk2:=Sig(lcp): k1 := Sigk1[2]: k2 := Sigk2[2]: # apply the criteria to determine whether this critical pair is # necessary (i.e., whether it is already normalized) # current implementation uses Faug\303\250re's formulation, # with Stegers' marked as "TS:" #TS: if Top_Reducible(u1*Sigk1[1],G[k1+1],ordering) then if Top_Reducible(u1*Sigk1[1],G_prev,ordering) then checked_pairs := checked_pairs union {k,l}: return {}: end if: #TS: if k2+1<=nops(G) and Top_Reducible(u2*Sigk2[1],G[k2+1],ordering) then if k2=i and Top_Reducible(u2*Sigk2[1],G_prev,ordering) then checked_pairs := checked_pairs union {k,l}: return {}: end if: return {[t,u1,kcp,u2,lcp]}: end: # Crit_Pair SPols considers each critical pair generated by Crit_Pair, checks whether it is rewritable according to a different rule, which means that a previous computation can be used to show that the corresponding S-polynomial reduces to zero, and thus can be skipped. If not, it computes the S-polynomial. Since this S-polynomial is necessary, it adds it to the Gr\303\266bner basis (by adding it to lps). # B is a set of critical pairs SPols:=proc(B) # DESCRIPTION OF LOCAL VARIABLES # F holds the indices of new polynomials added to the basis r # cps is the list of critical pairs under consideration # cp is the current critical pair under consideration # [t,u,k,v,l]=cp # c1,c2 are the leading coefficients of lps[k],lps[l], resp. # s is the S-polynomial # N is the index of the last polynomial added to the basis, # assuming that one has been added # heap_sorter is used to sort the critical pairs local F,cps,cp,t,u,k,v,l,c1,c2,s,N,heap_sorter: # initialize record keeping F:={}: # the ordering here is reversed because heap[extract] # selects the MAXIMUM elements of the heap, and we want # it to select the MINIMUM instead heap_sorter:=proc(a,b) return Groebner[TestOrder](b[1],a[1],ordering): end: # heap_sorter cps:=heap[new](heap_sorter,op(B)); # loop until no critical pairs remain while not heap[empty](cps) do # pick the smallest critical pair in the set cp := heap[extract](cps): # obtain information necessary to compute S-polynomial u := cp[2]: k := cp[3]: v := cp[4]: l := cp[5]: # is this S-polynomial rewritable? if so, skip it if (not Is_Rewritable(u,k)) and (not Is_Rewritable(v,l)) then c1 := Groebner[LeadingCoefficient](Poly(k),ordering): c2 := Groebner[LeadingCoefficient](Poly(l),ordering): s := simplify(c2 * u * Poly(k) - c1 * v * Poly(l)): if characteristic > 0 then s := s mod characteristic: end if: # check whether we actually need to add this S-polynomial # to the set against the rewriting records kept -- # if so, add it; otherwise, ignore it if s<>0 then Add_Labeled_Poly([Sig_Prod(u,k),s,{k,l}]): Add_Rule(): F := F union {last_lp}: end if: # s<>0 end if: # rewritable end do: # cp # we sort F in Reduction, so sorting is unnecessary here return F: end: # SPols Reduction takes the list of polynomials provided (indexed by F), and reduces them modulo the previous basis. The result is a list of indices. # F is a list indexing the polynomials added to the basis in this step # G_prev is the previous Gr\303\266bner basis # G_curr is the set of generators for the current ideal Reduction:=proc(F,G_prev,G_curr) # DESCRIPTION OF LOCAL VARIABLES # to_do is the list of polynomials yet to be reduced # completed is the set of polynomials that have been reduced # k indexes the current iteration through to_do # h is the reduced form of lps[k] # newly_completed contains the indices of polynomials completed on this iteration # redo contains the indices of polynomials that still need reduction # heap_sorter is used to sort to_do # reducers is the polynomials of the previous Groebner basis local to_do,completed,k,h,newly_completed,redo,j,heap_sorter,reducers: # Initialization reducers:=[seq(Poly(seqi),seqi in G_prev)]: completed := {}: # sort to_do by increasing signature heap_sorter := proc(a,b) return not Sig_Comp(a,b) end: to_do := heap[new](heap_sorter,op(F)): # loop until no polynomials left to reduce while not heap[empty](to_do) do # compute the normal form of the first polynomial in the list k := heap[extract](to_do): h := Groebner[NormalForm](expand(Poly(k)),reducers,ordering): lps[k] := [Sig(k),h,lps[k][3]]: # See if this polynomial can be topreduced newly_completed, redo := Top_Reduction(k, G_prev, G_curr union completed): # update records completed := completed union newly_completed: for j in redo do heap[insert](j,to_do): end do: # j end do: # dummy return completed: end: # Reduction Top_Reduction checks whether a polynomial LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSJrRidGL0YyL0YzUSdub3JtYWxGJy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYn in the current basis can be top-reduced modulo a polynomial in the previous basis. It returns a list of two sets. The first is {k} if no top-reduction took place, where k is the index of the polynomial; otherwise, it is {}. The second is a set of polynomials that need further reduction, in case LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSJrRidGL0YyL0YzUSdub3JtYWxGJy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYn was not top reduced. Whether the set contains only one element (k) or two (k and the polynomial that reduced it) depends on which has the smaller signature. # k is the index of a polynomial in the basis # G_prev is the previous Gr\303\266bner basis # G_curr is the generators of the current ideal Top_Reduction:=proc(k,G_prev,G_curr) # DESCRIPTION OF LOCAL VARIABLES # p is the unlabeled polynomial of lps[k] # J is a set of reductors for lps[k] (but this set only has one el, # unless it is empty) # j is the index of a reductor lps[j] for lps[k], if one exists # q is the unlabeled polynomial of the reductor lps[j] # u is a multiplier for lm(q) to obtain lm(p) local p,J,j,q,u: # the following should not happen with inputs # that are regular sequences if Poly(k)=0 then cps_of_zero_reductions := cps_of_zero_reductions union {lps[k][3]}: if warn_on_zero_reduction then printf("Warning, reduction to zero for pair (%d,%d): F5 may not terminate\134n", lps[k][3]): end if: return {},{}: end if: # initialization; find reductor of lps[k], if one exists p := Poly(k): J := Find_Reductor(k,G_prev,G_curr): if J = {} then # no reductor; make the leading coefficient 1 # nonzero characteristic? if characteristic > 0 then p := p * ((1 / Groebner[LeadingCoefficient](p,ordering)) mod characteristic) mod characteristic: else p := p / Groebner[LeadingCoefficient](p,ordering): lps[k] := [Sig(k),p,lps[k][3]]: end if: return {k},{}: else # a reductor lps[j]! we will top-reduce lps[k] over lps[j] # and notify of the need for additional checks as appropriate j := J[1]: q := Poly(j): u := Groebner[LeadingMonomial](p,ordering) / Groebner[LeadingMonomial](q,ordering): # nonzero characteristic? if characteristic > 0 then p := Groebner[LeadingCoefficient](q,ordering) * p - Groebner[LeadingCoefficient](p,ordering) * u * q mod characteristic: else p := Groebner[LeadingCoefficient](q,ordering)* p - Groebner[LeadingCoefficient](p,ordering) * u * q: end: if First_Sig_Min(Sig_Prod(u,j),Sig(k),ordering) then lps[k] := [Sig(k),p,lps[k][3]]: return {},{k}: else Add_Labeled_Poly([Sig_Prod(u,j),p,lps[k][3]]): Add_Rule(): return {},{k,last_lp}: end if: end if: end: # Top_Reduction Find_Reductor looks for a normalized polynomial in GPrev that will top-reduce . The result is empty if all the polynomials that can reduce satisfy the normalization or if the computation can be reused (rewriting criteria). # k is the index of the polynomial to reduce # G_prev is the previous Gr\303\266bner basis # G_curr is the generators of the current ideal Find_Reductor:=proc(k,G_prev,G_curr) # DESCRIPTION OF LOCAL VARIABLES # t is the leading monomial of lps[k] # j is the current iteration through G_curr # tprime is the leading monomial of lps[j] # mu_j is the multiplier for lps[j], nu_j is the syzygy basis element # u is the multiplier for tprime to obtain t local t,j,tprime,mu_j,nu_j,u: t := Groebner[LeadingMonomial](Poly(k),ordering): # loop through the polynomials of the running GB # if the lm of lps[j] divides the lm of lps[k], AND their syzygy # has not been disposed of, AND (u,j) is not rewritable, AND # u*lps[j] is not top-reducible, then return j for j in G_curr do tprime := Groebner[LeadingMonomial](Poly(j),ordering): if divide(t,tprime) then u := t/tprime: mu_j := Sig(j)[1]: nu_j := Sig(j)[2]: if Sig_Prod(u,j)<>Sig(k) and not Is_Rewritable(u,j) and not Top_Reducible(u*mu_j,G_prev,ordering) #TS: Sig_Prod(u,j)=Sig(k) or Is_Rewritable(u,j) #TS: or (nu_j<nops(G) and Top_Reducible(u*mu_j,G[nu_j+1],ordering)) then #TS: # do nothing; return on "else" return {j}: end if: end if: end do: # j # no reductor return {}: end: # Find_Reductor Add a rule for a newly generated polynomial. Add_Rule:=proc() # DESCRIPTION OF LOCAL VARIABLES # j is the index of a newly generated labeled polynomial # t is the multiplier for the signature of lps[j] # i is the syzygy index for the signature of lps[j] local j,t,i: j := last_lp: t := Sig(j)[1]: i := Sig(j)[2]: # okay, add the rule! last_Rules[i] := last_Rules[i] + 1: Rules[i][last_Rules[i]] := [t,j]: end: # Add_Rule Is there a rule that rewrites u\342\213\205lps[k] (aside from the rule(s) based on itself)? # u is a term # k is the index of a labeled polynomial Is_Rewritable:=proc(u,k) local kprime: kprime := Find_Rewriting(u,k): return evalb(k<>kprime): end: # Is_Rewritable Find the first rule that rewrites u\342\213\205lps[k]. The search proceeds through the rules for the incremental Gr\303\266bner basis that generated . # u is a term # k is the index of a labeled polynomial Find_Rewriting:=proc(u,k) # v is the monomial multiple of the syzygy for the kth labeled poly # l is the index in the signature for the kth labeled poly # j is a loop variable # t is a multiplier for the Rule of index j # jprime is the index in the Rule of index j local mu,nu,j,t,jprime: mu := Sig(k)[1]: nu := Sig(k)[2]: for j from 1 to last_Rules[nu] do t := Rule_Element(nu,j)[1]: jprime := Rule_Element(nu,j)[2]: if divide(u*mu,t) then return jprime: end if: end do: # j return k: end: # Find_Rewriting Add_Labeled_Poly was a kludge to workaround Maple's lack of linked lists. Now I'm using a table, so this does very little. An optimization might remove it altogether. # p is a labeled polynomial Add_Labeled_Poly:=proc(lp) last_lp := last_lp + 1: lps[last_lp]:=lp: end: # Add_Labeled_Poly This is a support function specified by Faug\303\250re. It gives the signature of a labeled polynomial. Sig:=proc(k::posint) return lps[k][1]: end: # Sig This is a support function specified by Faug\303\250re. It gives the unlabeled polynomial of a labeled polynomial. Poly:=proc(k::posint) return lps[k][2]: end: # Poly This is a support function implementing multiplication of a labeled polynomial's signature by a monomial. Sig_Prod:=proc(t::monomial,k::posint) return [t*lps[k][1][1],lps[k][1][2]]; end: # Sig_Prod This is a support function that compares two signatures and returns the smaller (used for sorting). Sig_Comp:=proc(i,j) return First_Sig_Min(Sig(i),Sig(j)); end: # Sig_Comp; This is a support function that returns true if a signature s1 is smaller than (or equal to) another signature s2, with respect to the term ordering tord. First_Sig_Min:=proc(s1,s2) local i,j: i := s1[2]: j := s2[2]: if i > j then return true: end if: if i = j and Groebner[TestOrder]( Groebner[LeadingMonomial](s1[1],ordering), Groebner[LeadingMonomial](s2[1],ordering),ordering) then return true: end if: return false: end: # First_Sig_Min This is a support function that returns true if the monomial t is divisible by the leading monomial of a polynomial indexed by G, with respect to tord. Top_Reducible:=proc(t,G) local i: for i in G do if divide(t,Groebner[LeadingMonomial](Poly(i),ordering)) then return true: end: # if end do: # i return false: end: # Top_Reducible For some reason I don't remember, I store the rules in a manner that is not quite the same as that specified by Faug\303\250re, and this function accesses the elements correctly. Rule_Element:=proc(rule_number, element_number) return Rules[rule_number][last_Rules[rule_number]-element_number+1]: end: # Rule_Element Diagnostic functions... Report_Number_Of_SPolys:=proc() return number_of_spolys: end: # Report_Number_Of_SPolys Report_Number_Of_Zero_Reductions:=proc() return nops(cps_of_zero_reductions): end: # Report_Number_Of_Zero_Reductions Report_Set_Of_Critical_Pairs:=proc() return critical_pairs_used: end: # Report_Set_Of_Critical_Pairs Report_Pairs_Giving_Zero_Reductions:=proc() return cps_of_zero_reductions: end: # Report_Pairs_Giving_Zero_Reductions Report_Rules:=proc() return Rules: end: # Report_Rules Report_Rule_Counts:=proc() return last_Rules: end: # Report_Rule_Counts Report_Labeled_Polys:=proc() return lps: end: # end: # f5_computation

f5_utility := module() export Generate_Cyclicn, is_gb, BC2: local sort_by_ascending_lcm: Generate_Cyclicn:=proc(n) local result, intermediate, facs, addem, base_set, i, j, k: result := []: base_set := [seq(x[seqctr],seqctr=1..n)]: for i from 1 to n-1 do addem := 0: for j from 1 to n do facs := 1: for k from j to i+j-1 do if k<=n then facs := facs * x[k]: else facs := facs * x[k-n]: end if: end do: addem := addem + facs: end do: result := [op(result), addem]: end do: result := [op(result), product(x[ctr],ctr=1..n)-1]: return result: end: # Generate_Cyclicn is_gb:=proc(F::list(polynom),tord) local i,j,Sij,Rij,bad_cps,good_cps,check_cps,lts,cp: lts := [seq(Groebner[LeadingMonomial](F[ctr],tord), ctr=1..nops(F))]: check_cps := [op(combinat[choose]({seq(ctr,ctr=1..nops(F))},2))]: check_cps := sort_by_ascending_lcm(check_cps,lts,tord): bad_cps := {}: good_cps := {}: for cp from 1 to nops(check_cps) do i := check_cps[cp][1]: j := check_cps[cp][2]: if gcd(lts[i],lts[j])<>1 and BC2(i,j,lts,good_cps)={} then #printf("Checking (%d,%d) %a\134n",i,j,lcm(lts[i],lts[j])); Sij:=Groebner[SPolynomial](F[i],F[j],tord): Rij:=Groebner[NormalForm](Sij,F,tord): if Rij<>0 then bad_cps := bad_cps union {{i,j}}; else good_cps := good_cps union {{i,j}}; end if: else good_cps := good_cps union {{i,j}}; end if: end do: # cp return evalb(bad_cps={}),bad_cps end: # is_gb sort_by_ascending_lcm:=proc(cps,lts,tord) local sorter: sorter := proc(a,b) return Groebner[TestOrder](lcm(lts[a[1]],lts[a[2]]), lcm(lts[b[1]],lts[b[2]]),tord): end: # sorter return sort(cps,sorter): end: # sort_by_ascending_lcm BC2 := proc(i::posint,j::posint,lts::list(monomial), cps::set(set(posint))) local other_cps,lt,lcm_ij,k: lcm_ij := lcm(lts[i],lts[j]): for k from 1 to nops(lts) do lt := lts[k]: if k<>i and k<>j and divide(lcm_ij,lt) and {i,k} in cps and {j,k} in cps then return {k}: end if: end do: # k return {}: end: # BC2 end: # f5_utility

tord:=tdeg(x,y,z,t); FEx:=f5_computation[Basis]([ y*z^3-x^2*t^2, x*z^2-y^2*t, x^2*y-z^2*t ],tord): nops(FEx); f5_utility[is_gb](FEx,tord); f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions]();

It takes a 2.8GHz Intel Core 2 Duo machine approximately 15 seconds to compute the basis. A7583 := Ore_algebra[poly_algebra](x[1],x[2],x[3],x[4],x[5], characteristic=7583): tord7583 := Groebner[MonomialOrder](A7583, tdeg(x[1],x[2],x[3],x[4],x[5])): start_time:=time(): Cyclic5:=f5_computation[Basis](f5_utility[Generate_Cyclicn](5),tord7583): total_time:=time()-start_time: trunc(total_time/60.), round(total_time) mod 60; nops(Cyclic5); Cyclic5:=Groebner[InterReduce](Cyclic5,tord7583): nops(Cyclic5); f5_utility[is_gb](Cyclic5,tord7583); f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions]();

The following computation took 14 minutes, 38 seconds on an Intel Core 2 Duo running Fedora Core 7. This is an abysmal 60-fold increase over Cyclic 5; compare to Stegers' 40-fold increase (from .179 seconds to 7.07). start_time:=time(): Cyclic6:=f5_computation[Basis](f5_utility[Generate_Cyclicn](6),tdeg(x[1],x[2],x[3],x[4],x[5],x[6])): total_time:=time()-start_time: trunc(total_time/60.), round(total_time) mod 60; nops(Cyclic6); Cyclic6:=Groebner[InterReduce](Cyclic6,tdeg(x[1],x[2],x[3],x[4],x[5],x[6])): nops(Cyclic6); f5_utility[is_gb](Cyclic6,tdeg(x[1],x[2],x[3],x[4],x[5],x[6])); Notice that here we get 16 zero reductions, just as Stegers reports on page 46 of his thesis. We seem to have generated a few more polynomials than Stegers: he reports 218 reductions, whereas we have 238 total. This seems odd, because we have generated only 198 S-polynomials. f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions](); # WARNING: Uncommenting the following line will give you # a lot of output! #f5_computation[Report_Set_Of_Critical_Pairs]();

Eco6 := [ x1*x2*x6 + x2*x3*x6 + x3*x4*x6 + x4*x5*x6 + x1*x6*h - h^3, x1*x3*x6 + x2*x4*x6 + x3*x5*x6 + x2*x6*h - 2*h^3, x1*x4*x6 + x2*x5*x6 + x3*x6*h - 3*h^3, x1*x5*x6 + x4*x6*h - 4*h^3, x5*x6 - 5*h^2, x1 + x2 + x3 + x4 + x5 + h ]; A2 := Ore_algebra[poly_algebra](x1,x2,x3,x4,x5,x6,h, characteristic=2); tord2 := Groebner[MonomialOrder](A2,tdeg(x1,x2,x3,x4,x5,x6,h)); A7583 := Ore_algebra[poly_algebra](x1,x2,x3,x4,x5,x6,h, characteristic=7583); tord7583 := Groebner[MonomialOrder](A7583, tdeg(x1,x2,x3,x4,x5,x6,h)); On a 2.8GHz Intel Core 2 Duo machine, this reports about 9 seconds. start_time:=time(): Eco6_Basis2 := f5_computation[Basis](Eco6,tord2): total_time:=time()-start_time: trunc(total_time/60.), round(total_time) mod 60; For both characteristic 2 and characteristic 7583 the implementation computes 16 zero reductions, the same as Stegers' implementation computes, but more than Faug\303\250re's 7. Faug\303\250re writes, [the number of useless critical pairs] does not depend on the implementation (at least for the F5 algorithm). Argh. f5_computation[Report_Number_Of_Zero_Reductions](), f5_computation[Report_Number_Of_SPolys](); f5_computation[Report_Pairs_Giving_Zero_Reductions](); nops(Eco6_Basis2); On a 2.8GHz Intel Core 2 Duo, the following timing has given about 16 seconds. start_time:=time(): Eco6_Basis7583 := f5_computation[Basis](Eco6,tord7583): total_time:=time()-start_time: trunc(total_time/60.), round(total_time) mod 60; f5_computation[Report_Number_Of_Zero_Reductions](), f5_computation[Report_Number_Of_SPolys](); Studying the pairs that give zero reductions... f5_computation[Report_Pairs_Giving_Zero_Reductions](); #lps:=f5_computation[Report_Labeled_Polys](): #seq(Groebner[LeadingMonomial](lps[ctr][2],tord7583),ctr in [8,23,47]); #seq(Groebner[LeadingMonomial](simplify(lps[ctr][2]/h^2),tord7583),ctr in [8,23,47]); rule_counts:=f5_computation[Report_Rule_Counts](); rewrite_rules:=f5_computation[Report_Rules](): If you want to verify that we really did compute a Gr\303\266bner basis, uncomment these lines & run them. Of course, you can do this with Eco6_Basis2 as well. Eco6_Interreduced:=Groebner[InterReduce](Eco6_Basis7583,tord7583): nops(Eco6_Interreduced); f5_utility[is_gb](Eco6_Interreduced,tord7583); Eco6_Interreduced2:=Groebner[InterReduce](Eco6_Basis2,tord2): nops(Eco6_Interreduced2); f5_utility[is_gb](Eco6_Interreduced2,tord2);

Weispfenning94 := [ y^4+x*y^2*z+x^2*h^2-2*x*y*h^2+y^2*h^2+z^2*h^2, x*y^4+y*z^4-2*x^2*y*h^2-3*h^5, -x^3*y^2+x*y*z^3+y^4*h+x*y^2*z*h-2*x*y*h^3 ]: A := Ore_algebra[poly_algebra](x,y,z,h,characteristic=7583); tord:=Groebner[MonomialOrder](A,tdeg(x,y,z,h)); On a 2.8GHz Intel Core 2 Duo, this timing reports about 38 seconds. start_time := time(): Weispfenning94_Basis:=f5_computation[Basis](Weispfenning94,tord): total_time := time() - start_time: trunc(total_time/60.), round(total_time) mod 60; nops(Weispfenning94_Basis); New code: 4 reductions to zero, 58 S-polynomials. I thought I was sorting the inputs, but I wasn't. Now we match Stegers. Old code: 6 reduction to zero, 59 S-polynomials. I replaced the function Min_CP(), used to find the smallest element in a set of critical pairs, with a Maple heap. This improved performance somewhat, and gave consistency. (Maple seems to handle sets of certain types inconsistently). Older code: Repeated runnings: 11,66; 12,66; 15,66; The difference in each run may be due to Maple's handling of sets, since the result is still a Gr\303\266bner basis. f5_computation[Report_Number_Of_Zero_Reductions](), f5_computation[Report_Number_Of_SPolys](); f5_computation[Report_Pairs_Giving_Zero_Reductions](); Weispfenning94_Interreduced := Groebner[InterReduce](Weispfenning94_Basis,tord): nops(Weispfenning94_Interreduced); f5_utility[is_gb](Weispfenning94_Interreduced,tord);

Gerdt93:=[ h*l-l^2-4*l*s+h*y, h^2*s-6*l*s^2+h^2*z, x*h^2-l^2*s-h^3 ]; A:=Ore_algebra[poly_algebra](h,l,s,x,y,z,characteristic=7583); tord:=Groebner[MonomialOrder](A,tdeg(h,l,s,x,y,z)); Gerdt93_Basis:=f5_computation[Basis](Gerdt93,tord): nops(Gerdt93_Basis); Gerdt93_Interreduced:=Groebner[InterReduce](Gerdt93_Basis,tord): nops(Gerdt93_Interreduced); f5_utility[is_gb](Gerdt93_Interreduced,tord); f5_computation[Report_Number_Of_Zero_Reductions]();

Liu := [ y*z - y*t - x*h + a*h, z*t - z*x - y*h + a*h, t*x - y*t - z*h + a*h, x*y - z*x - t*h + a*h ]; A := Ore_algebra[poly_algebra](x,y,z,t,a,h,characteristic=2); tord := Groebner[MonomialOrder](A,tdeg(x,y,z,t,a,h)); Liu_Basis := f5_computation[Basis](Liu,tord): nops(Liu_Basis); Liu_Interreduced := Groebner[InterReduce](Liu_Basis,tord): nops(Liu_Interreduced); f5_computation[Report_Number_Of_Zero_Reductions](); f5_utility[is_gb](Liu_Interreduced,tord);

Sym3_3:=[ y*z^3+h^3*x-2*h^4,x^3*z+h^3*y-2*h^4,x*y^3+h^3*z-2*h^4 ]; A:=Ore_algebra[poly_algebra](h,x,y,z,characteristic=7583): tord:=Groebner[MonomialOrder](A,tdeg(h,x,y,z)); Sym3_3_Basis:=f5_computation[Basis](Sym3_3,tord): nops(Sym3_3_Basis); Sym3_3_Interreduced:=Groebner[InterReduce](Sym3_3_Basis,tord): nops(Sym3_3_Interreduced); f5_utility[is_gb](Sym3_3_Interreduced,tord); f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions]();

tord:=tdeg(x0,x1,x2,x3,x4,x5,x6): NonRegular:=[ x0^3*x4^6+x4^6+x0^3*x4^3+x4^3+x0^3*x4^2+x4^2, x1^2*x2^5+x2^3+x2^2, x0^3*x3+x0^3+x3+1, x0^3*x1^2+x0^3+x1^2+1 ]: NonRegular:=f5_computation[Basis](NonRegular,tord): f5_utility[is_gb](NonRegular,tord); nops(NonRegular); f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions](); f5_computation[Report_Pairs_Giving_Zero_Reductions](); NonRegular := Groebner[InterReduce](NonRegular,tord): nops(NonRegular);

NonRegular_Generators := [ x0^3*x4^6+x4^6*x5^3+x0^3*x4^3*x5^3+x4^3*x5^6+x0^3*x4^2*x5^4+x4^2*x5^7, x1^2*x2^5+x2^3*x5^4+x2^2*x5^5, x0^3*x1^2+x0^3*x5^2+x1^2*x5^3+x5^5, x0^3*x3+x0^3*x5+x3*x5^3+x5^4 ]: tord:=tdeg(x0,x1,x2,x3,x4,x5): NonRegular:=f5_computation[Basis](NonRegular_Generators,tord): f5_utility[is_gb](NonRegular,tord); nops(NonRegular); f5_computation[Report_Number_Of_SPolys](), f5_computation[Report_Number_Of_Zero_Reductions](); f5_computation[Report_Pairs_Giving_Zero_Reductions]();

#regular: F:=[x^2*y-z*h^2,x*z^2-y*h^2,z^3-x*h^2]: #nonregular: #F:=[x^3*y-x^2*h^2,x*y^3,x^4-x*h^3]: #nonregular, nonhomogeneous: F:=[x^3*y-x^2,x*y^3,x^4-x]: tord:=tdeg(x,y,z,h): f5_computation[Basis](F,tord); f5_computation[Report_Number_Of_Zero_Reductions](), f5_computation[Report_Number_Of_SPolys]() ; f5_computation[Report_Pairs_Giving_Zero_Reductions](); lps:=f5_computation[Report_Labeled_Polys](): lps[7];lps[9]; rls:=f5_computation[Report_Rules](): rls[1]; The S-polynomial for lps[12] and lps[14] is interesting. The algorithm as published by Faug\303\250re and Stegers computes a critical pair, even though the pair is not normalized, because the builders of the S-polynomial have equal signatures! It does not reduce to zero, because Is_Rewritten returns true to SPols. This pair turns out to be a syzygy that is not a principal syzygy, as we see below. lps:=f5_computation[Report_Labeled_Polys](): lps[12]; lps[14]; e:=[seq(lps[ctr][2],ctr=1..14)]: simplify(e[12]-x*e[14]); simplify((-e[6]+y*e[11])+x*(y*e[1]-x*e[10]+e[7])); simplify(e[3]-x*e[5]),e[7]; simplify(-e[6]+y*e[11]+x*(y*e[1]-x*e[10]+e[7])); simplify(-(e[1]-y*e[5])+y*(e[10]-y*e[7])+x*(y*e[1]-x*e[10]+e[7])); simplify((x^3-y^3-1)*e[1] + (x^2*y-x)*e[2] + (-x^2*y+x)*e[3]); e[5],simplify(-x*e[1]+y*e[3]);

Insert the name of your save library between quotes. These have been commented out to avoid accidental execution (e.g., Edit->Execute->Worksheet) ## example from my system: ##savelibname := "/Users/JackPerry/Library/maplelibs/mylibs.lib"; #savelibname := "/Users/JackPerry/Library/maplelibs/mylibs.lib"; ## uncomment next line only if you do not have a library already #march('create',savelibname); #savelib('f5_computation'); #march('list',savelibname); ## uncomment next line if you wish to delete it from your system #march('delete',savelibname,"f5_computation.m"); Execute the following lines only if you wish to test that the module saved successfully. #restart; ## to avoid having to do this automatically, place it into your ## Maple initialization file (maple.ini, .mapleinit; see Maple's ## documentation for details) #savelibname := "/Users/JackPerry/Library/maplelibs/mylibs.lib"; #libname := libname, savelibname; #with(f5_computation);