## PolynomialIdeals v2 - August 2005 ## Primarily authored and maintained by ## Roman Pearce, Simon Fraser University ## contact rpearcea@cecm.sfu.ca ## Joint work with ## Michael Monagan ## Jeff Farr (Groebner Walk, VanishingIdeal) ## Jennifer De Kleine (documentation, testing) ## Juergen Gerhard # todo: optimize inter_reduce_gb # optimize rat_nf and mod_nf # make FGLM call fglm_system and XLA:-LUSolve if 0 < char < 2^25 # note to maintainer(s) : # # SPEED CRITICAL means that a procedure is used extensively as # a subroutine in some algorithm, and its execution time # drastically affects overall performance # please do not add argument checking to such procedures, # and do not export them for the user # Groebner basis routines : # call PolynomialIdeals:-GroebnerBasis with method='buchberger' or method='f4' # rather than use the procedures below # # PolynomialIdeals:-Buchberger : Buchberger algorithm (no RootOfs) # PolynomialIdeals:-f4_algorithm:-f4 : F4 algorithm (no RootOfs or parameters) # PolynomialIdeals:-FGLM : FGLM # PolynomialIdeals:-GroebnerWalk : Groebner Walk # # The Buchberger algorithm can express the Groebner basis in terms of the input. # My implementation of FGLM uses a custom linear solver. # Groebner Walk is my first implementation, and it can probably be improved. # Interesting and useful stuff # # PolynomialIdeals:-substitute_rootofs # - replace RootOfs with variables and return their minimal polynomials # along with forward and reverse substitutions. Handles nesting :) # # PolynomialIdeals:-substitute_everything # - substitutes for constants (Pi=_pi, etc) and correctly rewrites fractional extensions # ie: Pi^(1/3) + Pi^(1/2) -> _s^2 + _s^3 where _s = Pi^(1/6) # # PolynomialIdeals:-linear_solve # - incremental Gaussian elimination, treats polynomials as linear in their monomials # # PolynomialIdeals:-similar_orders # - select term orders most similar to a given order # # PolynomialIdeals:-gbwalk_algorithm:-tord2mat # - construct a matrix representation of a term order, usually with one line :) # # PolynomialIdeals:-gbwalk_algorithm:-initial_form # - compute the initial form of a polynomial w.r.t a weight vector # # PolynomialIdeals:-ideal_algorithms:-monomial_basis # - compute a basis of the quotient ring as a vector space # # PolynomialIdeals:-ideal_algorithms:-vectorize # - convert a polynomial to a vector of coefficients # # PolynomialIdeals:-ideal_algorithms:-qr_dim # - just compute the dimension of the quotient space # input types for various routines #$define POLY_T Or(radalgfun,list(radalgfun),set(radalgfun)) # type check for field extensions #$define ALG_EXT Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))) # SPEED CRITICAL - divide out integer content # in order to be inlined, this procedure needs to be global (not savelib'd) iprimpart := proc(f) option inline; `if`(f=0, 0, f/icontent(f)) end proc: # check for exact monomial division - used by F4 mdivide := proc(m1, m2) option inline; evalb(denom(m1/m2)=1) end proc: unprotect(PolynomialIdeals); PolynomialIdeals := module() option package, `Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2005`; description "a package for manipulating commutative polynomial ideals"; local gb_table, triset_table, termorder_procs, termorder_types, evaluate_termorder, tord_split, is_elim_order, elim_order, similar_orders, nextmon_procs, evaluate_nextmon, substitute_rootofs, collect_rootof_powers, substitute_everything, normalform_procs, spolynomial_procs, buchberger_algorithm, Buchberger, InterReduce, InterReduceGB, xratrecon, f4_algorithm, fglm_system, triset_algorithm, linear_solve, fglm_algorithm, gbwalk_algorithm, LeadingCoefficient, ideal_algorithms, vanishing_ideal, PointSort; global `type/PolynomialIdeal`, `print/POLYNOMIALIDEAL`, `convert/PolynomialIdeal`, `expand/POLYNOMIALIDEAL`, `mod/POLYNOMIALIDEAL`, `latex/POLYNOMIALIDEAL`; export _pexports, PolynomialIdeal, `<,>`, IdealInfo, Generators, XLA, TriangularSet, FGLM, GroebnerWalk, SuggestVariableOrder, GroebnerBasis, ModuleGB, LeadingMonomial, TrailingMonomial, InitialForm, WeightedDegree, TermOrderMatrix, NormalForm, UnivariatePolynomial, EliminationIdeal, Intersect, Saturate, Contract, Simplify, Homogenize, MaximalIndependentSet, HilbertDimension, IsZeroDimensional, NumberOfSolutions, MonomialBasis, IsProper, IdealMembership, IdealContainment, `in`, `subset`, RadicalMembership, Add, Multiply, Quotient, ZeroDimensionalDecomposition, EquidimensionalDecomposition, Radical, IsRadical, IsPrime, IsPrimary, IsMaximal, PrimeDecomposition, PrimaryDecomposition, VanishingIdeal, Operators; _pexports := () -> [op({exports(PolynomialIdeals)} minus {':-_pexports', ':-XLA', ':-TriangularSet', ':-FGLM', ':-GroebnerWalk', ':-SuggestVariableOrder', ':-GroebnerBasis', ':-ModuleGB', ':-LeadingMonomial', ':-TrailingMonomial', ':-InitialForm', ':-WeightedDegree', ':-TermOrderMatrix', ':-NormalForm', ':-MonomialBasis' })]; # Groebner basis table kludge # input : generating set, characteristic # output : unique table for storing Groebner bases # DON'T MODIFY OPTION REMEMBER gb_table := proc() option remember; table() end proc; triset_table := proc() option remember; table() end proc; # Ideal Constructor # input : generators (either directly, or in a list or set) # optional arguments (in the form of equations) # # optional arguments # variables=set(name) : specify the ring variables # parameters=set(name) : specify indeterminates of the coefficient field # this is converted to a 'variables'= # characteristic=nonnegint : specify the ring characteristic # known_groebner_bases=list : specify a list of known Groebner bases # list elements should be of the form # termorder=gb, ie: plex(x,y)=[y^2, x-1] # # default variables are all names not appearing in a RootOf or radical # default characteristic is zero # # output : POLYNOMIALIDEAL(sequence of generators, 'variables'=set(name), # 'characteristic'=nonnegint, 'known_groebner_bases'=gb_table::table) # # example : # PolynomialIdeal(x^2-1, y-z, characteristic=3, variables={x,y}, known_groebner_bases=[plex(x,y)=[y-z,x^2-1]) # # POLYNOMIALIDEAL(x^2-1, y-z, 'variables'={x,y}, 'characteristic'=3, # 'known_groebner_bases'=table([ plex(x,y) = [y-z, x^2-1] ])); # # additional note : # - procedures must be able to work with malformed ideals # - ordering of generators is Maple's set ordering # - ordering of options is Maple's table ordering PolynomialIdeal := proc() option system, remember; local gens, opts, defopts, KGB, newGB, vars, char, i, unknown; # kludge, continued # - if the only argument is a single ideal, reconcile its Groebner basis table with the system # - used for ideals which have been saved and then loaded into a different Maple session. if nargs = 1 and type(args[1], 'PolynomialIdeal') then if hasoption(select(type, [op(args[1])], 'equation'), 'known_groebner_bases'=table, 'KGB') then return procname(args[1], 'known_groebner_bases'=op(op(KGB))) else return procname(op(args[1])) end if; end if; opts, gens := selectremove(type, [args], 'equation'); gens := map(proc(a) `if`(type(a, 'Or'('list','set', 'PolynomialIdeal')), op(a), a) end proc, gens); defopts, gens := selectremove(type, gens, 'equation'); unknown := remove(type, gens, 'radalgfun'); if nops(unknown) > 0 then error "invalid arguments(s)", op(unknown) end if; vars := indets(gens,'name') minus indets(indets(gens,'Or'('RootOf','radical','And'('name','constant'))), 'name'); defopts := ['characteristic'=0, 'variables'=vars, op(remove(i->evalb(lhs(i)='known_groebner_bases'),defopts))]; newGB := `if`(hasoption(opts, 'known_groebner_bases'='Or'('list'('equation'),'set'('equation')), 'i', 'opts'), i, []); opts := table([op(defopts), op(opts), `if`(hasoption(opts, 'parameters'='set'('name'), 'i'), 'variables'=vars minus i, NULL)]); unassign('opts[parameters]'); char := opts['characteristic']; if not (char=0 or isprime(char)) then error "the characteristic must be zero or a prime" end if; vars := opts['variables']; gens := `if`(char=0, {seq(iprimpart(i), i=gens)}, {op(gens)} mod char) minus {0}; if not type(gens, 'set'('polynom'('radalgfun', vars))) then error "the generators are not polynomial in the variables" end if; # get whatever Groebner basis table the system knows about # and add any Groebner bases that the user explicity specified KGB := gb_table(gens, char); for i in newGB do KGB[lhs(i)] := rhs(i) end do; 'POLYNOMIALIDEAL'(op(gens), op(op(op(opts))), 'known_groebner_bases'=op(KGB)); end proc; # constructor shortcut `<,>` := eval(PolynomialIdeal,1): # this is an API to extract things from the ideal data structure # all programs should call these routines IdealInfo := module() option package; export Generators, Characteristic, Variables, NonVariables, Parameters, DefaultMonomialOrder, KnownGroebnerBases; Generators := proc(J::PolynomialIdeal, $) option system, remember; remove(type, {op(J)}, 'equation') minus {0} end proc; # returns characteristic. default is zero if missing or invalid Characteristic := proc(J::PolynomialIdeal, $) option system, remember; local c; `if`(hasoption(select(type,[op(J)], 'equation'), 'characteristic'='nonnegint', 'c'), c, 0) end proc; # returns the variables of the ideal, trusts 'variables'= # otherwise the default is every name not appearing in a RootOf or radical Variables := proc(J::PolynomialIdeal, $) option system, remember; local v; if hasoption(select(type,[op(J)], 'equation'), 'variables'='set'('name'), 'v') then v else indets(Generators(J), 'name') minus NonVariables(J) end if; end proc; # returns a set of the names which appear in RootOfs or radical expressions. NonVariables := proc(J::PolynomialIdeal, $) option system, remember; indets(indets(Generators(J), 'Or'('RootOf','radical', 'And'('name','constant'))), 'name'); end proc; # returns the parameters, basically any name which is not a variable Parameters := proc(J::PolynomialIdeal, $) option system, remember; indets(Generators(J), 'name') minus Variables(J); end proc; # suggests term order which will be fast to compute # you can optionally specify the type of term order type or the variables # it will select a term order for which the shortest (in bit length) # Groebner basis is known, falling back on SuggestVariableOrder # DON'T USE OPTION REMEMBER DefaultMonomialOrder := proc(J::PolynomialIdeal) local vars, tordtype, KGB, i, unknown; vars,unknown := selectremove(type, [args][2..-1], 'set'('name')); vars := [op(vars), IdealInfo:-Variables(J)][1]; tordtype,unknown := selectremove(member, unknown, termorder_types); tordtype := [op(tordtype), 'anything'][1]; # this procedure just has to ignore garbage arguments # otherwise we have to program all the cases in the # call from GroebnerBasis # if nops(unknown)>0 then # error "invalid argument(s)", op(unknown) # end if: if hasoption(select(type, [op(J)], 'equation'), 'known_groebner_bases'=table, 'KGB') then if tordtype = 'anything' then KGB := sort([seq(`if`(indets(i,'name')=vars, [op(i), length(KGB[op(i)])], NULL), i=indices(KGB))], proc(a,b) option inline; evalb(a[2] < b[2]) end proc) else KGB := sort([seq(`if`(op(0,i[1])=tordtype and indets(i,'name')=vars, [op(i), length(KGB[op(i)])], NULL), i=indices(KGB))], proc(a,b) option inline; evalb(a[2] < b[2]) end proc) end if; if nops(KGB) > 0 then return KGB[1][1] end if; end if; if member(tordtype, {'anything', 'wdeg', 'lexdeg', 'matrix', 'prod'}) then 'tdeg'(SuggestVariableOrder(J, vars)) else tordtype(SuggestVariableOrder(J, vars)) end if end proc; # returns a set of termorders for which Groebner bases are known. # call GroebnerBasis with the termorder to get the actual Groebner basis. # DON'T USE OPTION REMEMBER KnownGroebnerBases := proc(J::PolynomialIdeal, X::set(name), $) local KGB, i; if hasoption(select(type, [op(J)], 'equation'), 'known_groebner_bases'=table, 'KGB') then `if`(nargs >= 2, {seq(`if`(indets(i,'name')=X, op(i), NULL), i=indices(KGB))}, {seq(op(i), i=indices(KGB))}) else {} end if end proc; end module; # a shortcut to get the generators Generators := eval(IdealInfo:-Generators, 1): ## ## Termorder Procedures # SPEED CRITICAL (for the resulting procedures, not the root procedures in the table) # we describe term orders using functions, for example: plex(x,y,z) # to computations we use procedures which evaluate the leading term of a polynomial # and return the sequence (LEADING_COEFFICIENT, LEADING_MONOMIAL) # the evaluate_termorder command constructs these procedures from their descriptions # evaluate_termorder(plex(x,y,z)) calls termorder_procs[plex](x,y,z) # the indices of this table are assumed to be all of the valid term orders # furthermore it is assumed that the set of names appearing in a term order # description is precisely equal to the term order variables termorder_procs := table([ # pure lexicographic order 'plex' = proc() option remember; local VARS; if type([args],'list'('name')) then subs(VARS=[args], proc(f) local lm; lcoeff(f,VARS,'lm'), lm end proc ) else error "sequence of variable names expected" end if end proc, # graded lexicographic order 'grlex' = proc() option remember; local VARS_LIST, SUBS_SET, i, h; if type([args],'list'('name')) then subs({VARS_LIST=[args], SUBS_SET={seq(i=i*h,i={args})}}, proc(f) local lm; lcoeff(lcoeff(subs(SUBS_SET,f), h), VARS_LIST, 'lm'), lm end proc ) else error "sequence of variable names expected" end if end proc, # graded-reverse lexicographic order 'tdeg' = proc() option remember; local SRAV_LIST, SUBS_SET, i, h; if type([args],'list'('name')) then subs({SRAV_LIST=[seq(args[nargs-i+1],i=1..nargs)], SUBS_SET={seq(i=i*h,i={args})}}, proc(f) local lm; tcoeff(lcoeff(subs(SUBS_SET,f), h), SRAV_LIST, 'lm'), lm end proc ) else error "sequence of variable names expected" end if end proc, # graded-reverse elimination order # note : lexdeg([a,b,c],[x,y,z]) is equivalent to doing tdeg(a,b,c) then tdeg(x,y,z) # this could probably be done with a single substitution, but I'm lazy 'lexdeg' = proc() option remember; local LT, i; LT := [seq(evaluate_termorder(`if`(nops(i)=1, 'plex'(op(i)), 'tdeg'(op(i)))), i=[args])]; proc(f) local lc, lm, lmt, lt; lc := f; lm := 1; for lt in LT do lc, lmt := lt(lc); lm := lm*lmt end do; lc, lm end proc; end proc, # general product order. evaluates the leading term using other termorders, from left to right # typical usage is something like : prod(tdeg(x,y,z), plex(u,v), grlex(t,w)) 'prod' = proc() option remember; local LT, i; LT := [seq(evaluate_termorder(i), i=[args])]; proc(f) local lc, lm, lmt, lt; lc := f; lm := 1; for lt in LT do lc, lmt := lt(lc); lm := lm*lmt end do; lc, lm end proc; end proc, # weighted-degree order, breaks ties using inverse lexicographic order 'wdeg' = proc(W::list(rational), vars::list(name)) option remember; local weights, SUBS_SET, SRAV_LIST, i, w; if nops(W) <> nops(vars) then error "there must be exactly one weight per variable" else weights := lcm(seq(denom(i), i=W))*W; subs({SUBS_SET={seq(vars[i]=vars[i]*w^weights[i], i=1..nops(vars))}, SRAV_LIST=[seq(vars[-i], i=1..nops(vars))]}, proc(f) local lm; tcoeff(lcoeff(subs(SUBS_SET, f), w), SRAV_LIST, 'lm'), lm end proc) end if end proc, # matrix order, breaks ties with successive rows # and with lexicographic order at the very end 'matrix' = proc() termorder_procs['mat'](args) end proc, 'mat' = proc(weights, vars::list(name)) option remember; local SUBS_SET, VARS, WVARS, W, w, i, j, a, b; if type(weights, 'listlist') then W := weights; if {nops(vars)} <> {seq(nops(i), i=W)} then error "the number of weights does not match the number of variables" end if; elif type(weights, 'Matrix') then a,b := op(1,weights); if b <> nops(vars) then error "the number of weights does not match the number of variables" end if; W := [seq([seq(weights[i,j], j=1..b)], i=1..a)] elif type(weights, 'matrix') then WARNING("this type of matrix is incompatible with remember tables, please use a Matrix or a list of lists instead"); W := weights; a,b := seq(op(2,i),i=op(1..2,op(W))); if b <> nops(vars) then error "the number of weights does not match the number of variables" end if; W := [seq([seq(W[i,j], j=1..b)], i=1..a)] else error "first argument must be a list of lists or a Matrix" end if; if not type(W, 'listlist'('rational')) then error "the matrix must have rational entries" else W := [seq(lcm(seq(denom(i),i=W[j]))*W[j], j=1..nops(W))] end if; subs({SUBS_SET={seq(vars[i]=vars[i]*mul(w[j]^W[j][i], j=1..nops(W)), i=1..nops(vars))}, VARS=vars, WVARS=[seq(w[j],j=1..nops(W))]}, proc(f) local lm; tcoeff(lcoeff(subs(SUBS_SET, f), WVARS), VARS, 'lm'), lm end proc) end proc ]); # set of all supported term order names termorder_types := map(op,{indices(termorder_procs)}); # returns a procedure for evaluating the leading term of a polynomial # the leading term is returned as the sequence (LEADING_COEFFCIENT, LEADING_MONOMIAL) evaluate_termorder := proc(tord::function) option system, remember; local lt; lt := termorder_procs[op(0,tord)](op(tord)); if not type(lt, 'procedure') then error "invalid monomial order" else eval(lt) end if end proc; # characterizes a termorder's elimination properties like this # prod(grlex(x,y), plex(z,w,t), tdeg(a,b)) -> [{x,y}, {z}, {w}, {t}, {a,b}] tord_split := proc(f::function) option system, remember; local i; if op(0,f)='plex' then [seq({i}, i=op(f))] elif op(0,f)='lexdeg' then [seq({op(i)}, i=[op(f)])] elif op(0,f)='prod' then [seq(op(procname(i)), i=[op(f)])] else [indets(f,'name')] end if end proc; # returns true if U \ X >> X for the term order is_elim_order := proc(f::function, X::set(name)) local Xt, L, v; L := ListTools:-Reverse(tord_split(f)); Xt := {}; for v in L do Xt := Xt union v; if not Xt subset X then return false elif Xt = X then return true end if; end do; false end proc; # remove variables from a term order, keeps X elim_order := proc(tord::function, X::set(name)) local T, i; if member(op(0,tord), {'plex', 'grlex', 'tdeg'}) then T := subs({seq(i=NULL, i=indets(tord, 'name') minus X)}, [op(tord)]); if nops(T)=0 then NULL elif nops(T)=1 then 'plex'(op(T)) else op(0,tord)(op(T)) end if; elif op(0, tord)='lexdeg' then T := [op(subs({seq(i=NULL, i=indets(tord, 'name') minus X)}, tord))]; T := remove(member, T, {[]}); if nops(T)=0 then NULL elif nops(T)=1 then 'tdeg'(op(T[1])) else 'lexdeg'(op(T)) end if elif op(0, tord)='prod' then T := map(procname, [op(tord)], X); if nops(T)=0 then NULL elif nops(T)=1 then T[1] else 'prod'(op(T)) end if else NULL end if; end proc; # return the term orders most similar to the given order # there is a real algorithm (Bayer ?) for comparing term orders using matrices # if anyone wants to try, you can use gbwalk_algorithm:-tord2mat similar_orders := proc(tordlist::{list(function), set(function)}, tord::function) local L, tsplit, vars, i, val; vars := indets(tord, 'name'); L := [op(select(proc(a) evalb(indets(a, 'name')=vars) end proc, tordlist))]; if member(tord, L) then return [tord] elif nops(L) <= 1 then return L end if; tsplit := tord_split(tord); L := map(proc(a) [a, mul(2^`if`(is_elim_order(a, `union`(op(tsplit[i..-1]))), i, 0), i=1..nops(tsplit))] end proc, L); val := max(seq(i[2], i=L)); [seq(`if`(i[2]=val, i[1], NULL), i=L)]; end proc; ## ## Next Monomial Procedures # SPEED CRITICAL (for the resulting procedures) # not for FGLM, but if you want to compute a basis for # the quotient ring as a vector space... # the structure is the same as for term orders, except these procedures # compute the next monomial not divisible by a set of monomials # if no such monomial exists, they return FAIL # evaluate_nextmon(plex(x,y,z)) calls nextmon_procs[plex](x,y,z) # the indicies of this table are assumed to be the term orders # for which a next monomial can be computed # matrix and wdeg term orders are not supported # because I can't figure out how to do it efficiently nextmon_procs := table([ 'plex' = proc() option remember; local VARS; if nargs=0 then proc() FAIL end proc elif type([args],'list'('name')) then subs('VARS'=[args], proc(m1, M) local m, i; if member(1,M) then return FAIL end if; m := m1*VARS[-1]; i := 1; while member(true, map2(divide, m, M)) do m := `if`(i < nops(VARS), lcoeff(m, VARS[-i..-1])*VARS[-i-1], FAIL); i := i+1; end do; m end proc ) else error "sequence of variable names expected" end if end proc, 'grlex' = proc() option remember; local VARS, NOPSVARS; if nargs=0 then proc() FAIL end proc elif type([args],'list'('name')) then subs({'VARS'=[args], 'NOPSVARS'=nargs}, proc(m::monomial, M::set(monomial)) local d, q, s, s1, i, R, max_deg; if member(1,M) then return FAIL end if; R := map(proc(a) `if`(nops(indets(a))=1 and member(op(indets(a)), VARS), a, NULL) end proc, M); max_deg := `if`(indets(R)=indets(VARS), add(degree(i)-1, i=R), infinity); s1 := m; while (member(true, map2(divide, s1, M)) or s1=m) and degree(s1) <= max_deg do d := degree(s1); s := VARS[-1]^(d+1); for i from 1 to NOPSVARS-1 do if divide(s1, VARS[-i]) then s := lcoeff(s1, VARS[-i..-1], 'q'); s := VARS[-i-1]*s*VARS[-1]^(degree(q)-1); break; end if; end do; s1 := s; end do; `if`(degree(s1) > max_deg, FAIL, s1) end proc ) else error "sequence of variable names expected" end if end proc, 'tdeg' = proc() option remember; local VARS, SRAV, NOPSVARS, i; if nargs=0 then proc() FAIL end proc elif type([args],'list'('name')) then subs({VARS=[args], SRAV=[seq(args[-i], i=1..nargs)], NOPSVARS=nargs}, proc(m::monomial, M::set(monomial)) local d, d2, s, s1, i, R, max_deg; if member(1,M) then return FAIL end if; R := map(proc(a) `if`(nops(indets(a))=1 and member(op(indets(a)), VARS), a, NULL) end proc, M); max_deg := `if`(indets(R)=indets(VARS), add(degree(i)-1, i=R), infinity); s1 := m; while (member(true, map2(divide, s1, M)) or s1=m) and degree(s1) <= max_deg do d := degree(s1); s := VARS[-1]^(d+1); for i from 1 to NOPSVARS-1 do if divide(s1, VARS[-i]) then d2 := degree(s1, VARS[1]); s := VARS[-i-1]^(d2+1)*s1/(VARS[-i]*VARS[1]^d2); end if; end do; s1 := s; end do; `if`(degree(s1) > max_deg, FAIL, s1) end proc ) else error "sequence of variable names expected" end if end proc, 'lexdeg' = proc() option remember; local i; eval(evaluate_nextmon('prod'(seq(`if`(nops(i) > 0, 'tdeg'(op(i)), NULL), i=args)))) end proc, 'prod' = proc() option remember; local nm, vars, i; if nargs=0 then return proc(m,M) FAIL end proc else nm := [seq([i,indets(i,'name')], i=[args])]; nm := remove(proc(a) evalb(nops(a[2])=0) end proc, nm); vars := map(a->a[2], nm); nm := map(a->evaluate_nextmon(a[1]), nm); end if; proc(m, M) local lc, lm, i, m1, r; if member(1, M) then return FAIL end if; m1 := m; r := FAIL; i := nops(vars); while i > 0 and (r=FAIL or member(true, map2(divide, m1, M))) do lc := lcoeff(m1, `union`(op(i..-1,vars)), 'lm'); lcoeff(lm, vars[i], 'lm'); r := lc*nm[i](lm, M); if r=FAIL then i := i-1 else m1 := r; i := nops(vars) end if; end do; r end proc end proc ]); # returns a procedure for computing the next monomial # which is not divisible by a given set of monomials # or FAIL if no such monomial exists evaluate_nextmon := proc(tord::function) option system, remember; local nm; nm := nextmon_procs[op(0,tord)](op(tord)); if not type(nm, 'procedure') then error "a next monomial procedure is not available for this monomial order" else eval(nm) end if end proc; ## ## Rootofs substitution # returns substitutions for all RootOfs and radical expressions # the idea is to convert each RootOf or radical into a variable # and add its minimal polynomial to the generating set # we perform computations using a product term order in which the RootOf # and radical variables are strictly less than the original variables # this procedure returns a sequence of two sets and a list; # a set of forward substitutions, a set of reverse substitutions, # and a list of minimal polynomials # if a RootOf and a radical expression are equivalent, the # reverse substitution will return the radical expression substitute_rootofs := proc(f) local Rofs, MinimalPolys, ForwardRadicalSubs, ForwardRootOfSubs, ReverseRadicalSubs, _z, i, fsubs, rsubs; global _Z; ForwardRadicalSubs := {`if`(has(f,interface('imaginaryunit')), interface('imaginaryunit')=RootOf(_z^2+1,'index'=1), NULL), seq(i=convert(i,'RootOf'), i=indets(f,'radical'))}; ReverseRadicalSubs := map(proc(a) a=convert(a,'radical') end proc, indets(ForwardRadicalSubs, 'RootOf')); Rofs := sort([op(indets(f, 'RootOf') union map(lhs, ReverseRadicalSubs))], length); ForwardRootOfSubs := []; MinimalPolys := []; for i from 1 to nops(Rofs) do ForwardRootOfSubs := [op(ForwardRootOfSubs), Rofs[i]=cat(_z,i)]; MinimalPolys := [op(MinimalPolys), subs({'_Z'=cat(_z,i), op(ForwardRootOfSubs)}, op(1, Rofs[i]))]; end do; fsubs := {op(subs(ForwardRootOfSubs, ForwardRadicalSubs)), op(subs(ReverseRadicalSubs, ForwardRootOfSubs)), op(ForwardRootOfSubs)}; rsubs := {op(subs(ReverseRadicalSubs, [seq(rhs(i)=lhs(i), i=ForwardRootOfSubs)])), op(subs(ForwardRootOfSubs, ReverseRadicalSubs))}; fsubs, rsubs, MinimalPolys; end proc; collect_rootof_powers := proc(f) local F, R, M, T, rofsvars, rewrite_vars, m, lm, p, t, i, j, k, _s, d, sp; F, R, M := substitute_rootofs(f); rofsvars := indets(M, 'name') minus indets(f, 'name'); rewrite_vars := map(lhs, select(type, R, 'name'='name'^'rational')); T := table(); for m in M do p := lcoeff(m, rewrite_vars, 'lm')*lm; t := p-m; if assigned(T[t]) then T[t] := T[t] union {p} else T[t] := {p} end if; end do; j := 1; for i in [indices(T)] do t := op(i); p := T[t]; if type(t, 'name') and not member(t, indets(p,'name')) then # bugfix M := remove(has, M, indets(p,'name') intersect rewrite_vars); d := mul(degree(k),k=p); sp := t=cat(_s,j)^d; if not type(t,'name') or type(t, 'name'^'integer') then M := [op(M), cat(_s,j)^d - t]; end if; F := map(a->lhs(a)=subs(sp, rhs(a)), subs({seq(op(indets(k))=cat(_s,j)^(d/degree(k)), k=p)}, F)) union {sp}; R := remove(has, R, indets(p,'name')) union {cat(_s,j) = t^(1/d)}; j := j+1; end if; end do; F, R, subs(F, M), indets(subs(F,M),'name') intersect rofsvars; end proc; substitute_everything := proc(f) local C, F, R, M, V, FC, RC; global constants; C := {constants} intersect indets(f,'name'); if nops(C)=0 then return collect_rootof_powers(f) end if; FC := map(a->a=`tools/gensym`(a), C); RC := map(rhs=lhs, FC); F, R, M, V := collect_rootof_powers(subs(FC, f)); F := subs(RC, F); R := subs(RC, R); F, R, M, V; end proc; ## ## Normalform Procedures # SPEED CRITICAL # this table contains procedures for computing normal forms, # ie : dividing a polynomial by a list of polynomials. # there are optimized versions for specific cases normalform_procs := table([ # integer coefficients # basis format : [lc, lm, poly, sugar] # input : polynom, sugar, basis, leadterm_proc, characteristic, variables 'rat_int_cof' = proc(f, fs, G, lt) local p, r, rs, g, qm, qc, lc, lm; p := f; r := 0; rs := fs; while p <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := lc/g[1]; p := expand(denom(qc)*p - numer(qc)*qm*g[3]); r := denom(qc)*r; rs := max(rs, g[4]+degree(qm)); break end if; end do; if qm = 0 then r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; r := iprimpart(r); lc, lm := lt(r); `if`(signum(lc) = -1, [-lc, lm, -r, rs], [lc, lm, r, rs]) end proc, # integer coefficients - no scaling # basis format : [lc, lm, poly] # input : polynom, basis, leadterm_proc # used by FGLM and NormalForm # r/s is equal to the remainder with no scaling 'rat_int_noscale' = proc(f, G, lt) local p, r, s, g, qm, qc, lc, lm; p := expand(f); qc := icontent(f); p, s := op(`if`(qc=0 or qc=1, [p,1], [p/qc, 1/qc])); r := 0; while p <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := lc/g[1]; p := expand(denom(qc)*p - numer(qc)*qm*g[3]); r := denom(qc)*r; s := denom(qc)*s; break end if; end do; if qm = 0 then r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; qc := icontent(r); `if`(qc=0 or qc=1, [r,s], [r/qc, s/qc]) end proc, # faster version - maps over sums # has a bug 'rat_int_noscale2' = proc(f, G, lt) local p, p1, vars, REDUCE, qc; REDUCE := proc(p) local lc, lm, g, qm; lc := lcoeff(p, vars, 'lm'); for g in G do if divide(lm, g[2], 'qm') then return p-lc/g[1]*qm*g[3] end if; end do; procname(p) := p; p end proc; p := expand(f); p1 := 0; vars := indets(lt(mul(i,i=indets(p,'name')))[2], 'name'); while p <> p1 do p1 := p; p := `if`(type(p,`+`), expand(map(REDUCE, p)), REDUCE(p)); end do; qc := icontent(p); `if`(qc=0 or qc=1, [p,1], [p/qc, 1/qc]) end proc, # rational function coefficients # basis format : [lc, lm, poly, sugar] # input : polynom, sugar, basis, leadterm_proc, characteristic, variables 'rat_fcn_cof' = proc(f, fs, G, lt) local vars, p, r, rs, s, g, qm, qc, lc, lm; p := f; r := 0; rs := fs; vars := args[6]; while p <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := normal(lc/g[1]); p := expand(denom(qc)*p - numer(qc)*qm*g[3]); r, rs := denom(qc)*r, max(rs, g[4]+degree(qm)); break end if; end do; if qm = 0 then r := r + lc*lm; p := expand(p - lc*lm) end if; end do; if r=0 then [0, 1, 0, rs] else qc := content(r, vars); if qc <> 0 then r := normal(r/qc) end if; lc, lm := lt(r); s := signum(lcoeff(lc)); [s*lc, lm, s*r, rs] end if; end proc, # rational function coefficients - no scaling # basis format : [lc, lm, poly] # input : polynom, basis, leadterm_proc, variables # used by FGLM and NormalForm # r/s is equal to the remainder with no scaling 'rat_fcn_noscale' = proc(f, G, lt) local vars, p, r, s, g, qm, qc, lc, lm; vars := `if`(nargs >= 5, args[5], indets(f,'name') union indets(G,'name')); p := normal(f); qc := content(p, vars, 'p'); s := `if`(qc=0, 1, 1/qc); r := 0; while p <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := normal(lc/g[1]); p := normal(denom(qc)*p - numer(qc)*qm*g[3]); r := denom(qc)*r; s := denom(qc)*s; break end if; end do; if qm = 0 then r := r + lc*lm; p := normal(p - lc*lm) end if; end do; qc := content(r, vars, 'r'); `if`(qc=0 or qc=1, [r, normal(s)], [r, normal(s/qc)]) end proc, # faster version - maps over sums 'rat_fcn_noscalex' = proc(f, G, lt) local p, p1, vars, REDUCE, qc, d; REDUCE := proc(p) local lc, lm, g, qm; lc := lcoeff(p, vars, 'lm'); for g in G do if divide(lm, g[2], 'qm') then return p-lc/g[1]*qm*g[3] end if; end do; procname(p) := p; p end proc; p := normal(f, 'expanded'); p, d := numer(p), denom(p); p1 := 0; vars := indets(lt(mul(i,i=indets(p,'name')))[2], 'name'); while p <> p1 do p1 := p; p := `if`(type(p,`+`), normal(map(REDUCE, p), 'expanded'), normal(REDUCE(p))); p, d := numer(p), d*denom(p); end do; qc := content(p, vars, 'p'); `if`(qc=0 or qc=1, [p,1/d], [p, 1/d/qc]) end proc, # rational function coefficients, keeps track of divisions # polynomial format : [lc, lm, poly, sugar, LIST] : where LIST expresses the element in terms of the input # input : polynom, sugar, coefficient_vector, basis, leadterm_proc, characteristic, variables, lazy_evaluation 'rat_trk_cof' = proc(f, fs, fL, G, lt) local vars, p, r, rs, rl, s, qm, qc, lc, lm, i, GL, FL; vars := args[7]; p, r, rs := normal(f), 0, fs; GL := [seq(cat(`GL`,i), i=1..nops(G))]; rl := `FL`; while p <> 0 do lc, lm := lt(p); for i from 1 to nops(G) do if divide(lm, G[i][2], 'qm') then qc := normal(lc/G[i][1]); p := normal(denom(qc)*p - numer(qc)*qm*G[i][3]); rl := denom(qc)*rl - numer(qc)*qm*GL[i]; r, rs := denom(qc)*r, max(rs, G[i][4]+degree(qm)); break end if; end do; if i > nops(G) then r := r + lc*lm; p := normal(p - lc*lm) end if; end do; r := normal(r); # skip work for r=0 since result is discarded anyway if r=0 and nargs > 7 and args[8]=true then [0, 1, 0, rs, [seq(0,i=1..nops(fL))]] elif r=0 then [0, 1, 0, rs, map(a->normal(a), expand(subs({seq(GL[i]=G[i][5], i=1..nops(G)), `FL`=fL}, rl)))] else qc := content(r, vars); r := normal(r/qc); lc, lm := lt(r); s := signum(lcoeff(lc)); rl := collect(rl, {op(GL),`FL`}, 'distributed'); [s*lc, lm, s*r, rs, map(a->normal(s*a/qc), expand(subs({seq(GL[i]=G[i][5], i=1..nops(G)), `FL`=fL}, rl)))] end if end proc, # integer coefficients mod p # basis format : [lc, lm, poly, sugar] # input : polynom, sugar, basis, leadterm_proc, characteristic, variables 'mod_int_cof' = proc(f, fs, G, lt, char) local p, r, rs, g, qm, lc, lm; p := f; r := 0; rs := fs; while lc <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then p := Expand(p - lc/g[1]*qm*g[3]) mod char; rs := max(rs, g[4]+degree(qm)); break end if; end do; if qm = 0 then r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; r := Expand(r) mod char; lc := lt(r)[1]; if lc <> 0 then r := r/lc mod char; end if; [lt(r), r, rs] end proc, 'mod_int_noscale' = proc(f, G, lt, char) local p, r, g, qm, lc, lm; p := Expand(f) mod char; r := 0; while lc <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then p := Expand(p - lc/g[1]*qm*g[3]) mod char; break end if; end do; if qm = 0 then r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; [Expand(r) mod char, 1] end proc, # faster version - maps over sums # has a bug 'mod_int_noscale2' = proc(f, G, lt, char) local p, p1, vars, REDUCE, qc; REDUCE := proc(p) local lc, lm, g, qm; lc := lcoeff(p, vars, 'lm'); for g in G do if divide(lm, g[2], 'qm') then return p-lc/g[1]*qm*g[3] end if; end do; procname(p) := p; p end proc; p := expand(f) mod char; p1 := 0; vars := indets(lt(mul(i,i=indets(p,'name')))[2], 'name'); while p <> p1 do p1 := p; p := `if`(type(p,`+`), expand(map(REDUCE, p)) mod char, REDUCE(p) mod char); end do; [expand(p) mod char, 1] end proc, # rational function coefficients mod p # basis format : [lc, lm, poly, sugar] # input : polynom, sugar, basis, leadterm_proc, characteristic, variables 'mod_fcn_cof' = proc(f, fs, G, lt, char, vars) local p, r, rs, g, c, qm, qc, lc, lm; p := f; r := 0; rs := fs; while lc <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := Normal(lc/g[1]) mod char; p := Normal(denom(qc)*p - numer(qc)*qm*g[3]) mod char; r, rs := denom(qc)*r, max(rs, g[4]+degree(qm)); break end if; end do; if qm = 0 then r := r + lc*lm; p := Normal(p - lc*lm) mod char end if; end do; r := Primpart(r, vars) mod char; c := lcoeff(lt(r)[1]); if c <> 0 then r := r/c mod char; end if; [lt(r), r, rs] end proc, # reduce polynomial with rational function coefficients mod p, with no scaling # basis format : [lc, lm, poly, sugar] # input : polynom, sugar, basis, leadterm_proc, characteristic, variables 'mod_fcn_noscale' = proc(f, G, lt, char, vars) local p, r, s, g, qm, qc, lc, lm; p := Normal(f) mod char; p, qc := numer(p), denom(p); qc := (Content(p, vars, 'p') mod char)/qc; s := `if`(qc=0, 1, 1/qc); r := 0; while p <> 0 do lc, lm := lt(p); qm := 0; for g in G do if divide(lm, g[2], 'qm') then qc := Normal(lc/g[1]) mod char; p := Normal(denom(qc)*p - numer(qc)*qm*g[3]) mod char; r := denom(qc)*r; s := denom(qc)*s; break end if; end do; if qm = 0 then r := r + lc*lm; p := Normal(p - lc*lm) mod char end if; end do; r := Normal(r) mod char; qc, r := 1/denom(r), numer(r); qc := qc*Content(r, vars, 'r') mod char; `if`(qc=0 or qc=1, [r, Normal(s) mod char], [r, Normal(s/qc) mod char]) end proc, # rational function coefficients mod p, keeps track of divisions # polynomial format : [lt, lm, poly, sugar, LIST] : where LIST expresses the element in terms of the input # input : polynom, sugar, coefficient_vector, basis, leadterm_proc, characteristic, variables, lazy_evaluation 'mod_trk_cof' = proc(f, fs, fL, G, lt, char, vars) local p, r, rs, rl, qm, qc, lc, lm, s, i, GL, FL; p := f; r := 0; rs := fs; GL := [seq(cat(`GL`,i), i=1..nops(G))]; rl := `FL`; lc, lm := lt(p); while lc <> 0 do for i from 1 to nops(G) do if divide(lm, G[i][2], 'qm') then qc := Normal(lc/G[i][1]) mod char; p := Normal(denom(qc)*p - numer(qc)*qm*G[i][3]) mod char; rl := denom(qc)*rl - numer(qc)*qm*GL[i]; r, rs := denom(qc)*r, max(rs, G[i][4]+degree(qm)); lc, lm := lt(p); break end if; end do; if i > nops(G) then r := r + lc*lm; p := Normal(p - lc*lm) mod char; lc, lm := lt(p) end if; end do; # skip work for r=0 since result is discarded anyway if r=0 and nargs > 7 and args[8]=true then [0, 1, 0, rs, [seq(0,i=1..nops(fL))]] elif r=0 then [0, 1, 0, rs, map(a->Normal(a) mod char, expand(subs({seq(GL[i]=G[i][5], i=1..nops(G)), `FL`=fL}, rl)))] else qc := Content(r, vars) mod char; r := Normal(r/qc) mod char; lc, lm := lt(r); s := signum(lcoeff(lc)); rl := collect(rl, {op(GL),`FL`}, 'distributed'); [s*lc, lm, s*r, rs, map(a->Normal(s*a/qc) mod char, expand(subs({seq(GL[i]=G[i][5], i=1..nops(G)), `FL`=fL}, rl)))] end if end proc ]); ## ## S-polynomial procedures # just like for normal forms, we have separate procedures # for computing syzygy polynomials. spolynomial_procs := table([ 'rat_int_cof' = proc(f, g) ## compute S-polynomial with integer coefficients option inline; denom(f[1]/g[1])*expand(lcm(f[2],g[2])/f[2]*f[3]) -numer(f[1]/g[1])*expand(lcm(f[2],g[2])/g[2]*g[3]), max(f[4] + degree(lcm(f[2],g[2])/f[2]), g[4] + degree(lcm(f[2],g[2])/g[2])) end proc, 'rat_fcn_cof' = proc(f, g) ## compute S-polynomial with rational function coefficients local mf, mg, c; gcd(f[2],g[2],'mg','mf'); c := normal(f[1]/g[1]); expand(denom(c)*mf*f[3] - numer(c)*mg*g[3]), max(f[4] + degree(mf), g[4] + degree(mg)) end proc, 'rat_trk_cof' = proc(f, g) ## compute S-polynomial with rational function coefficients ## keeps track of divisions, polynomial format : [lt, lm, poly, LIST] local mf, mg, c; gcd(f[2],g[2],'mg','mf'); c := normal(f[1]/g[1]); expand(denom(c)*mf*f[3] - numer(c)*mg*g[3]), max(f[4] + degree(mf), g[4] + degree(mg)), expand(denom(c)*mf*f[5] - numer(c)*mg*g[5]) end proc, 'mod_int_cof' = proc(f, g, char) ## compute S-polynomial with integer coefficients mod p option inline; denom(f[1]/g[1])*Expand(lcm(f[2],g[2])/f[2]*f[3]) -numer(f[1]/g[1])*Expand(lcm(f[2],g[2])/g[2]*g[3]) mod char, max(f[4] + degree(lcm(f[2],g[2])/f[2]), g[4] + degree(lcm(f[2],g[2])/g[2])) end proc, 'mod_fcn_cof' = proc(f, g, char) ## compute S-polynomial with rational function coefficients mod p local mf, mg; gcd(f[2],g[2],'mg','mf'); Expand(g[1]*mf*f[3] - f[1]*mg*g[3]) mod char, max(f[4] + degree(mf), g[4] + degree(mg)) end proc, 'mod_trk_cof' = proc(f, g, char) ## compute S-polynomial with rational function coefficients mod p ## keeps track of divisions, polynomial format : [lt, lm, poly, LIST] local mf, mg; gcd(f[2],g[2],'mg','mf'); Expand(g[1]*mf*f[3] - f[1]*mg*g[3]) mod char, max(f[4] + degree(mf), g[4] + degree(mg)), expand(g[1]*mf*f[5] - f[1]*mg*g[5]) end proc ]); ## ## Groebner Basis Algorithm # The method implemented is the Gebauer & Moller variant of the Buchberger algorithm # the selection strategy is the normal strategy with ties broken by lowest sugar. # there are no apologies for this, the strategy is optimized for tdeg. # in most cases lexicographic and elimination orders should be handled by conversion. # furthermore, the strategy is coded to use ONE termorder evaluation per pair selected, # which is VERY VERY good if your termorder is costly to evaluate. buchberger_algorithm := module() local preprocess, update, G, P, lt, ltdeg, vars, char, coftype; export gbasis, inter_reduce, inter_reduce_gb, STRATEGY, KEEP_BASIS_REDUCED, DISCARD_CRITERION1, DISCARD_CRITERION2, DISCARD_REDUCTION_PAIRS, MAKE_ALL_PAIRS, PRE_REDUCE, POST_REDUCE; # default options for Groebner basis computations, these are all reset by BuchbergerCF. STRATEGY := 'normal_sugar'; # use the normal selection strategy, break ties with sugar KEEP_BASIS_REDUCED := true; # remove redundant polynomials from the partial basis DISCARD_CRITERION1 := true; # remove pair s(f,g) if lt(f) and lt(g) are relatively prime DISCARD_CRITERION2 := true; # remove pair s(f,g) if there exist h with lt(h) dividing lt(s(f,g)) # and s(f,h) and s(g,h) have already been considered. DISCARD_REDUCTION_PAIRS := false; # discard "reduction critical pairs" from a paper by Francois Boulier, off by default. MAKE_ALL_PAIRS := false; # make syzygies with ALL polynomials, including those already removed from the basis PRE_REDUCE := true; # inter-reduce the input polynomials before running Buchberger's algorithm POST_REDUCE := true; # inter-reduce after, obtaining a reduced Groebner basis # WARNING: DON'T TURN THIS OFF ! Some programs expect reduced bases ! # inter-reduce if desired, otherwise just sort the input polynomials preprocess := proc(gens) if PRE_REDUCE then inter_reduce_gb(gens, lt, coftype, vars, char) else sort(gens, proc(a,b) local lm; lm := lt(a[2]+b[2])[2]; evalb((lm <> b[2]) or (lm = a[2] and length(a[3]) > length(b[3]))) end proc) end if end proc; # program for inter-reducing polynomials. # maintains an ascending sorted list of polynomials for maximum efficiency inter_reduce := proc(gens, lt, coftype, vars, char) local F, i, j, r, sp; sp := proc(a,b) option remember; evalb( lt(a[2]+b[2])[2] <> a[2]) end proc; F := sort(gens, sp); i := 1; while i <= nops(F) do r := F[i]; F := [seq(F[j], j=1..i-1), seq(F[j], j=i+1..nops(F))]; r := normalform_procs[coftype](seq(r[j], j=3..nops(r)), F, lt, char, vars, true); if r[1] <> 0 then i := 1; while i <= nops(F) and sp(F[i], r) do i := i+1; end do; F := [seq(F[j], j=1..i-1), r, seq(F[j], j=i..nops(F))]; i := i+1 end if; end do; F end proc; # program for inter-reducing a Groebner basis. # we can do this operation in one pass on a sorted list of polynomials # this can produce significant savings in time # this is also a better pre-processor for Buchberger's algorithm # assumes no leading monomial is divisible by another (ie: minimal Groebner basis) # safe for Gebauer & Moller and for the Groebner walk # SPEED CRITICAL (especially for Groebner Walk) inter_reduce_gb := proc(gens, lt, coftype, vars, char) local F, R, i, j, k; F := sort(gens, proc(a,b) evalb(lt(a[2]+b[2])[2] = b[2]) end proc); R := table(); j := 1; for i from 1 to nops(F) do R[j] := normalform_procs[coftype](op(3..-1, F[i]), [seq(R[k], k=1..j-1)], lt, char, vars, true); j := `if`(R[j][1]=0, j, j+1); end do; [seq(R[i], i=1..j-1)] end proc; # Buchberger's algorithm # a pretty straightforward implementation of the Gebauer & Moller's installation # the work of keeping track of pairs and all is done in the update module below gbasis := proc(generators, ltproc, coefficienttype, variables, characteristic, zap_pairs:=[]) local G1, r, s, i, j, k, st, stats; global NoName; #stats : [# of zero-red, zero-red time, # other red, other reduction time, total time] stats := array(1..5, [0,0,0,0,time()]); userinfo(2, {'_buchberger', 'buchberger'}, 'NoName', `-> Buchberger algorithm`, `if`(member(coefficienttype,{'rat_trk_cof', 'mod_trk_cof'}), `(extended)`, ``)); update:-initialize(); lt := ltproc; vars := variables; char := characteristic; coftype := coefficienttype; # this is used ONLY by the degree selection strategy ltdeg := evaluate_termorder('grlex'(op(vars))); if not type(update:-select_pair[STRATEGY], 'procedure') then error "invalid strategy", STRATEGY end if; userinfo(3, {'_buchberger', 'buchberger'}, 'NoName', `if`(characteristic=0, ` rational coefficients`, ` modular coefficients`)); userinfo(3, {'_buchberger', 'buchberger'}, 'NoName', ` strategy:`,STRATEGY); userinfo(4, {'_buchberger', 'buchberger'}, 'NoName', ` preprocessing generators`); r := `if`(zap_pairs=[], preprocess(generators), generators); userinfo(4, {'_buchberger', 'buchberger'}, 'NoName', ` adding generators`); for i in r do G1 := update:-add_basis_element(i); end do; # remove pair (i,j) if {i,j} is a subset of any zap_pairs[k] for k in zap_pairs do update:-zap_pairs(k) end do; userinfo(4, {'_buchberger', 'buchberger'}, 'NoName', ` processing syzygies`); s := update:-select_pair[STRATEGY](); while s <> FAIL do userinfo(5, {'_buchberger', 'buchberger'}, 'NoName', ` current pair`,[s[1..2]]); st := time(); r := normalform_procs[coftype]( spolynomial_procs[coftype](G[s[1]],G[s[2]], char), G1, lt, char, vars, true); st := time()-st; update:-remove_pair(s); if r[1] = 0 then stats[1], stats[2] := stats[1]+1, stats[2]+st; userinfo(5, {'_buchberger', 'buchberger'}, 'NoName', ` zero reduction`,st,`sec`) else stats[3], stats[4] := stats[3]+1, stats[4]+st; userinfo(5, {'_buchberger', 'buchberger'}, 'NoName', ` new polynomial`,st,`sec`); userinfo(5, {'_buchberger', 'buchberger'}, 'NoName', ` lead monomial`,r[2]); G1 := update:-add_basis_element(r) end if; s := update:-select_pair[STRATEGY](); end do; userinfo(4, {'_buchberger', 'buchberger'}, 'NoName', ` postprocessing generators`); G1 := `if`(POST_REDUCE, inter_reduce_gb(update:-get_basis(), lt, coftype, vars, char), update:-get_basis()); stats[5] := time()-stats[5]; userinfo(2,{'_buchberger', 'buchberger'}, 'NoName',` total reductions :`,stats[3]+stats[1]); userinfo(3,{'_buchberger', 'buchberger'}, 'NoName',` zero reductions :`,stats[1]); userinfo(3,{'_buchberger', 'buchberger'}, 'NoName',` time wasted (sec)`,stats[2]); userinfo(2,{'_buchberger', 'buchberger'}, 'NoName',` total time (sec)`,stats[5]); userinfo(3,{'_buchberger', 'buchberger'}, 'NoName', `-----------------------------`); G1 end proc; # this module does all of the work # the pairs and everything are kept separate to simplify things update := module() local NP, B, n, f, a, b, c; export initialize, add_basis_element, remove_pair, zap_pairs, select_pair, get_basis; # initialize the internal variables # P - pair table, indexed by [i,j], stores the leading monomial of Spoly(i,j) # G - table of all basis elements found, in the order they were found # B - set of integers denoting which elements of G are in the current basis # n - index for the table G, holds nops(G). initialize := proc() P := table('symmetric'); G := table(); B := []; n := 0; f := 0; NULL end proc; # adds a polynomial r to the basis, and updates all pairs add_basis_element := proc(r) local i, j, p, q, B2; global NoName; n := n+1; G[n] := r; # constructs new pairs involving r, note : uses ALL of G NP := [seq([i, n, lcm(G[i][2], r[2]), evalb(lcm(G[i][2], r[2])=G[i][2]*r[2])], i=1..n-1)]; B2 := B; # remove basis elements which r can eliminate : lm(r) | lm(G[i]) B := [seq(`if`(divide(G[i][2], r[2]), NULL, i), i=B)]; # remove any old pairs which are rendered superfluous by a new pair involving r # doesn't clash with the MAKE_ALL_PAIRS option below (think about the leading monomials) for p in indices(P) do i,j := op(p); if DISCARD_CRITERION2 and divide(P[i,j], NP[i][3], 'q') and q <> 1 and divide(P[i,j], NP[j][3], 'q') and q <> 1 then P[i,j] := evaln(P[i,j]) elif DISCARD_REDUCTION_PAIRS and not (member(i,B) or member(j,B) or divide(G[j][2],G[i][2])) then P[i,j] := evaln(P[i,j]) end if; end do; # sort the new pairs, removing extra ones if desired NP := sort(`if`(MAKE_ALL_PAIRS, NP, [seq(NP[i],i=B2)]), proc(a,b) option inline; evalb(degree(a[3]) < degree(b[3])) end proc); # apply criterion to eliminate some of the new pairs i := 1; while i <= nops(NP) do if not NP[i][4] and DISCARD_CRITERION2 then q := true; for j from 1 to i-1 do if divide(NP[i][3], NP[j][3]) then NP := subsop(i=NULL, NP); q := false; break; end if; end do; if q then i := i+1; end if else i := i+1 end if; end do; # apply criterion to select only useful new pairs for q in NP do if not (DISCARD_CRITERION1 and q[4]) and (op(0,G[q[1]][3])=`+` or op(0,G[q[2]][3])=`+`) then i,j := q[1],q[2]; P[i,j] := q[3] end if; end do; # add r to the basis, and rebuild the selection polynomial B := [op(B),n]; f := add(a^p[1]*b^p[2]*c^(max(G[p[1]][4]-degree(G[p[2]][2]), G[p[2]][4]-degree(G[p[1]][2])) + degree(P[op(p)]))/P[op(p)], p=indices(P)); userinfo(4, {'_buchberger', 'buchberger'}, 'NoName', ` basis / pairs`, nops(B),`if`(type(f,'`+`'),nops(f),`if`(f=0,0,1))); # return the basis, either all of G, or a top-reduced subset of G `if`(KEEP_BASIS_REDUCED, [seq(G[i], i=B)], [seq(G[i], i=1..n)]) end proc; # get the basis, always gets a top-reduced subset get_basis := proc() local i; [seq(G[i],i=B)] end proc; # consign a pair to oblivion # hopefully you got it first, using select_pair below. remove_pair := proc(i,j) if type(f, '`+`') then f := remove(proc(t) evalb(degree(t,a)=i and degree(t,b)=j) end proc, f) elif degree(f,a)=i and degree(f,b)=j then f := 0 end if; P[i,j] := evaln(P[i,j]) end proc; # remove pairs involving elements of S zap_pairs := proc(S) local i,j; if type(f, '`+`') then f := remove(proc(t) evalb({degree(t,a),degree(t,b)} subset S) end proc, f) elif {degree(f,a),degree(f,b)} subset S then f := 0 end if; for i in S do for j in S do P[i,j] := evaln(P[i,j]) end do; end do; end proc; # get the highest priority pair, according to the strategy # stores the pair list as a polynomial, so we can do ONE termorder evaluation per selection :) # returns the indices of each polynomial, and the sugar of the pair (used by normalform) select_pair := table([ # normal selection strategy 'normal' = proc() local t; t := lt(f)[1]; if t = 0 then return(FAIL) elif type(t,'`+`') then tcoeff(t,[b,a],'t') end if; degree(t,a), degree(t,b), degree(t,c) end proc, # normal selection strategy, breaks ties using lowest sugar 'normal_sugar' = proc() local t, sug; t := tcoeff(lt(f)[1], c, 'sug'); if t = 0 then return(FAIL) elif type(t,'`+`') then tcoeff(t,[b,a],'t') end if; degree(t,a), degree(t,b), degree(sug,c) end proc, # sugar strategy, breaks ties with the normal strategy 'sugar' = proc() local t, sug; t := lt(tcoeff(f, c, 'sug'))[1]; if t = 0 then return(FAIL) elif type(t,'`+`') then tcoeff(t,[b,a],'t') end if; degree(t,a), degree(t,b), degree(sug,c) end proc, # degree strategy, added for the Groebner Walk # used to get around a limitation in Maple # Maple can not compute leading term of a polynomial w.r.t a matrix # termorder, when the term order has large integers (300 digits) # and the polynomial contains negative powers of the variables # this is how pairs are selected 'degree' = proc() local t, sug; t := tcoeff(ltdeg(f)[1], c, 'sug'); if t=0 then return FAIL elif type(t,'`+`') then tcoeff(t,[b,a],'t') end if; degree(t,a), degree(t,b), degree(sug,c) end proc ]); end module; end module; # Buchberger's algorithm with minimal overhead # SPEED IS CRITICAL (for Groebner Walk) # not for export, since it doesn't handle RootOfs # you should really handle RootOfs at the top level anyway # F : list of polynomials # tord : term order, ie: tdeg(x,y,z) # char : characteristic # optional arguments # # track=set(posint) : set of indicies of F. We will express the # resulting Groebner basis in terms of these F[i] (in set) # ignoring contributions from any other F[j] (not in set) # ex : track={1,2,...,nops(F)} expresses the Groebner basis in terms of the input # while track={1} expresses the Groebner basis in terms of only F[1] # (useful for computing in quotient rings) # # strategy=name : selection strategy : 'normal', 'normal_sugar', 'sugar', 'degree' # # zap_pairs=list(set) : ignore syzygies (i,j) if {i,j} are common to any set # ex: to compute a GB of [op(G1),op(G2)] where G1 and G2 are Groebner bases # use zap_pairs=[{seq(i,i=1..nops(G1))}, {seq(i,i=nops(G1)+1..nops(G1)+nops(G2))}] Buchberger := proc(F1, tord, char) option system, remember; local F, gb, lt, vars, i, j, k, TRACK, opts; lt := evaluate_termorder(tord); vars := indets(tord,'name'); opts := table(['track'={}, 'strategy'='normal_sugar', 'zap_pairs'=[], op(select(type, [args], 'equation'))]); buchberger_algorithm:-STRATEGY := opts['strategy']; TRACK := opts['track']: if char=0 then if nops(TRACK) > 0 then F := normal(F1); gb := buchberger_algorithm:-gbasis([seq([lt(F[i]), F[i], degree(F[i], vars), `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'rat_trk_cof', vars, 0, opts['zap_pairs']); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := expand(F1); gb := buchberger_algorithm:-gbasis([seq([lt(i),i,degree(i)], i=F)], lt, 'rat_int_cof', indets(tord,'name'), 0, opts['zap_pairs']); [seq(i[3], i=gb)] else F := normal(F1); gb := buchberger_algorithm:-gbasis([seq([lt(i),i,degree(i,vars)], i=F)], lt, 'rat_fcn_cof', indets(tord,'name'), 0, opts['zap_pairs']); [seq(i[3], i=gb)] end if elif nops(TRACK) > 0 then F := Normal(F1) mod char; gb := buchberger_algorithm:-gbasis([seq([lt(F[i]), F[i], degree(F[i], vars), `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'mod_trk_cof', vars, char, opts['zap_pairs']); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := Expand(F1) mod char; gb := buchberger_algorithm:-gbasis([seq([lt(i),i,degree(i)], i=F mod char)], lt, 'mod_int_cof', indets(tord,'name'), char, opts['zap_pairs']); [seq(i[3], i=gb)] else F := Normal(F1) mod char; gb := buchberger_algorithm:-gbasis([seq([lt(i),i,degree(i,vars)], i=F mod char)], lt, 'mod_fcn_cof', indets(tord,'name'), char, opts['zap_pairs']); [seq(i[3], i=gb)] end if end proc; # fully inter-reduce a list of polynomials # note that Groebner:-inter_reduce does not actually do this # not used by anything, but may be useful InterReduce := proc(F1, tord, char) option system, remember; local F, gb, lt, vars, i, j, k, TRACK, opts; lt := evaluate_termorder(tord); vars := indets(tord,'name'); opts := table(['track'={}, op(select(type, [args], 'equation'))]); TRACK := opts['track']: if char=0 then if nops(TRACK) > 0 then F := normal(F1); gb := buchberger_algorithm:-inter_reduce([seq([lt(F[i]), F[i], 0, `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'rat_trk_cof', vars, 0); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := expand(F1); gb := buchberger_algorithm:-inter_reduce([seq([lt(i),i,0], i=F)], lt, 'rat_int_cof', indets(tord,'name'), 0); [seq(i[3], i=gb)] else F := normal(F1); gb := buchberger_algorithm:-inter_reduce([seq([lt(i),i,0], i=F)], lt, 'rat_fcn_cof', indets(tord,'name'), 0); [seq(i[3], i=gb)] end if elif nops(TRACK) > 0 then F := Normal(F1) mod char; gb := buchberger_algorithm:-inter_reduce([seq([lt(F[i]), F[i], 0, `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'mod_trk_cof', vars, char); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := Expand(F1) mod char; gb := buchberger_algorithm:-inter_reduce([seq([lt(i),i,0], i=F)], lt, 'rat_int_cof', indets(tord,'name'), char); [seq(i[3], i=gb)] else F := Normal(F1) mod char; gb := buchberger_algorithm:-inter_reduce([seq([lt(i),i,0], i=F)], lt, 'rat_fcn_cof', indets(tord,'name'), char); [seq(i[3], i=gb)] end if end proc; # /fast/ inter-reduction of a minimal Groebner basis # SPEED IS CRITICAL (for Groebner Walk) InterReduceGB := proc(F1, tord, char) option system, remember; local F, gb, lt, vars, i, j, k, TRACK, opts; lt := evaluate_termorder(tord); vars := indets(tord,'name'); opts := table(['track'={}, op(select(type, [args], 'equation'))]); TRACK := opts['track']: if char=0 then if nops(TRACK) > 0 then F := normal(F1); gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(F[i]), F[i], 0, `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'rat_trk_cof', vars, 0); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := expand(F1); gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(i),i,0], i=F)], lt, 'rat_int_cof', indets(tord,'name'), 0); [seq(i[3], i=gb)] else F := normal(F1); gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(i),i,0], i=F)], lt, 'rat_fcn_cof', indets(tord,'name'), 0); [seq(i[3], i=gb)] end if elif nops(TRACK) > 0 then F := Normal(F1) mod char; gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(F[i]), F[i], 0, `if`(member(i,TRACK,'k'), [seq(0, j=1..k-1), 1, seq(0, j=k+1..nops(TRACK))], [seq(0, j=1..nops(TRACK))])], i=1..nops(F))], lt, 'mod_trk_cof', vars, char); [seq(i[3], i=gb)], [seq(i[5], i=gb)] elif indets(F1,'name') subset vars then F := Expand(F1) mod char; gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(i),i,0], i=F)], lt, 'mod_int_cof', indets(tord,'name'), char); [seq(i[3], i=gb)] else F := Normal(F1) mod char; gb := buchberger_algorithm:-inter_reduce_gb([seq([lt(i),i,0], i=F)], lt, 'mod_fcn_cof', indets(tord,'name'), char); [seq(i[3], i=gb)] end if end proc; # almost fraction free rational reconstruction xratrecon := proc(f, vars, N) local d; d := iratrecon(tcoeff(f,vars), N); `if`(d=FAIL, FAIL, iratrecon(denom(d)*f, N)) end proc; # XLA library for multimodular linear algebra # all operations are performed in place when possible XLA := module() export n, m, dtype, char, images, primes, A, B, N, p, mode, sig, rank, C, Cmod, Cdone, R, Ri, Rj, Rcol, subspec, RowReduceREF, Initialize, NumberOfImages, AddImages, ReplaceImage, ChineseRemainder, RationalReconstruction, Transpose, IsIndependent, LUSolve, RowReduce, RREF, ReduceRow, Pivots, APivots, NewPivots; # init images, primes, etc Initialize := proc() local opts, i, N1, TIMER; opts := table(['datatype'='float'[8], 'numprimes'=-1, 'characteristic'=0, 'mode'="default", op(select(type, [args], 'equation'))]); dtype := opts['datatype']; mode := opts['mode']; char := opts['characteristic']; N := 0; if dtype='float'[8] then # 5 primes = 38 digits N1 := 5; p := 2^25-1; elif dtype='integer'[8] then # 4 primes = 39 digits N1 := 4; p := 2^32-1; elif dtype='integer'[4] then # 8 primes = 39 digits N1 := 8; p := 2^16-1; end if; A := select(type,[args],{'Matrix','Vector'}); if nops(A)=0 then n := select(type, [args], 'posint'); if nops(n) > 1 then n,m := n[1],n[2]; elif nops(n)=1 then n,m := n[1],n[1]; else error "can not determine the dimensions of the matrix"; end if; A := rtable(1..n,1..m,'subtype'='Matrix', 'datatype'='integer', 'order'='C_order'); else if type(A[1],'Matrix') then A := A[1]; else error "only matrices are supported at this time" end if; n,m := op(1, A); end if; if char=0 then N1 := `if`(opts['numprimes']=-1, N1, opts['numprimes']) elif char > p then error "characteristic %1 too large for datatype %2", char, dtype; elif not isprime(char) then error "characteristic must be prime" else N1 := 1; p := char+1; end if; # used in ReplaceImage so must be initialized first subspec := []; C, Cmod, Cdone := 0, 1, {}; R, Ri, Rj, Rcol := 0, 1, 1, true; sig, rank := -1, -1; NumberOfImages(N1); NULL end proc; NumberOfImages := proc(N1::nonnegint) local L, pvec, ref, r, s; if nargs=0 then return N elif N1 < N or N1=0 then # remove some images while N > N1 do primes[N] := evaln(primes[N]); images[N] := evaln(images[N]); N := N-1; end do elif mode="default" then # possibly add some images while N1 > N and p > 2 do N,p := N+1,prevprime(p); primes[N] := p; images[N] := LinearAlgebra:-Modular:-Mod(p, A, dtype); end do; if N1 > N then error "failed to add images: ran out of primes"; end if; elif mode="RowReduce" or mode="RREF" then # add images of row reduced A ref := evalb(mode="RREF"); N := N+1; while N1 >= N do try p := prevprime(p); catch: error "failed to add images: ran out of primes"; end try; primes[N] := p; try images[N] := LinearAlgebra:-Modular:-Mod(p, A, dtype); LinearAlgebra:-Modular:-RowReduce(p, images[N], n, m, m, 0, 0, 'r', 's', 0, ref); if r=rank and s=sig then N := N+1; elif r <= rank then ; else mode := "default"; NumberOfImages(0); NumberOfImages(N1); RowReduceREF(); return NULL end if; catch: end try; end do; N := N-1; elif mode="LUSolve" then # add images of a solution Bx=A N := N+1; pvec := rtable(1..n-1,subtype='Vector'['column']); while N1 >= N and p > 2 do p := prevprime(p); primes[N] := p; try images[N] := LinearAlgebra:-Modular:-Mod(p, A, dtype); L := LinearAlgebra:-Modular:-Mod(p, B, dtype); LinearAlgebra:-Modular:-LUDecomposition(p, L, pvec, 0); LinearAlgebra:-Modular:-LUApply(p, L, pvec, images[N]); N := N+1; catch: end try; end do; N := N-1; if N1 > N then error "failed to add images: ran out of primes"; end if; else error "unknown mode"; end if; N; end proc; AddImages := proc(k::integer) NumberOfImages(N+k) end proc; ReplaceImage := proc(i::posint) if p = 2 then error "can not replace image: ran out of primes" end if; p := prevprime(p); primes[i] := p; images[i] := LinearAlgebra:-Modular:-Mod(p, A, dtype); if member(i, Cdone) then # chinese remainder is no longer valid C, Cmod, Cdone := 0, 1, {}; R, Ri, Rj, Rcol := 0, 1, 1, true; end if; NULL; end proc; # C is matrix, Cmod is modulus, Cdone is finished images ChineseRemainder := proc() local i, Ctodo; if subspec <> [args] then subspec := [args]; if nops(select(type, subspec, posint)) > 0 then error "only matrices are supported, please use ranges in the subspec" end if; C, Cmod, Cdone := 0, 1, {}; R, Ri, Rj, Rcol := 0, 1, 1, true; end if; if subspec=[] then if Cmod=1 then C := LinearAlgebra:-Modular:-ChineseRemainder(0, 'Cmod', images[1], primes[1]); if N > 1 then LinearAlgebra:-Modular:-ChineseRemainder(C, 'Cmod', [seq(images[i],i=2..N)], [seq(primes[i], i=2..N)]); end if; Cdone := {seq(i,i=1..N)}; else Ctodo := {seq(i,i=1..N)} minus Cdone; if nops(Ctodo) > 0 then LinearAlgebra:-Modular:-ChineseRemainder(C, 'Cmod', [seq(images[i], i=Ctodo)], [seq(primes[i], i=Ctodo)]); end if; Cdone := {seq(i,i=1..N)}; end if; else if Cmod=1 then C := LinearAlgebra:-Modular:-ChineseRemainder(0, 'Cmod', LinearAlgebra:-Modular:-Copy(primes[1], images[1], op(subspec)), primes[1]); if N > 1 then LinearAlgebra:-Modular:-ChineseRemainder(C, 'Cmod', [seq(LinearAlgebra:-Modular:-Copy(primes[i], images[i], op(subspec)),i=2..N)], [seq(primes[i], i=2..N)]); end if; Cdone := {seq(i,i=1..N)}; else Ctodo := {seq(i,i=1..N)} minus Cdone; if nops(Ctodo) > 0 then LinearAlgebra:-Modular:-ChineseRemainder(C, 'Cmod', [seq(LinearAlgebra:-Modular:-Copy(primes[i], images[i], op(subspec)), i=Ctodo)], [seq(primes[i], i=Ctodo)]); end if; Cdone := {seq(i,i=1..N)}; end if; end if; NULL end proc; # returns true if successful RationalReconstruction := proc() local i, j, r, a, b; ChineseRemainder(args); a, b := op(1,C); if R=0 then R := rtable(1..a, 1..b, 'subtype'='Matrix'); end if; i, j := Ri, Rj; if Rcol then # go down each column r := iratrecon(C[i,j], Cmod); to a*b while r <> FAIL do R[i,j] := r; i,j := modp(i,a)+1, `if`(i=a,modp(j,b)+1,j); r := iratrecon(C[i,j], Cmod); end do; if r=FAIL then Ri,Rj := i,j; false; else true; end if; else # go across each row r := iratrecon(C[i,j], Cmod); to a*b while r <> FAIL do R[i,j] := r; j,i := modp(j,a)+1, `if`(j=b,modp(i,a)+1,i); r := iratrecon(C[i,j], Cmod); end do; if r=FAIL then Ri,Rj := i,j; false; else true; end if; end if; end proc; Transpose := proc() local i; if n=m then LinearAlgebra:-LA_Main:-Transpose(A,'inplace'=true, outputoptions = []); for i from 1 to N do LinearAlgebra:-Modular:-Transpose(primes[i], images[i], 'inplace') end do; Ri, Rj, Rcol := Rj, Ri, not Rcol; else A := LinearAlgebra:-LA_Main:-Transpose(A,'inplace'=false, outputoptions = []); for i from 1 to N do images[i] := LinearAlgebra:-Modular:-Transpose(primes[i], images[i]) end do; Ri, Rj, Rcol := Rj, Ri, not Rcol; end if; end proc: # test if a row is independent, the row is always added to A as row k1 IsIndependent := proc(V, k1) local i, k, bad_images, is_dep; if nargs > 1 then A[k1, 1..m] := V; k := k1; else k := V; end if; bad_images := {}; is_dep := true; for i from 1 to N do try LinearAlgebra:-Modular:-Mod(primes[i], A, k, 1..m, images[i], k, 1..m); if k > 1 then LinearAlgebra:-Modular:-RowReduce(primes[i], images[i], k, 1..m, m, 0, 0, 'r', 0, 0, false); end if; if is_dep and rtable_scanblock(images[i], [k, 1..m], 'HasNonZero') then is_dep := false; bad_images := {seq(j, j=1..i-1)}; end if catch: bad_images := bad_images union {i}; end try; end do; while nops(bad_images) > 0 do i := bad_images[1]; ReplaceImage(i); try LinearAlgebra:-Modular:-RowReduce(primes[i], images[i], 1..k, 1..m, m, 0, 0, 'r', 0, 0, false); if is_dep and rtable_scanblock(images[i], [k, 1..m], 'HasNonZero') then is_dep := false; bad_images := {seq(j,j=1..N)} minus {i}; else bad_images := bad_images minus {i}; end if; catch: end try: end do; not is_dep; end proc; RowReduceREF := proc() local sigtable, s, r, i, j, x, ref; ref := evalb(mode="RREF"); sigtable := table('sparse'); rank := -1; i := 1: while i <= N do # row reduce all the images try LinearAlgebra:-Modular:-RowReduce(primes[i], images[i], n, m, m, 0, 0, 'r', 's', 0, ref); if r = rank then sigtable[i] := s; i := i+1; elif r < rank then ReplaceImage(i) elif i=1 then rank := r; sigtable[i] := s; i := i+1; else rank := r; for j from 1 to i-1 do ReplaceImage(j) end do; i := 1; end if; catch: ReplaceImage(i) end try; end do; # now all of the images have the same maximal rank # determine which signature is most common and therefore most likely to be correct if nops({entries(sigtable)}) = 1 then sig := s; else x := evaln(x); sig := [member(maxnorm(f), [coeffs(add(x^sigtable[i],i=1..N),x,'s')],'j'), degree(s[j])][2]; for i to N do while sigtable[i] <> sig do ReplaceImage(i); try LinearAlgebra:-Modular:-RowReduce(primes[i], images[i], n, m, m, 0, 0, 'r', 's', 0, ref); if r = rank and s = sig then sigtable[i] := s; elif r <= rank then ; else return FAIL; end if; catch: end try; end do; end do; end if; NULL end proc: RowReduce := proc() mode := "RowReduce"; RowReduceREF(args); end proc: RREF := proc() mode := "RREF"; RowReduceREF(args); end proc: # returns the list of pivot columns corresponding to each row # assumes that the rows are sorted so that row i appears before # row j if and only the pivot of row i appears in a column # to the left (or the same) as that of row j Pivots := proc() local P, i, j; if N < 1 then error "no images to use" end if; P := table(); j := 1; for i from 1 to rank do while images[1][i,j]=0 do j := j+1 end do; P[i] := j end do; op(op(P)); end proc: APivots := proc() local P, i, j; P := table(); i,j := 1,1; while i <= n and j <= m do while j <= m and A[i,j]=0 do j := j+1 end do; if j <= m then P[i] := j; i := i+1; end if; end do; op(op(P)); end proc: # old_pivots should be a sparse table with old_pivots[j]=1 if j is an old pivot NewPivots := proc(old_pivots::table) local P, i, j; if N < 1 then error "no images to use" end if; P := table(); j := 1; for i from 1 to rank do while images[1][i,j]=0 do j := j+1 end do; if old_pivots[j]=0 then P[i] := j end if; end do; op(op(P)); end proc; # note the images of X in A1 X = B1 are being stored LUSolve := proc(A1, B1) B := A1; Initialize(B1,args[3..-1],'mode'="LUSolve"); NULL end proc: end module: # F4 Algorithm # supports both commutative and non-commutative domains # does not support parameters, ie: rational function coefficients f4_algorithm := module() local G, B, P, n, lt, vars, listvars, INITIAL_IMAGES, IMAGE_MULTIPLIER, KEEP_BASIS_REDUCED, DISCARD_CRITERION1, DISCARD_CRITERION2, PRE_REDUCE, POST_REDUCE, msort_proc, mgcd, update_pairs; export f4; INITIAL_IMAGES := 3; # number of images initially used for row reduction IMAGE_MULTIPLIER := 1.7; # multiplier for the number of images if rational reconstruction fails (takes ceiling) KEEP_BASIS_REDUCED := true; # remove redundant polynomials from the partial basis DISCARD_CRITERION1 := true; # remove pair s(f,g) if lt(f) and lt(g) are relatively prime DISCARD_CRITERION2 := true; # remove pair s(f,g) if there exist h with lt(h) dividing lt(s(f,g)) # and s(f,h) and s(g,h) have already been considered. PRE_REDUCE := true; # inter-reduce the input polynomials before running Buchberger's algorithm POST_REDUCE := true; # inter-reduce after, obtaining a reduced Groebner basis # WARNING: DON'T TURN THIS OFF ! Some programs expect reduced bases ! # return a procedure to sort a list of monomials into descending order msort_proc := proc(tord) local S, rlistvars, h; if nops(listvars) = 1 or op(0,tord)='plex' then proc() `if`(nargs > 1, [op(sort(`+`(args), listvars))], [args]) end proc; elif op(0,tord)='grlex' then proc() `if`(nargs > 1, [op(sort(`+`(args), listvars, 'tdeg'))], [args]) end proc; elif op(0,tord)='tdeg' then h := evaln(h); rlistvars := [h,op(ListTools:-Reverse(listvars))]; proc() local f; if nargs > 1 then f := `+`(args); f := [coeffs(subs([seq(i=i/h,i=listvars)], expand(h^degree(f,vars)*f)), h)]; [seq(`if`(type(i,`+`), op(sort(i, rlistvars, 'plex', 'ascending')), i), i=sort(f, (a,b)->evalb(degree(a) > degree(b))))] else [args] end if; end proc; else proc() sort([args], proc(a,b) evalb(lt(a+b)[2]=a) end proc) end proc; end if; end proc; # monomial gcd, equivalent but faster than evalb(gcd(m1,m2,'q')=1) mgcd := proc(m1, m2, q) local f, g; f := m1+m2; g := mul(i^ldegree(f,i),i=listvars); q := m1/g; evalb(g=1); end proc; # many of these improvements should be incorporated into the Buchberger code update_pairs := proc(f) local lc, lm, r, NP, zap, R, q, i, j, k, L, TIMER; n := n+1; G[n] := f; lc, lm, r := op(f); # eliminate old pairs using new pairs involving r NP := table(): zap := table(): R := table('sparse'): for k to nops(P) do i,j,L := op(1..3, P[k]); if R[j]=0 then NP[j] := [mgcd(G[j][2],lm,'q'), lm*q, degree(lm*q)]; R[j] := 1; end if; if divide(L, NP[j][2], 'q') and q <> 1 then if R[i]=0 then NP[i] := [mgcd(G[i][2],lm,'q'), lm*q, degree(lm*q)]; R[i] := 1; end if; if divide(L, NP[i][2], 'q') and q <> 1 then zap[k] := 1; end if; end if; end do; zap := [indices(zap)]; if DISCARD_CRITERION2 and nops(zap) > 0 then P := subsop(seq(op(i)=NULL, i=zap), P) end if; # construct new pairs and sort by degree of lcm # credit for the "sorting with attributes" trick goes to Joe Riel :) NP := [seq(`if`(R[i]=0, [i, n, mgcd(G[i][2],lm,'q'), lm*q, degree(lm*q)], [i, n, op(NP[i])]), i=B)]: NP := map(attributes, sort([seq(setattribute(SFloat(i[5],0),i),i=NP)], `<`)); # apply criterion 2 to eliminate some new pairs if DISCARD_CRITERION2 then zap := table(); for i to nops(NP) do for j to i-1 do if mdivide(NP[i][4], NP[j][4]) then zap[i] := 1; break end if; end do; end do; zap := [indices(zap)]; if nops(zap) > 0 then NP := subsop(seq(op(i)=NULL, i=zap), NP) end if; end if; # update list of pairs (unsorted) if DISCARD_CRITERION1 then P := [op(P), seq(`if`(i[3], NULL, subsop(3=NULL,i)), i=NP)]; else P := [op(P), seq(subsop(3=NULL,i), i=NP)]; end if; # remove basis elements which r can eliminate if KEEP_BASIS_REDUCED then B := [seq(`if`(mdivide(G[i][2],lm),NULL,i),i=B), n]; else B := [op(B), n]; end if; NULL; end proc; f4 := proc(F, tord) local opts, char, ore_case, F1, d, Pd, M, Mdone, A, an, sk, Apiv, c, m, q, A1, i, j, k, NewRows, r, rd, zap, S, msort, N, stats, TIMER, TIMER2; # initialize G := table(): n := 0: # nops(G) B := []: # list of i with G[i] in the current basis P := []: # list of pairs sorted by degree # [1] total time [2] update pairs [3] symbolic proc. [4] build matrix # [5] row reduce [6] chinese rem. [7] extract rows [8] rational reconstruction # [9] interreduce # [10] debug stats := array([time(), 0, 0, 0, 0, 0, 0, 0, 0, 0]); userinfo(2, {'f4', '_f4'}, 'NoName', `-> F4 algorithm`); if type(tord, MonomialOrder) then ore_case := true; lt := tord["leading_term"]; listvars := tord["order_indets"]; vars := {op(listvars)}; msort := proc() ListTools:-Reverse(sort([args], tord["order"])) end proc; char := tord["algebra"]["commutation"]["ground_ring"]["characteristic"]; else ore_case := false; lt := evaluate_termorder(tord); vars := indets(tord,'name'); listvars := sort([op(vars)], (a,b)->evalb(lt(a+b)[2]=a)); msort := msort_proc(tord); char := `if`(type(args[-1], posint), args[-1], 0); end if; if not indets(F,'name') subset vars then error "parameters are not supported" elif char > 2^25 then error "characteristic must be smaller then 2^25" elif char<>0 and not isprime(char) then error "characteristic must be zero or prime" end if; if PRE_REDUCE and not ore_case then TIMER := time(): F1 := buchberger_algorithm:-inter_reduce_gb([seq([lt(i),i,0], i=F)], lt, `if`(char=0,'rat_int_cof','mod_int_cof'), vars, char); stats[9] := stats[9] + time() - TIMER; TIMER := time(): for i in F1 do update_pairs(i[1..3]) end do; stats[2] := stats[2] + time() - TIMER: else TIMER := time(): for i in F do update_pairs([lt(i),i]) end do; stats[2] := stats[2] + time() - TIMER: end if; # format of P : [i, j, lcm, deg] while nops(P) > 0 do # select pairs of minimal degree d := min(seq(i[4],i=P)); Pd, P := selectremove(proc(a) evalb(a[4]=d) end proc, P); userinfo(4, {'f4', '_f4'}, 'NoName', cat(` basis: `, n, ` degree: `, d, ` pairs: `, nops(Pd), ` of `, nops(Pd)+nops(P))); # construct rows corresponding to syzygies (two rows per pair) # mark the leading monomials of these rows as done TIMER := time(): Mdone := table([seq(i[3]=1, i=Pd)], 'sparse'); A := table(): an := 1: for i in Pd do A[an] := [i[1], i[3]/G[i[1]][2]]; A[an+1] := [i[2], i[3]/G[i[2]][2]]; an := an+2; end do; an := an-1; Apiv := table(): # symbolic preprocessing sk := 0; # number of processed A[i] while sk < an do sk := sk+1; k, q := op(A[sk]); # row is q*G[k] Apiv[sk] := q*G[k][2]; TIMER2 := time(): # extract coefficients to form the row if ore_case then A[sk] := [coeffs(tord["algebra"]["product"](q, G[k][3]) , vars, 'm')], [m] else A[sk] := [coeffs(expand(q*G[k][3]), vars, 'm')], [m]; end if; stats[10] := stats[10] + time() - TIMER2; # add a row to the queue for each new monomial for i in [m] do if Mdone[i]=0 then Mdone[i] := 1; # mark the monomial as done for j in B do if divide(i,G[j][2],'q') and j <> i then # add q*G[j] an := an+1; A[an] := [j, q]; break; end if; end do; end if; end do; end do; stats[3] := stats[3] + time() - TIMER; TIMER := time(): # sort the monomials corresponding to matrix columns M := msort(seq(op(i),i=[indices(Mdone)])); Mdone := table([seq(M[i]=i, i=1..nops(M))]); # build matrix A1 := rtable(1..an, 1..nops(M), 'storage'='rectangular', 'subtype'='Matrix'); # sparse strategy: put sparse rows first in the matrix # credit for the "sorting with attributes" trick goes to Joe Riel :) S := map(attributes, sort([seq(setattribute(SFloat(nops(A[i][2]),0),i), i=1..an)], `<`)); k := 1: for i in S do A1[k, [seq(Mdone[j], j=A[i][2])]] := rtable(A[i][1]); k := k+1 end do; # make a table of the current pivots Apiv := table([seq(Mdone[op(i)]=1, i=[entries(Apiv)])], 'sparse'); stats[4] := stats[4] + time() - TIMER; N := INITIAL_IMAGES; userinfo(4, {'f4', '_f4'}, 'NoName', cat(` matrix `, an, ` x `, nops(M))); userinfo(5, {'f4', '_f4'}, 'NoName', cat(` images `, N)); # row reduce matrix mod p and identify rows with new pivots TIMER := time(): XLA:-Initialize(A1, 'numprimes'=N, 'characteristic'=char); XLA:-RREF(); NewRows := XLA:-NewPivots(Apiv); stats[5] := stats[5] + time() - TIMER; userinfo(4, {'f4', '_f4'}, 'NoName', cat(` new pivots `, nops(NewRows))); # extract new rows using chinese remaindering and rational reconstruction for k in NewRows do i,j := op(k); TIMER := time(): XLA:-ChineseRemainder(i..i, j..-1); stats[6] := stats[6] + time() - TIMER; TIMER := time(): r := add(rhs(o)*M[[lhs(o)][2]+j-1], o=rtable_elems(XLA:-C)); stats[7] := stats[7] + time() - TIMER; TIMER := time(): if char=0 then r := xratrecon(r, vars, XLA:-Cmod); end if; stats[8] := stats[8] + time() - TIMER; while r=FAIL do N := ceil(N*IMAGE_MULTIPLIER); userinfo(5, {'f4', '_f4'}, 'NoName', cat(` images `, N)); TIMER := time(): XLA:-NumberOfImages(N); stats[5] := stats[5] + time() - TIMER; TIMER := time(): XLA:-ChineseRemainder(i..i, j..-1); stats[6] := stats[6] + time() - TIMER; TIMER := time(): r := add(rhs(o)*M[[lhs(o)][2]+j-1], o=rtable_elems(XLA:-C)); stats[7] := stats[7] + time() - TIMER; TIMER := time(): r := xratrecon(r, vars, XLA:-Cmod); stats[8] := stats[8] + time() - TIMER; end do; TIMER := time(): if char=0 then r := iprimpart(r) end if; stats[8] := stats[8] + time() - TIMER; TIMER := time(): update_pairs([lt(r),r]); stats[2] := stats[2] + time() - TIMER; end do; end do; # end of main loop TIMER := time(): if POST_REDUCE then # sort the Groebner basis G := sort([seq(G[i], i=B)], (a,b)->evalb(lt(a[2]+b[2])[2]=b[2])); # remove redundant Groebner basis elements zap := table('sparse'); for i from 2 to nops(G) do for j to i-1 do if divide(G[i][2], G[j][2]) then zap[i] := 1; break end if; end do; end do; if ore_case then G := [seq(i[3], i=subsop(seq(op(i)=NULL, i=[indices(zap)]), G))]; G := Groebner[InterReduce](G, tord) else G := [seq([op(i),0], i=subsop(seq(op(i)=NULL, i=[indices(zap)]), G))]; G := [seq(j[3], j=buchberger_algorithm:-inter_reduce_gb(G, lt, `if`(char=0,'rat_int_cof','mod_int_cof'), vars, char))]; end if; else G := [seq(G[i][3], i=B)]; end if; stats[9] := stats[9] + time() - TIMER; # clear memory used by XLA XLA:-Initialize(1,1); # output crap # [1] total time [2] update pairs [3] symbolic proc. [4] build matrix # [5] row reduce [6] chinese rem. [7] extract rows [8] rational reconstruction stats[1] := time() - stats[1]; userinfo(3, {'f4', '_f4'}, `NoName`, ` update pairs (sec)`, evalf(stats[2],3)); userinfo(3, {'f4', '_f4'}, `NoName`, ` symbolic proc. (sec)`, evalf(stats[3],3)); userinfo(3, {'f4', '_f4'}, `NoName`, ` build matrix (sec)`, evalf(stats[4],3)); userinfo(3, {'f4', '_f4'}, `NoName`, ` row reduce (sec)`, evalf(stats[5],3)); if char=0 then userinfo(3, {'f4', '_f4'}, `NoName`, ` chinese remainder (sec)`, evalf(stats[6],3)); userinfo(3, {'f4', '_f4'}, `NoName`, ` extract rows (sec)`, evalf(stats[7],3)); else userinfo(3, {'f4', '_f4'}, `NoName`, ` extract rows (sec)`, evalf(stats[6]+stats[7],3)); end if; if char=0 then userinfo(3, {'f4', '_f4'}, `NoName`, ` rational reconst. (sec)`, evalf(stats[8],3)); end if; userinfo(3, {'f4', '_f4'}, `NoName`, ` interreduce (sec)`, evalf(stats[9],3)); userinfo(2, {'f4', '_f4'}, `NoName`, ` total time (sec)`, evalf(stats[1],3)); userinfo(3, {'f4', '_f4'}, `NoName`, ` --------------------------------`); G; end proc; end module: # construct the linear system A x = B arising from an FGLM order change # returns the sequence (A, B, Ma, Mb) where the columns of A are the independent normal forms, # the columns of B are dependent, and if A x = B then Mb - Ma.x is the desired Groebner basis fglm_system := module() local ModuleApply, basis_lookup, nf_cache, rat_nf, mod_nf, char; # compute the remainder of f divided by G and put the # coefficients into row i of A (assumed to be zero) # we employ a caching scheme to remember divisors, eliminating an inner loop # fglm uses previous normalforms, for example: nf(x^3) = nf(x * nf(x^2)) # so we encounter the same monomials over and over rat_nf := proc(f, G, lt, A, i:=NULL) local p, c, lc, lm, lmd, j, qm, qc, g; p := expand(f); c := icontent(p); if c = 1 then ; elif c = 0 then c := 1 else p := p/c end if; lc, lm := lt(p); while p <> 0 do j := nf_cache[lm]; # check the cache # cache says the monomial is not divisible by any G[j] if j = FAIL then A[i,basis_lookup[lm]] := lc*c; p := expand(p - lc*lm) # cache says monomial is divisible by G[i] elif j <> 0 then g := G[j]; qc := lc/g[1]; divide(lm, g[2], 'qm'); p := expand(denom(qc)*p - numer(qc)*qm*g[3]); c := c/denom(qc) # cache doesn't know - we have to check each element of G else qm := 0; for g in G do if divide(lm, g[2], 'qm') then member(g, G, 'j'); nf_cache[lm] := j; # cache divisor qc := lc/g[1]; p := expand(denom(qc)*p - numer(qc)*qm*g[3]); c := c/denom(qc); break end if; end do; # no divisor found if qm=0 then nf_cache[lm] := FAIL; # cache "no divisor" A[i,basis_lookup[lm]] := lc*c; p := expand(p - lc*lm); end if end if; lc, lm := lt(p); end do; NULL end proc; mod_nf := proc(f, G, lt, A, i:=NULL) local p, lc, lm, lmd, j, qm, g; p := expand(f) mod char; lc, lm := lt(p); while p <> 0 do j := nf_cache[lm]; # check the cache # cache says the monomial is not divisible by any G[j] if j = FAIL then A[i,basis_lookup[lm]] := lc; p := expand(p - lc*lm) # cache says monomial is divisible by G[i] elif j <> 0 then divide(lm, G[j][2], 'qm'); p := expand(p - lc/G[j][1]*qm*G[j][3]) mod char; # cache doesn't know - we have to check each element of G else qm := 0; for g in G do if divide(lm, g[2], 'qm') then member(g, G, 'j'); nf_cache[lm] := j; # cache divisor p := expand(p - lc/g[1]*qm*g[3]) mod char; break end if; end do; # no divisor found if qm=0 then nf_cache[lm] := FAIL; # cache "no divisor" A[i,basis_lookup[lm]] := lc; p := expand(p - lc*lm) mod char; end if end if; lc, lm := lt(p); end do; NULL end proc; # compute an FGLM system AX=B # input: Groebner basis, GB term order, desired order ModuleApply := proc(G1::list, tord::function, tord2::function) local lt, G, vars, nextmon, mbasis, nf_proc, n, m, b, A, B, MA, MB, Alookup, ai, bi, v, i, j, stats, TIMER; # stores the time for various parts of the algorithm # [1] total time [2] normal form [3] linear algebra stats := rtable([time(), 0, 0]); # userinfo(2, {'_fglm', 'fglm'}, `NoName`, `-> FGLM Linear System`); char := `if`(type(args[-1],'posint'), args[-1], 0); if G1=[1] or G1=[] then # degenerate case return FAIL; elif char > 2^25 then error "characteristic must be less than 2^25" elif char = 0 then nf_proc := rat_nf elif not isprime(char) then error "characteristic must be zero or prime" else nf_proc := mod_nf end if; # sort the variables in the new order lt := evaluate_termorder(tord2): vars := sort([op(indets(tord2,'name'))], (a,b)->lt(a+b)[2]=a); lt := evaluate_termorder(tord): G := [seq([lt(i),i], i=G1)]; nextmon := evaluate_nextmon(tord2); mbasis := table(); n := 0; # n = nops(mbasis) m := 1; # current monomial b := {seq(i[2], i=G)}; # border monomials of G (for now) # compute a monomial basis for the quotient (assumed finite dimensional) while m <> FAIL do n := n+1; mbasis[n] := m; m := nextmon(m,b) end do; userinfo(3, {'_triset', '_fglm', 'fglm'}, `NoName`, cat(` the quotient has `,n,` monomials`)); # construct a lookup table for monomials (used by nf_proc) mbasis := [seq(mbasis[i],i=1..n)]; basis_lookup := table(zip(`=`,mbasis,[seq(i,i=1..n)])); # initialize XLA to build an n by n matrix A # the independent remainders are the rows A := rtable(1..n,1..n,'subtype'='Matrix','order'='C_order'); XLA:-Initialize(A, 'characteristic'=char); ai := 1; if char = 0 then userinfo(3, {'_triset', '_fglm', 'fglm'}, `NoName`, cat(` dependency check with `, XLA:-N, ` primes`)); userinfo(3, {'_triset', '_fglm', 'fglm'}, `NoName`, ` probability of error`, evalf((1-(1-1/2^25)^n)^XLA:-N, 2)) end if; # stores the dependent remainders as rows in a table B := table(); bi := 0; # nops B # independent monomials (whose remainders are in A) MA := rtable(1..n); # dependent monomials (whose remainders are in B) MB := table(); # lookup table for A: monomial -> row Alookup := table([1=1]): # cache used in normal form computations nf_cache := table([seq(i=FAIL, i=mbasis)], 'sparse'): # run FGLM mod p to determine which remainders are independent # note: # - for monomials m1, m2, if m1 < m2 then the remainder of m1 has already been computed # - we choose the largest variable (wrt tord2) which divides m and compute # normalform(i*normalform(m/i)) = normalform(m) m := 1; # current monomial b := {}; # border monomials while ai <= n and m <> FAIL do userinfo(4, {'_triset', '_fglm', 'fglm'}, `NoName`, ` current monomial `, m); TIMER := time(): v := select(member, vars, indets(m,'name')); # find variable dividing m v := remove(a->evalb(Alookup[m/a]=0), v); if nops(v)=0 then nf_proc(m, G, lt, A, ai) # if no such variable, compute remainder of m else i := Alookup[m/v[1]]; # otherwise compute remainder of v*normalform(m/v) nf_proc(v[1]*add(A[i,j]*mbasis[j], j=1..n), G, lt, A, ai) end if; stats[2] := stats[2] + time() - TIMER; TIMER := time(): if XLA:-IsIndependent(ai) then MA[ai] := m; Alookup[m] := ai; ai := ai+1; else bi := bi+1; B[bi] := rtable(A[ai,1..-1]); for i from 1 to n do # zero out A A[ai,i] := 0 end do; MB[bi] := m; b := b union {m}; end if; stats[3] := stats[3] + time() - TIMER; m := nextmon(m, b) end do; while m <> FAIL do # the remaining rows are all dependent rows userinfo(4, {'_triset', '_fglm', 'fglm'}, `NoName`, ` current monomial `, m); bi := bi+1; B[bi] := rtable(1..n); v := select(member, vars, indets(m,'name')); v := remove(a->evalb(Alookup[m/a]=0), v); TIMER := time(): if nops(v)=0 then nf_proc(m, G, lt, B[bi], NULL) else i := Alookup[m/v[1]]; nf_proc(add(v[1]*A[i,j]*mbasis[j], j=1..n), G, lt, B[bi], NULL); end if; stats[2] := stats[2] + time() - TIMER; MB[bi] := m; b := b union {m}; m := nextmon(m, b) end do; stats[1] := time()-stats[1]; userinfo(3, {'_fglm', 'fglm'}, `NoName`, ` normalforms (sec)`, evalf(stats[2],3)); userinfo(3, {'_fglm', 'fglm'}, `NoName`, ` linear algebra (sec)`, evalf(stats[3],3)); userinfo(2, {'_fglm', 'fglm'}, `NoName`, ` total time (sec)`, evalf(stats[1],3)); userinfo(3, {'_fglm', 'fglm'}, `NoName`, `----------------------------------`); # clear memory used by XLA XLA:-Initialize(1,1); LinearAlgebra:-LA_Main:-Transpose(A, 'inplace'=true, outputoptions = []), rtable(1..n, 1..bi, [seq(convert(B[i],'list'),i=1..bi)], 'subtype'='Matrix', 'transpose'=true), rtable(1..n, [seq(MA[i], i=1..n)], 'subtype'='Vector[row]'), rtable(1..bi, [seq(MB[i], i=1..bi)], 'subtype'='Vector[row]'); end proc; end module: # triangular set algorithm triset_algorithm := proc(G, tord, vars::list(name), char::nonnegint:=0) local A, B, Ma, Mb, X, INIT, MULT, T, g, dg, keep, i, j, n, N, r, stats, TIMER; # [1] total time [2] normalforms [3] linear algebra [4] chinese rem. # [5] polynomial ops [6] rational reconstruction stats := rtable([time(), 0, 0, 0, 0, 0]); userinfo(2, {'_triset'}, `NoName`, `-> Triangular Set`); INIT := 13; # number of initial images MULT := 1.7: # multiplier for when more images are needed if char > 2^25 then error "the characteristic must be less than 2^25" elif char > 0 and not isprime(char) then error "the characteristic must be zero or prime" end if; TIMER := time(): A, B, Ma, Mb := fglm_system(G, tord, 'plex'(op(vars)), char); stats[2] := stats[2] + time() - TIMER; # initial solve TIMER := time(): XLA:-LUSolve(A, B, 'characteristic'=char, 'numprimes'=INIT); stats[3] := stats[3] + time() - TIMER; TIMER := time(): XLA:-ChineseRemainder(); stats[4] := stats[4] + time() - TIMER; if char > 0 then # output triangular set TIMER := time(): T := convert(Mb - LinearAlgebra:-VectorMatrixMultiply(Ma, XLA:-C), 'list'); g := T[1]; if g=1 then T := [1]; else dg := diff(g, vars[-1]) mod char; Divide(dg, Gcd(dg, g) mod char, 'dg') mod char; T := [g, seq(expand(Rem(dg*i, g, vars[-1]) mod char) mod char, i=T[2..-1])]; T := [seq(i/lcoeff(i,vars) mod char, i=T)]; end if; stats[5] := stats[5] + time() - TIMER; else n := op(1,Ma); # number of rows N := op(1,Mb); # number of columns i := 1; # current column T := rtable(1..N); userinfo(4, {'_triset', '_fglm'}, `NoName`, cat(` starting batch of `, XLA:-N, ` primes`)); while i = 1 do # reconstruct first polynomial TIMER := time(): r := xratrecon(Mb[i] - add(XLA:-C[j, 1]*Ma[j], j=1..n), vars, XLA:-Cmod); stats[6] := stats[6] + time() - TIMER; if r = FAIL then TIMER := time(): XLA:-NumberOfImages(ceil(XLA:-N * MULT)); stats[3] := stats[3] + time() - TIMER; userinfo(4, {'_triset', '_fglm'}, `NoName`, cat(` increase batch to `, XLA:-N, ` primes`)); TIMER := time(): XLA:-ChineseRemainder(); stats[4] := stats[4] + time() - TIMER; else TIMER := time(): T[1] := iprimpart(r); g := T[1]; dg := diff(g, vars[-1]); divide(dg, gcd(dg,g), 'dg'); stats[5] := stats[5] + time() - TIMER; userinfo(4, {'_triset', '_fglm'}, `NoName`, ` reconstruct poly `,1); i := 2; end if; end do; while i <= N do # reconstruct the rest TIMER := time(): r := Mb[i] - add(XLA:-C[j, i]*Ma[j], j=1..n); r := expand(Rem(dg*r, g, vars[-1]) mod XLA:-Cmod); r := r/lcoeff(r, vars) mod XLA:-Cmod; stats[5] := stats[5] + time() - TIMER; TIMER := time(): r := xratrecon(r, vars, XLA:-Cmod); stats[6] := stats[6] + time() - TIMER; if r = FAIL then TIMER := time(): XLA:-NumberOfImages(ceil(XLA:-N * MULT)); userinfo(4, {'_triset', '_fglm'}, `NoName`, cat(` increase batch to `, XLA:-N, ` primes`)); stats[3] := stats[3] + time() - TIMER; TIMER := time(): XLA:-ChineseRemainder(); stats[4] := stats[4] + time() - TIMER; else TIMER := time(): T[i] := iprimpart(r); stats[5] := stats[5] + time() - TIMER; userinfo(4, {'_triset', '_fglm'}, `NoName`, ` reconstruct poly `,i); i := i+1; end if; end do; T := convert(T, 'list'); end if; # clear memory used by XLA XLA:-Initialize(1,1); stats[1] := time() - stats[1]; userinfo(3, {'_triset', '_fglm'}, `NoName`, ` normalforms (sec)`, evalf(stats[2],3)); userinfo(3, {'_triset', '_fglm'}, `NoName`, ` linear algebra (sec)`, evalf(stats[3],3)); userinfo(3, {'_triset', '_fglm'}, `NoName`, ` chinese remainder (sec)`, evalf(stats[4],3)); userinfo(3, {'_triset', '_fglm'}, `NoName`, ` polynomial ops (sec)`, evalf(stats[5],3)); userinfo(3, {'_triset', '_fglm'}, `NoName`, ` rational reconst. (sec)`, evalf(stats[6],3)); userinfo(2, {'_triset' }, `NoName`, ` total time (sec)`, evalf(stats[1],3)); userinfo(3, {'_triset' }, `NoName`, `----------------------------------`); T, dg; end proc; TriangularSet := proc(J::{list, set, PolynomialIdeal}, X::{list(name), name, function}, Y) local G, tord, vars, char, T, g, rofsub, rofsF, rofsR, rofsmpolys, rofsvars, rtord, rvars, lt; global infolevel; if nargs = 0 then error "insufficient arguments" elif type(J, 'PolynomialIdeal') then char := IdealInfo:-Characteristic(J) else char := `if`(type(args[-1],'posint'), args[-1], 0) end if; if nargs=1 or type(X, 'name') then G := GroebnerBasis(J, 'tord', char); vars := [SuggestVariableOrder(J, indets(tord,'name'))]; if nargs > 1 then X := vars end if; elif type(X, 'function') then if type(J, 'set') then error "a Groebner basis must be given as an ordered list" end if; tord := X; G := `if`(type(J, 'PolynomialIdeal'), GroebnerBasis(J, tord, char), J); if nargs > 2 and type(Y, 'list'(name)) then vars := Y else vars := [SuggestVariableOrder(J, indets(tord,'name'))]; if nargs > 2 and type(Y, 'name') then Y := vars end if; end if; else error "invalid arguments" end if; if not (indets(vars,'name') subset indets(tord,'name')) then error "some variables do not appear in the monomial order" end if; # set infolevels infolevel['_fglm'] := 1; if assigned(infolevel[':-TriangularSet']) then infolevel['_triset'] := infolevel[':-TriangularSet']; else infolevel['_triset'] := 1; end if; # substitute for RootOfs and radicals if nops(indets(G, 'Or'('RootOf','radical','nonreal'))) > 0 then rofsub := true; G := evala(Normal(G), {'expanded'}); rofsF, rofsR, rofsmpolys, rofsvars := substitute_everything(G); rofsvars := [SuggestVariableOrder(rofsmpolys, rofsvars)]; rvars := [op(rofsvars), op(vars)]; rtord := 'prod'(tord, 'tdeg'(op(rofsvars))); G := InterReduceGB(simplify([op(rofsmpolys), op(subs(rofsF, G))], 'symbolic'), rtord, char); else rofsub := false; rofsF := {}; rofsR := {}; rofsvars := []; rvars := vars; rtord := tord; end if; if not ideal_algorithms:-isfinite(G, rtord) then error "the ideal is not finite dimensional" elif not (indets(G,'name') subset indets(rvars,'name')) then error "parameters are not supported" elif char > 2^25 then error "the characteristic must be less than 2^25" elif char <> 0 and not isprime(char) then error "the characteristic must be prime" end if; T, g := triset_algorithm(G, rtord, rvars, char); if rofsub then lt := evaluate_termorder(rtord); rofsvars := {op(rofsvars)}; rofsmpolys, T := selectremove(proc(a) evalb(indets(lt(a)[2]) subset rofsvars) end proc, T); T := subs(rofsR, T); g := subs(rofsR, g); if char=0 then T := [seq(evala(Normal(i),{'expanded'}), i=T)]; g := evala(Normal(g), {'expanded'}); else T := [seq(evala(Normal(i),{'expanded'}), i=T)] mod char; g := evala(Normal(g), {'expanded'}) mod char; end if; end if; if nops(T)=0 then T, g := [1], 1; end if; T, g; end proc; # incremental linear solve - used by FGLM linear_solve := module() local L, LM, n, coftype, vars, cofvars, char; export initialize, reduce, add_equation, reduce_equation; # initialize solver # returns the coefficient type initialize := proc(vars1::Or(set(name),list(name)), cofvars1::Or(set(name),list(name)), char1::nonnegint) # all termorders are equivalent for linear equations # the order of the variables can matter vars := [op(vars1)]; cofvars := {op(cofvars1)}; char := char1; L := table(); n := 0; LM := []; if char=0 and nops(cofvars)=0 then coftype := 'rat_int' elif char=0 then coftype := 'rat_fcn' elif nops(cofvars)=0 then coftype := 'mod_int' else coftype := 'mod_fcn' end if end proc; # adds an equation to the system # gaussian elimination is performed on the left hand side # while the right hand side can be used to track dependencies # in the FGLM algorithm, we use remainder = monomial # so that when a dependency is found we get the polynomial automatically # input can be an equation or an expression, in which case rhs=0 is assumed # returns true if the equation is linearly dependent # the rhs for the dependency is optionally assigned to dep # there is no point in using a modular heuristic here # because you have to solve the system over the rationals anyway add_equation := proc(eqn, dep::name) local f, fm, r; f, fm := `if`(type(eqn, 'equation'), op(eqn), op([eqn,0])); r := reduce[coftype](f, fm); if nargs > 1 then dep := r[4]; end if; if r[1] = 0 then true else n := n+1; L[n] := r; LM := [op(LM), r[2]]; false end if end proc; # reduces an equation but does not add it to the system # returns whatever format was given # not used, but good to have reduce_equation := proc(eqn) local r; if type(eqn, 'equation') then r := reduce[coftype](op(eqn)); r[3]=r[4] else reduce[coftype](eqn, 0)[3] end if end proc; # its a linear system, you can do this in one pass # SPEED CRITICAL reduce := table([ # linear reduction - rational coefficients # basis format : [lc, lm, poly, monomials] # input : polynom, monomial 'rat_int' = proc(f, fm) local p, r, rm, qc, lc, lm, i; p := f; r := 0; rm := fm; while p <> 0 do lc := lcoeff(p, vars, 'lm'); if member(lm, LM, 'i') then qc := lc/L[i][1]; p := expand(denom(qc)*p - numer(qc)*L[i][3]); rm := expand(denom(qc)*rm - numer(qc)*L[i][4]); r := denom(qc)*r else r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; qc := icontent(r); if qc <> 0 then [lcoeff(r/qc, vars, 'lm'), lm, r/qc, rm/qc] else [lcoeff(r, vars, 'lm'), lm, r, rm] end if end proc, # linear reduction - rational function coefficients # basis format : [lc, lm, poly, monomials] # input : polynom, monomial 'rat_fcn' = proc(f, fm) local p, r, rm, qc, lc, lm, i; p := f; r := 0; rm := fm; while p <> 0 do lc := lcoeff(p, vars, 'lm'); if member(lm, LM, 'i') then qc := normal(lc/L[i][1]); p := normal(denom(qc)*p - numer(qc)*L[i][3]); rm := denom(qc)*rm - numer(qc)*L[i][4]; r := denom(qc)*r else r := r + lc*lm; p := normal(p - lc*lm) end if; end do; qc := content(r, vars, 'r'); if qc <> 0 then [lcoeff(r, vars, 'lm'), lm, r, normal(rm/qc)] else [lcoeff(r, vars, 'lm'), lm, r, normal(rm)] end if end proc, # linear reduction - integer coefficients mod p # basis format : [lc, lm, poly, monomials] # input : polynom, monomial 'mod_int' = proc(f, fm) local p, r, rm, qc, lc, lm, i; p := f; r := 0; rm := fm; while p <> 0 do lc := lcoeff(p, vars, 'lm'); if member(lm, LM, 'i') then qc := lc/L[i][1] mod char; p := expand(p - qc*L[i][3]) mod char; rm := expand(rm - qc*L[i][4]) mod char else r := expand(r + lc*lm); p := expand(p - lc*lm) end if; end do; [lcoeff(r, vars, 'lm'), lm, r, rm] end proc, # linear reduction - rational function coefficients mod p # basis format : [lc, lm, poly, monomials] # input : polynom, monomial 'mod_fcn' = proc(f, fm) local p, r, rm, qc, lc, lm, i; p := f; r := 0; rm := fm; while p <> 0 do lc := lcoeff(p, vars, 'lm'); if member(lm, LM, 'i') then qc := Normal(lc/L[i][1]) mod char; p := Normal(denom(qc)*p - numer(qc)*L[i][3]) mod char; rm := denom(qc)*rm - numer(qc)*L[i][4]; r := denom(qc)*r else r := r + lc*lm; p := Normal(p - lc*lm) mod char end if; end do; r := Primpart(r, vars, 'qc') mod char; if r <> 0 then [lcoeff(r, vars, 'lm'), lm, r, Normal(rm/qc) mod char] else [lcoeff(r, vars, 'lm'), lm, r, Normal(rm) mod char] end if end proc ]); end module; # FGLM algorithm - uses the 'noscale' normalform procedures # and the linear solver above # gb : starting Groebner basis # tord : term order for gb # tord2 : desired term order # char : characteristic # optional arguments # # degree_bound=posint : don't consider monomials above this degree # stop_proc=procedure : halt algorithm if current monomial satisfies this condition # please note that this algorithm performs poorly on dense systems of high degree, # especially over the rational numbers. You may want to try the Groebner Walk. # if N is the number of monomials in a vector space basis of the quotient # then FGLM effectively performs Gaussian elimination on an N x N matrix # with entries comparable to the coefficients in the desired Groebner basis fglm_algorithm := proc(gb1, tord, tord2, char) local lt, nextmon, vars, cofvars, degreebound, stopproc, opts, nf, gb, G, N, R, M, m, r, mprev, rprev, q, i, stats, st, A, B, Ma, Mb; global NoName; userinfo(2, {'_fglm', 'fglm'}, 'NoName', `-> FGLM algorithm`); # stats : [# of monomials, reduction time, linear solve time, total time] stats := array([0,0,0,time()]); opts := select(type, [args], 'equation'); vars := indets(tord, 'name'); cofvars := indets(gb1, 'name') minus vars; if char=0 then userinfo(3, {'_fglm', 'fglm'}, 'NoName', ` rational coefficients`) elif nops(cofvars) > 0 or opts <> [] or char > 2^25 then # run the old code userinfo(3, {'_fglm', 'fglm'}, 'NoName', ` modular coefficients`) else userinfo(3, {'_fglm', 'fglm'}, 'NoName', ` machine size modular coefficients`); A, B, Ma, Mb := fglm_system(gb1, tord, tord2, char); XLA:-LUSolve(A, B, 'characteristic'=char); XLA:-ChineseRemainder(); G := convert(Mb - LinearAlgebra:-VectorMatrixMultiply(Ma, XLA:-C), 'list') mod char; return G; end if; # nice naming trick, initialize returns the coefficient type nf := cat(linear_solve:-initialize(vars, cofvars, char), "_noscale"); lt := evaluate_termorder(tord); nextmon := evaluate_nextmon(tord2); opts := table(['degree_bound'=infinity, 'stop_proc'=proc(m) option inline; false end proc, op(opts)]); degreebound := opts['degree_bound']; stopproc := eval(opts['stop_proc']); gb := [seq([lt(i),i,0], i=gb1)]; G := table(); N := 0; M := {}; m := 1; # idea : rather than reduce a monomial directly # use a previous remainder (when possible) to speed up computation # for lexicographic order we can just try to use the previous remainder # for general orders, we may have to search # separate optimization for plex if op(0, tord2)='plex' then mprev := 1; rprev := [1, 1]; while m <> FAIL and not stopproc(m) do userinfo(4, {'_fglm', 'fglm'}, 'NoName', ` current monomial`,m); stats[1] := stats[1] + 1; st := time(): # check if we can use the previous remainder if divide(m, mprev, 'q') then r := normalform_procs[nf](q*rprev[1]/rprev[2], gb, lt, char, vars) else r := normalform_procs[nf](m, gb, lt, char, vars) end if; st := time()-st; stats[2] := stats[2]+st; userinfo(5, {'_fglm', 'fglm'}, 'NoName', ` reduce sec`,st); mprev := m; rprev := r; st := time(): # true if r is linearly dependent if linear_solve:-add_equation(r[1]=r[2]*m, 'q') then N := N + 1; G[N] := iprimpart(q); M := M union {m} elif degree(m) >= degreebound then M := M union {m} end if; st := time()-st; stats[3] := stats[3]+st; userinfo(5, {'_fglm', 'fglm'}, 'NoName', ` linear sec`,st); m := nextmon(m, M); end do; # general term order # you might think this is an awful lot of duplicity for just a little speedup, # but the general case requires double the memory of lexicographic order # on a large system, that can easily be all the memory in your machine else # R is a table of the remainders, indexed by monomial R := table(); while m <> FAIL and not stopproc(m) do userinfo(4, {'_fglm', 'fglm'}, 'NoName', ` current monomial`,m); stats[1] := stats[1] + 1; st := time(): # search for a large divisor of the current monomial r := 'r'; for i in indets(m) do if assigned(R[m/i]) then r := normalform_procs[nf](R[m/i]*i, gb, lt, char, vars); break end if; end do; if not assigned(r) then r := normalform_procs[nf](m, gb, lt, char, vars) end if; st := time()-st; stats[2] := stats[2]+st; userinfo(5, {'_fglm', 'fglm'}, 'NoName', ` reduce sec`,st); R[m] := r[1]/r[2]; st := time(): # true if r is linearly dependent if linear_solve:-add_equation(r[1]=r[2]*m, 'q') then N := N + 1; G[N] := iprimpart(q); M := M union {m} elif degree(m) >= degreebound then M := M union {m} end if; st := time()-st; stats[3] := stats[3]+st; userinfo(5, {'_fglm', 'fglm'}, 'NoName', ` linear sec`,st); m := nextmon(m, M); end do; end if; # IMPORTANT - free some memory linear_solve:-initialize(vars, cofvars, char); stats[4] := time()-stats[4]; userinfo(2, {'_fglm', 'fglm'}, 'NoName', ` total monomials :`, stats[1]); userinfo(3, {'_fglm', 'fglm'}, 'NoName', ` normalforms (sec)`, evalf(stats[2],3)); userinfo(3, {'_fglm', 'fglm'}, 'NoName', ` linearsolve (sec)`, evalf(stats[3],3)); userinfo(2, {'_fglm', 'fglm'}, 'NoName', ` total time (sec)`, evalf(stats[4],3)); userinfo(3, {'_fglm', 'fglm'}, 'NoName', `-----------------------------`); if char=0 and nops(cofvars) > 0 then expand([seq(primpart(G[i],vars), i=1..N)]) elif nops(cofvars) > 0 then expand([seq(Primpart(numer(G[i]),vars) mod char, i=1..N)]) mod char else [seq(G[i], i=1..N)] end if end proc; ## User command for FGLM, handles RootOfs FGLM := proc(igb::list, itord1, itord2) local gb, tord1, tord2, vars, char, rofsF, rofsR, rofsmpolys, rofsvars, rofsub, lt, i; if type(args[-1],'posint') then char := args[-1]; gb := expand(igb) mod char else char := 0; gb := expand(igb) end if; vars := indets(itord1, 'name'); if nops(indets(igb, 'Or'('RootOf','radical','nonreal'))) > 0 then rofsub := true; gb := evala(Normal(gb), {'expanded'}); rofsF, rofsR, rofsmpolys, rofsvars := substitute_everything(gb); rofsvars := [SuggestVariableOrder(rofsmpolys, rofsvars)]; tord1 := 'prod'(itord1, 'tdeg'(op(rofsvars))); tord2 := 'prod'(itord2, 'tdeg'(op(rofsvars))); # simplify/symbolic is safe here because RootOfs with different indices are given distinct variables gb := InterReduceGB(simplify([op(rofsmpolys), op(subs(rofsF, gb))], 'symbolic'), tord1, char); else rofsub := false; rofsF, rofsR, rofsvars := {}, {}, []; tord1 := itord1; tord2 := itord2; end if; gb := fglm_algorithm(gb, tord1, tord2, char, args[4..-1]); if rofsub then lt := evaluate_termorder(tord2); rofsvars := {op(rofsvars)}; gb := remove(proc(a) evalb(indets(lt(a)[2],'name') subset rofsvars) end proc, gb); gb := subs(rofsR, gb); if char=0 then gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))]; else gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))] mod char; end if; end if; gb; end proc; ## ## Groebner Walk algorithm gbwalk_algorithm := module() local multideg_product, linsolve; export initial_form, tord2mat, gbwalk, PERTURB_INITIAL, PERTURB_TARGET, INTER_REDUCE, INTER_REDUCE_COUNT; PERTURB_INITIAL := false; # perturb the initial vector PERTURB_TARGET := true; # perturb the target vector INTER_REDUCE := true; # inter-reduce intermediate Groebner bases INTER_REDUCE_COUNT := 2: # number of steps to do before inter_reducing # convert a term order, ie: tdeg(x,y,z) into a matrix representation # optionally reorder the columns of the matrix to match VARS # you are not expected to follow the lexdeg branch :) # SPEED IS IMPORTANT tord2mat := proc(tord, VARS) local vars, M, M2, T, i, j, k, m, n, zap_nested_lists, flatten_termorder; if op(0,tord)='plex' then vars := [op(tord)]; n := nops(vars); M := [seq([seq(0,i=1..j), 1, seq(0,i=j+2..n)], j=0..n-1)] elif op(0,tord)='grlex' then vars := [op(tord)]; n := nops(vars); if n=1 then M := [[1]] else M := [[seq(1,j=1..n)], seq([seq(0,i=1..j), 1, seq(0,i=j+2..n)], j=0..n-2)] end if elif op(0,tord)='tdeg' then vars := [op(tord)]; n := nops(vars); if n=1 then M := [[1]] else M := [[seq(1,j=1..n)], seq([seq(0,i=j+2..n), -1, seq(0,i=1..j)], j=0..n-2)] end if elif op(0,tord)='lexdeg' then vars := [op(tord)]; n := nops(vars); # nobody is expected to follow this M := [seq(seq( [seq(seq(0,m=1..nops(vars[k])), k=1..i-1), op(j), seq(seq(0,m=1..nops(vars[k])), k=i+1..n)], j=procname('tdeg'(op(vars[i])))), i=1..n)] elif op(0,tord)='prod' then flatten_termorder := proc(a) local i; `if`(op(0,a)='prod', seq(procname(i),i=[op(a)]), a) end proc: M2 := table(): vars := table(): T := [flatten_termorder(tord)]; for i from 1 to nops(T) do M2[i] := procname(T[i],vars[i]); end do; vars := [seq(vars[i],i=1..nops(T))]; n := nops(vars); M := [seq(seq( [seq(seq(0,m=1..nops(vars[k])), k=1..i-1), op(j), seq(seq(0,m=1..nops(vars[k])), k=i+1..n)], j=M2[i]), i=1..n)] elif op(0,tord)='wdeg' then T, vars := op(tord); n := nops(vars); if n=1 then M := [[1]] else M := [[seq(T[j],j=1..n)], seq([seq(0,i=j+2..n), -1, seq(0,i=1..j)], j=0..n-2)] end if elif op(0,tord)='mat' or op(0,tord)='matrix' then M := [op(op(1,tord))]; vars :=op(2,tord) end if; # swap columns of the matrix so the variable order # corresponds with VARS zap_nested_lists := proc(a) option inline; `if`(type(a,'list'), op(map(procname, a)), a) end proc: vars := map(zap_nested_lists, vars); if nargs=1 then return M elif type(VARS, 'name') then VARS := vars; return M elif not {op(VARS)} subset {op(vars)} then error "some variables do not appear in the monomial order" end if; M2 := table([]); for i from 1 to nops(vars) do M2[vars[-i]] := [seq(M[j][-i], j=1..nops(M))]; end do; [seq([seq(M2[VARS[i]][j], i=1..nops(VARS))], j=1..nops(M))] end proc: # computes the initial form of a polynomial wrt a vector # SPEED IMPORTANT initial_form := proc(f, W::list(rational), vars::list(name)) local weights, w, i; weights := lcm(seq(denom(i),i=W))*W; if nops(weights) <> nops(vars) then error "there must be exactly one weight per variable" end if; lcoeff(subs({seq(vars[i]=vars[i]*w^weights[i], i=1..nops(vars))}, f), w) end proc: # computes the dot product of a monomial with a vector # outputs a set if input is a polynomial # SPEED CRITICAL multideg_product := proc(f, vect, vars) option inline; `if`(op(0,f)=`+`, op(map(procname, [op(f)], vect, vars)), inner(vect, map2(degree, f, vars))) end proc: # solves a linear expression in t # SPEED CRITICAL linsolve := proc(eqn, t) option inline; `if`(type(eqn, 'rational'), NULL, -tcoeff(eqn,t)/lcoeff(eqn,t)) end proc: gbwalk := proc(initgb::list, itord::function, ftord::function, char::nonnegint) local V, V1, IM, FM, target, gb, iform, vars, tord, Ci, t, d, d2, min_t, S, T, i, j, k, n, TIMER, STEP_TIMER, stats, loop_count; global NoName; userinfo(2, {'_walk', 'walk'}, 'NoName', `-> Groebner Walk`); if char=0 then userinfo(3, {'_walk', 'walk'}, 'NoName', ` rational coefficients`) else userinfo(3, {'_walk', 'walk'}, 'NoName', ` modular coefficients`) end if; # stats [ # of steps, gb time, interreduce time, compute vector time, total time] stats := array(1..5, [0,0,0,0,time()]); FM := tord2mat(ftord, 'vars'); # vars is assigned IM := tord2mat(itord, vars); n := nops(vars); d2 := max(seq(degree(i,{op(vars)}), i=initgb)); d := (d2^2+2*d2)^(2^(n-1)); # initial vector must refine initial term order V := expand(add(d^(n-i)*IM[i], i=1..n)); if PERTURB_INITIAL and nops(vars) > 1 then userinfo(4, {'_walk', 'walk'}, 'NoName', ` perturbing initial vector`); V := V + 1/(d2+1)*IM[2] else V := V end if; tord := itord; gb := initgb; if PERTURB_TARGET then userinfo(4, {'_walk', 'walk'}, 'NoName', ` perturbing target vector`); target := expand(add(d^(n-i)*FM[i], i=1..n)) else target := FM[1] end if; userinfo(5, {'_walk', 'walk'}, 'NoName', ` walk to`, evalf(V,2)); STEP_TIMER := time(): TIMER:=time(): iform := map(initial_form, gb, V, vars); V1 := expand((1-t)*V + t*target); # SPEED CRITICAL min_t := 1; for i from 1 to nops(gb) do S := [multideg_product(iform[i], V1, vars)]; T := [multideg_product(gb[i]-iform[i], V1, vars)]; min_t := min(min_t, op(select(type, [seq(seq(linsolve(j-k,t), k=S), j=T)], 'positive'))); end do; V := eval(V1,t=min_t); TIMER := time()-TIMER: stats[4] := stats[4] + TIMER; userinfo(5, {'_walk', 'walk'}, 'NoName', ` find next vector`, TIMER); userinfo(5, {'_walk', 'walk'}, 'NoName', ` walking to`, evalf(V,3)); for loop_count to infinity do iform := map(initial_form, gb, V, vars); tord := 'mat'([V,op(FM)], vars); TIMER:=time(): Ci := Buchberger(iform, tord, char, 'track'={seq(i, i=1..nops(iform))}, 'strategy'='degree')[2]; TIMER := time()-TIMER: userinfo(5, {'_walk', 'walk'}, 'NoName', ` lifting step `, TIMER); stats[2] := stats[2] + TIMER; if INTER_REDUCE and irem(loop_count, INTER_REDUCE_COUNT)=0 then TIMER:=time(): gb := InterReduceGB([seq(inner(Ci[i],gb), i=1..nops(Ci))], tord, char); TIMER := time()-TIMER; userinfo(5, {'_walk', 'walk'}, 'NoName', ` interreduce `, TIMER); stats[3] := stats[3] + TIMER elif char=0 then gb := [seq(expand(add(Ci[i][j]*gb[j], j=1..nops(gb))), i=1..nops(Ci))] else gb := [seq(expand(add(Ci[i][j]*gb[j], j=1..nops(gb))), i=1..nops(Ci))] mod char end if; STEP_TIMER := time()-STEP_TIMER: userinfo(4, {'_walk', 'walk'}, 'NoName', ` step time (sec) `, STEP_TIMER); if V = target then break end if; STEP_TIMER := time(): TIMER:=time(): iform := map(initial_form, gb, V, vars); V1 := expand((1-t)*V + t*target); # SPEED CRITICAL min_t := 1; for i from 1 to nops(gb) do S := [multideg_product(iform[i], V1, vars)]; T := [multideg_product(gb[i]-iform[i], V1, vars)]; min_t := min(min_t, op(select(type, [seq(seq(linsolve(j-k,t), k=S), j=T)], 'positive'))); end do; V := eval(V1,t=min_t); TIMER := time()-TIMER: stats[4] := stats[4] + TIMER; userinfo(5, {'_walk', 'walk'}, 'NoName', ` find next vector`, TIMER); userinfo(5, {'_walk', 'walk'}, 'NoName', ` walking to`, evalf(V,3)); end do; if INTER_REDUCE and irem(loop_count,INTER_REDUCE_COUNT)=0 then ; else TIMER := time(): gb := InterReduceGB(gb, tord, char); TIMER := time()-TIMER; stats[3] := stats[3] + TIMER end if; stats[1] := loop_count; stats[5] := time()-stats[5]; userinfo(2, {'_walk', 'walk'}, 'NoName', ` total walk steps :`, stats[1]); userinfo(3, {'_walk', 'walk'}, 'NoName', ` solve vector (sec)`, evalf(stats[4],3)); userinfo(3, {'_walk', 'walk'}, 'NoName', ` lifting step (sec)`, evalf(stats[2],3)); userinfo(3, {'_walk', 'walk'}, 'NoName', ` interreduce (sec)`, evalf(stats[3],3)); userinfo(2, {'_walk', 'walk'}, 'NoName', ` total time (sec)`, evalf(stats[5],3)); userinfo(3, {'_walk', 'walk'}, 'NoName', `-----------------------------`); gb end proc: end module; ## User command for Groebner walk, handles RootOfs GroebnerWalk := proc(igb::list, itord1, itord2) local gb, tord1, tord2, vars, char, rofsF, rofsR, rofsmpolys, rofsvars, rofsub, lt; if type(args[-1],'posint') then char := args[-1]; gb := expand(igb) mod char else char := 0; gb := expand(igb) end if; vars := indets(itord1, 'name'); if nops(indets(igb, 'Or'('RootOf','radical','nonreal'))) > 0 then rofsub := true; gb := evala(Normal(gb), {'expanded'}); rofsF, rofsR, rofsmpolys, rofsvars := substitute_everything(gb); rofsvars := [SuggestVariableOrder(rofsmpolys, rofsvars)]; tord1 := 'prod'(itord1, 'tdeg'(op(rofsvars))); tord2 := 'prod'(itord2, 'tdeg'(op(rofsvars))); # simplify/symbolic is safe here because RootOfs with different indices are given distinct variables gb := InterReduceGB(simplify([op(rofsmpolys), op(subs(rofsF, gb))], 'symbolic'), tord1, char); else rofsub := false; rofsF, rofsR, rofsvars := {}, {}, []; tord1 := itord1; tord2 := itord2; end if; gb := gbwalk_algorithm:-gbwalk(gb, tord1, tord2, char, args[4..-1]); if rofsub then lt := evaluate_termorder(tord2); rofsvars := {op(rofsvars)}; gb := remove(proc(a) evalb(indets(lt(a)[2],'name') subset rofsvars) end proc, gb); gb := subs(rofsR, gb); if char=0 then gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))]; else gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))] mod char; end if; end if; gb; end proc; # returns a sequence of variables in an order that (hopefully) results # in a fast Groebner basis computation - mainly for lexicographic or elimination orders # note : if a set is a Groebner basis with respect to some lexicographic order then # SuggestVariableOrder should return that order SuggestVariableOrder := proc(J::Or(PolynomialIdeal,list,set)) option system, remember; local G, X, V, W, i; if type(J, 'PolynomialIdeal') then G := [op(IdealInfo:-Generators(J))]; X := `if`(nargs >= 2, args[2], IdealInfo:-Variables(J)) elif type(J, 'Or'('list','set')) then G := expand([op(J)]); X := `if`(nargs >= 2, args[2], indets(J,'name') minus indets(indets(J,'Or'('radical','RootOf')),'name')) end if; if nops(G) = 0 or nops(X)=0 then return op(X) end if; X := [op(X)]; W := [seq(map(degree, G, i), i=X)]; V := [seq([X[i],max(op(W[i])),`+`(op(W[i])),length(map(lcoeff,G,X[i]))], i=1..nops(X))]; seq(i[1], i=sort(V, proc(a,b) evalb(a[2] < b[2] or a[2]=b[2] and (a[3] < b[3] or a[3]=b[3] and a[4] < b[4])) end proc)); end proc: # compute a Groebner basis # handles RootOfs and radicals # SPEED IMPORTANT GroebnerBasis := proc(J, tord1) global infolevel; local gens, tord, vars, badv, char, KGB, opts, meth, gb, cb, itord, igb, rofsub, rofsF, rofsR, rofsmpolys, rofsvars, rtord, lt, i, savegb; # get generators, set up term order if type(J, 'PolynomialIdeal') then gens := [op(IdealInfo:-Generators(J))]; char := IdealInfo:-Characteristic(J); if type(tord1, 'name') then tord := IdealInfo:-DefaultMonomialOrder(J, args[3..-1]); tord1 := tord else tord := tord1 end if; vars := indets(tord, 'name'); badv := IdealInfo:-NonVariables(J); if not hasoption(select(type, [op(J)], 'equation'), 'known_groebner_bases'=table, 'KGB') then KGB := table(); end if elif type(J, 'Or'('list','set')) then gens := [op(J)]; char := `if`(type(args[-1], 'posint'), args[-1], 0); badv := indets(indets(gens, 'Or'('RootOf','radical','nonreal')), 'name'); if type(tord1, 'name') then vars := select(type, [args][3..-1], 'set'('name')); vars := `if`(nops(vars) > 0, vars[1], indets(gens, 'name') minus badv); if nargs > 2 and member(args[3], {'plex', 'grlex'}) then tord := (args[3])(SuggestVariableOrder(gens, vars)) else tord := 'tdeg'(SuggestVariableOrder(gens, vars)) end if; tord1 := tord else tord := tord1; vars := indets(tord, 'name') end if; # uses too much memory # KGB := gb_table({op(gens)}, char) KGB := table() else error "expects its 1st argument, J, to be of type Or(PolynomialIdeal, list, set), but received %1", J end if; if nops(vars intersect badv) > 0 then error "term order variables cannot appear inside a RootOf or radical" end if; opts := table(['method'='default', op(select(type, [args][3..-1], 'equation'))]); meth := StringTools:-LowerCase(convert(opts['method'], 'string')); opts['method'] := evaln(opts['method']); savegb := not (assigned(opts['stop_proc']) or assigned(opts['degree_bound'])); opts := op(op(op(opts))); # convert back to a sequence if meth="default" and assigned(KGB[tord]) then return KGB[tord] end if; # set up infolevels # the underscore ones are inherited from infolevel['GroebnerBasis'] # we set infolevel['_buchberger'] := 1 for Groebner Walk if assigned(infolevel[':-GroebnerBasis']) then infolevel['_f4'] := infolevel[':-GroebnerBasis']; infolevel['_buchberger'] := infolevel[':-GroebnerBasis']; infolevel['_fglm'] := infolevel[':-GroebnerBasis']; infolevel['_walk'] := infolevel[':-GroebnerBasis'] else infolevel['_f4'] := 1; infolevel['_buchberger'] := 1; infolevel['_fglm'] := 1; infolevel['_walk'] := 1 end if; # if we are converting, get a Groebner basis # prefer similar term orders for which a Groebner basis is known # for the default conversion method, use this basis to determine whether # the ideal is zero-dimensional and choose FGLM or Groebner Walk accordingly # assign term order to 'itord' and Groebner basis to 'igb' if meth="fglm" or meth="walk" or meth="default" and member(op(0,tord), {'plex', 'lexdeg', 'prod'}) then itord := similar_orders(map(op, [indices(KGB)]), tord); if nops(itord)=0 then itord := 'tdeg'(SuggestVariableOrder(gens, vars)); KGB[itord] := procname(gens, itord, opts, 'method'='direct', char) elif nops(itord)=1 then itord := itord[1] else itord := sort(map(proc(a) [a, length(KGB[a])] end proc, itord), proc(a,b) option inline; evalb(a[2] < b[2]) end proc)[1][1] end if; igb := KGB[itord]; if meth="default" then meth := `if`(ideal_algorithms:-isfinite(igb, itord), "fglm", "walk") end if; elif meth="default" then meth := "direct"; itord := 'tdeg'(); igb := [] end if; # substitute for RootOfs and radicals # we add their minimal polynomials to the generating set # and extend the term orders to handle them # if we are converting, we interreduce the modified initial Groebner basis # to ensure that it is correct if nops(indets(gens, 'Or'('RootOf','radical','nonreal'))) > 0 then rofsub := true; gens := evala(Normal(gens), {'expanded'}); rofsF, rofsR, rofsmpolys, rofsvars := substitute_everything(gens); rofsvars := [SuggestVariableOrder(rofsmpolys, rofsvars)]; gens := [op(rofsmpolys), op(subs(rofsF, gens))]; rtord := 'prod'(tord, 'tdeg'(op(rofsvars))); itord := 'prod'(itord, 'tdeg'(op(rofsvars))); if meth="fglm" or meth="walk" then # simplify/symbolic is safe here because RootOfs with different indices are given distinct variables igb := InterReduceGB(simplify([op(rofsmpolys), op(subs(rofsF, igb))], 'symbolic'), itord, char); end if else rofsub := false; rofsF := {}; rofsR := {}; rofsvars := []; rtord := tord; end if; # call the appropriate method cb := []; if (meth="direct" and indets(gens, 'name') subset indets(rtord,'name') and (char=0 or char < 2^25)) or (meth="F4" or meth="f4") then gb := f4_algorithm:-f4(gens, rtord, char) elif meth="direct" or meth="bb" or meth="buchberger" then gb := Buchberger(gens, rtord, char, opts, 'track'={}) elif meth="fglm" then gb := fglm_algorithm(igb, itord, rtord, char, opts) elif meth="walk" then infolevel['_buchberger'] := 1; gb := gbwalk_algorithm:-gbwalk(igb, itord, rtord, char, opts) elif meth="extended" then gb, cb := Buchberger(gens, rtord, char, opts, 'track'={seq(i, i=1..nops(gens))}) elif meth="qr" then gb, cb := Buchberger(gens, rtord, char, opts, 'track'={1}); cb := map(op, cb) else error "unknown method", meth end if; if rofsub then lt := evaluate_termorder(rtord); rofsvars := {op(rofsvars)}; rofsmpolys, gb := selectremove(proc(a) evalb(indets(lt(a)[2]) subset rofsvars) end proc, gb); gb := subs(rofsR, gb); if char=0 then if nops(cb) = 0 then gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))]; else cb := subs(rofsR, [seq(i[nops(rofsmpolys)+1..-1], i=cb[nops(rofsmpolys)+1..-1])]); gb := map(proc(a) local c; [evala(Normal(evala(Primpart(a,vars,'c'))), {'expanded'}),c] end proc, gb); cb := [seq(`if`(gb[i][2]=1, cb[i], [seq(evala(Normal(j/gb[i][2]), {'expanded'}), j=cb[i])]), i=1..nops(gb))]; gb := [seq(i[1], i=gb)]; end if; else if nops(cb) = 0 then gb := [seq(evala(Normal(i),{'expanded'}), i=map(evala@Primpart, gb, vars))] mod char; else cb := subs(rofsR, [seq(i[nops(rofsmpolys)+1..-1], i=cb[nops(rofsmpolys)+1..-1])]) mod char; gb := map(proc(a) local c; [evala(Normal(evala(Primpart(a,vars,'c'))), {'expanded'}),c] end proc, gb); cb := [seq(`if`(gb[i][2]=1, cb[i], [seq(evala(Normal(j/gb[i][2]), {'expanded'}), j=cb[i])]), i=1..nops(gb))] mod char; gb := [seq(i[1], i=gb)] mod char; end if; end if end if; gb := `if`(nops(gb)=0, [1], gb); if savegb or (meth <> "fglm") then KGB[tord] := gb end if; gb, `if`(nops(cb)=0, NULL, cb) end proc; # compute a Groebner basis for a module wrt TOP or POT order ModuleGB := proc(M, tord, ordertype, char::nonnegint:=0) local F, n, e, i, j, V, G, mtord; n := nops(M[1]); V := [seq('e'[i], i=1..n)]; F := [seq(seq('e'[i]*'e'[j], j=1..i), i=1..n), op(map(inner, M, V))]; if ordertype='POT' or ordertype='pot' then mtord := 'prod'('plex'(op(V)), tord); elif ordertype='TOP' or ordertype='top' then mtord := 'prod'(tord, 'plex'(op(V))); else error "only TOP or POT order is supported" end if; G := GroebnerBasis(F, mtord, op(select(type,[args],'equation')), char); G := remove(a->degree(a,{op(V)}) > 1, G); G := map(a->map2(coeff, a, V), G); end proc; # returns the leading monomial of a polymomial wrt a term order # note : if the input is an ideal, this returns the ideal of leading monomials LeadingMonomial := proc(f, lt, LC) option system, remember; local lc, lm; if type(f, 'Or'('list','set')) then map(procname, f, lt) elif type(f, 'PolynomialIdeal') then PolynomialIdeal(map(procname, GroebnerBasis(f,lt), lt), select(type, [op(f)], 'equation')) else if type([args][-1], 'posint') then lc, lm := evaluate_termorder(lt)(normal(f) mod args[-1]) else lc, lm := evaluate_termorder(lt)(normal(f)) end if; if nargs >= 3 and type(LC, 'name') then LC := lc end if; lm end if end proc; # compute the leading coefficient, optionally assign the leading monomial to LM LeadingCoefficient := proc(f, lt, LM) local lc, lm; lc := NULL; lm := LeadingMonomial(f, lt, 'lc'); if nargs > 2 and type(LM, 'name') then LM := lm end if; lc end proc; # compute the trailing monomial, optionally assign trailing coefficient to a third argument TrailingMonomial := proc(f, tord, TC) local v; 1/LeadingMonomial(subs({seq(v=1/v, v=indets(tord,'name'))}, f), tord, args[3..-1]) end proc: # compute the initial form of a polynomial wrt a weight vector InitialForm := proc(f, W::list(rational), vars::list(name), $) if type(f, {list,set}) then map(procname, f, W, vars) elif not type(f, 'polynom'('anything', vars)) then error "argument must be a polynomial in the variables", vars else gbwalk_algorithm:-initial_form(f, W, vars) end if; end proc: # compute the weighted degree of a polynomial WeightedDegree := proc(f, W::list(rational), vars::list(name), $) local M, i; if type(f, {list,set}) then map(procname, f, W, vars) elif not type(f, 'polynom'('anything', vars)) then error "argument must be a polynomial in the variables", vars else max(seq(inner(W, map2(degree,i,vars)), i=[[coeffs(expand(f), vars, 'M')],[M]][2])) end if; end proc; # construct a matrix for a term order, optionally order the columns wrt a list of variables TermOrderMatrix := proc(tord::function, X1::list(name), $) local X; if nargs=1 then X := sort([op(indets(tord,'name'))], (a,b)->evalb(LeadingMonomial(a+b,tord)=a)); else X := X1; end if; Matrix(gbwalk_algorithm:-tord2mat(tord, X)); end proc: # computes the normal form of a polynomial, either with respect # to an ideal (using a Groebner basis), or using a list of polynomials NormalForm := proc(f1, J, tord, Q1) option system, remember; local f, G, char, vars, cofvars, lt, r, s, i, j, rofsvars, rofsF, rofsR, rofsmpolys, rtord; # get a Groebner basis if the input was an ideal if type([J,tord], ['PolynomialIdeal','function']) then char := IdealInfo:-Characteristic(J); G := GroebnerBasis(J, tord) elif type([J,tord], ['list','function']) then char := `if`(nargs > 2 and type(args[-1], 'posint'), args[-1], 0); G := `if`(char = 0, normal(J), Normal(J) mod char) elif type(tord, 'name') then error "please specify a monomial order" elif type(J, 'set') then error "please use a list since the result may depend on the order" else error "cannot determine generating set" end if; vars := indets(tord, 'name'); if nops(indets({f1,G}, 'Or'('RootOf','radical','nonreal'))) > 0 then f := evala(Normal(f1), {'expanded','independent'}); G := evala(Normal(G), {'expanded','independent'}); rofsF, rofsR, rofsmpolys, rofsvars := substitute_everything({f,G}); rofsvars := [SuggestVariableOrder(rofsmpolys, rofsvars)]; if nops(select(member, vars, rofsvars)) > 0 then error "monomial order variables cannot appear inside a RootOf or radical" end if; rtord := 'prod'(tord, 'tdeg'(op(rofsvars))); G := subs(rofsF, G); f := subs(rofsF, f) else rofsF := {}; rofsR := {}; rofsvars := []; rofsmpolys := []; rtord := tord; f := f1 end if; vars := indets(rtord, 'name'); cofvars := indets([f,G], 'name') minus vars; lt := evaluate_termorder(rtord); if char=0 then if nargs > 3 and type(Q1, 'name') then f := f, 0, [1, seq(0, i=1..nops(G))]; G := [seq([lt(i),i,0,[0, seq(0, j=1..nops(G))]], i=rofsmpolys), seq([lt(G[i]),G[i],0,[0, seq(0, j=1..i-1), 1, seq(0, j=i+1..nops(G))]], i=1..nops(G))]; r := normalform_procs['rat_trk_cof'](f, G, lt, 0, vars); Q1 := evala(Normal(subs(rofsR, -r[5][2..-1]/r[5][1])), {'expanded','independent'}); r, s := r[3], r[5][1] elif nops(cofvars)=0 then r, s := op(normalform_procs['rat_int_noscale'](f, [seq([lt(i),i], i=[op(rofsmpolys),op(G)])], lt)) else r, s := op(normalform_procs['rat_fcn_noscale'](f, [seq([lt(i),i], i=[op(rofsmpolys),op(G)])], lt, 0, vars)) end if else if nargs > 3 and type(Q1, 'name') then f := f, 0, [1, seq(0, i=1..nops(G))]; G := [seq([lt(i),i,0,[0, seq(0, j=1..nops(G))]], i=rofsmpolys), seq([lt(G[i]),G[i],0,[0, seq(0, j=1..i-1), 1, seq(0, j=i+1..nops(G))]], i=1..nops(G))]; r := normalform_procs['mod_trk_cof'](f, G, lt, char, vars); Q1 := evala(Normal(subs(rofsR, -r[5][2..-1]/r[5][1])), {'expanded','independent'}); r, s := r[3], r[5][1] elif nops(cofvars)=0 then r, s := op(normalform_procs['mod_int_noscale'](f, [seq([lt(i),i], i=[op(rofsmpolys),op(G)])], lt, char)) else r, s := op(normalform_procs['mod_fcn_noscale'](f, [seq([lt(i),i], i=[op(rofsmpolys),op(G)])], lt, char, vars)) end if end if; evala(Normal(subs(rofsR, r/s)), {'expanded'}); end proc; # compute the smallest degree univariate polynomial in an ideal # this will be a generator for k[v] intersect J # uses a lexicographic Groebner Basis and resultants for ideals # uses an elimination Groebner Basis for non-ideals UnivariatePolynomial := proc(v, J, X) option system, remember; local G, char, vars, tord; # handle reversed 1st and 2nd argument if type(J, 'name') and type(v, 'Or'('PolynomialIdeal','list','set')) then return procname(J, v, args[3..-1]) end if; # determine default variables, characteristic if type(J, 'PolynomialIdeal') then vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J) elif type(J, 'Or'('list','set')) then vars := (indets(J,'name') minus indets(indets(J,'Or'('RootOf','radical')), 'name')); char := `if`(type(args[-1], 'posint'), args[-1], 0) else error "expects its 2nd argument, J, to be of type Or(PolynomialIdeal, list, set), but received %1", J end if; # override with 3rd argument variables if nargs >= 3 and type(args[3], 'Or'('list'('name'),'set'('name'))) then vars := {op(args[3])} elif nargs >= 3 and not type(args[3],'nonnegint') then error "expects its 3rd argument, X, to be of type set(name), but received %1", X end if; # see if the result has been computed already # first check for a compatible Groebner basis # otherwise search all Groebner bases hoping to get lucky if type(J, 'PolynomialIdeal') then tord := select(is_elim_order, IdealInfo:-KnownGroebnerBases(J, vars), {v}); if nops(tord) > 0 then tord := tord[1]; G := remove(has, GroebnerBasis(J, tord), vars minus {v}); return `if`(nops(G)=0, 0, G[1]) end if; G := map2(GroebnerBasis, J, IdealInfo:-KnownGroebnerBases(J, vars)); G := map2(remove, has, G, vars minus {v}) minus {[]}; if nops(G) > 0 then return G[1][1] end if end if; tord := 'lexdeg'(remove(member, [SuggestVariableOrder(J, vars)], {v}), [v]); vars := vars minus {v}; G := remove(has, GroebnerBasis(J, tord, 'stop_proc'=proc(m) has(m, vars) end proc, char), vars); # note for the future : you should really try to store this result # it would be useful for prime/primary decomposition # return 0 if no univariate polynomial exists `if`(nops(G)=0, 0, G[1]); end proc; # eliminate variables from an ideal # computes the intersection of J with the subring k[X] EliminationIdeal := proc(J::PolynomialIdeal, X::set(name), $) option system, remember; local X2, U, G, tord, i; X2, U := selectremove(member, IdealInfo:-Variables(J) intersect indets(IdealInfo:-Generators(J), 'name'), X); if nops(U) = 0 then return J elif nops(X2) = 1 then G := UnivariatePolynomial(op(X2), J); return PolynomialIdeal(G, 'characteristic'=IdealInfo:-Characteristic(J), 'variables'=X2, 'known_groebner_bases'=['plex'(op(X2))=[G], 'tdeg'(op(X2))=[G]]) end if; tord := select(is_elim_order, IdealInfo:-KnownGroebnerBases(J, X2 union U), X2); if nops(tord) > 0 then G := [seq(elim_order(i, X2) = remove(has, GroebnerBasis(J, i), U), i=tord)]; return PolynomialIdeal(rhs(G[1]), 'characteristic'=IdealInfo:-Characteristic(J), 'variables'=X2, 'known_groebner_bases'=G) end if; tord := 'lexdeg'(selectremove(member, [SuggestVariableOrder(J, X2 union U)], U)); G := remove(has, GroebnerBasis(J, tord, 'stop_proc'=proc(m) has(m, U) end proc), U); PolynomialIdeal(G, 'characteristic'=IdealInfo:-Characteristic(J), 'variables'=X2, 'known_groebner_bases'=[elim_order(tord, X2)=G]) end proc; # intersect two or more ideals # always uses the original generators # characteristic taken from the first ideal Intersect := proc(J::PolynomialIdeal, K::PolynomialIdeal) option system, remember; local vars, pars, char, gensJ, gensK, gensBoth, Z, Jz, G, T, S, B, n, i, j, unknown; global NoName; unknown := remove(type, [args], 'PolynomialIdeal'); if nops(unknown)>0 then error "invalid argument(s)", op(unknown) end if: char := map(IdealInfo:-Characteristic, {args}); if nops(char) > 1 then error "ideals with different characteristics can not be intersected" end if; char := [op(char),0][1]; vars := `union`(seq(IdealInfo:-Variables(i), i=[args])); pars := `union`(seq(IdealInfo:-NonVariables(i), i=[args])); if nops(vars intersect pars) > 0 then error "the variables of one ideal are parameters of another" end if; # algorithm for two ideals taken from Cox, Little, and O'Shea # improvement : we remove any common generators first if nargs = 2 then gensJ := Generators(J); gensK := Generators(K); if (nops(gensJ)=1 or nops(gensK)=1) then gensBoth := gensJ intersect gensK; if nops(gensBoth)=1 then return PolynomialIdeal(gensBoth, 'variables'=vars, 'characteristic'=char) end if else gensBoth, gensJ := selectremove(member, gensJ, gensK); gensK := remove(member, gensK, gensBoth) end if; Jz := PolynomialIdeal(seq(Z*i, i=gensJ), seq((1-Z)*i, i=gensK), 'characteristic'=IdealInfo:-Characteristic(J), 'variables'=vars union {Z}); PolynomialIdeal(EliminationIdeal(Jz, vars), gensBoth, 'variables'=vars) # for more than two ideals, we select a pair of ideals to intersect at each iteration # we prefer pairs with the most number of polynomials in common, since these are "cheap" # then we sort, preferring pairs which eliminate the most expensive pairs # finally we select the largest pair by total length, hoping to isolate nasty intersections early on # this is the best strategy, don't change it unless you know what you are doing # we must be able to quickly intersect prime and primary decompositions with many (20-50+) components elif nargs > 2 then gensBoth := `intersect`(seq(Generators(i), i=args)); G := table([seq(Generators(i) minus gensBoth, i=[args])]); n := nargs; B := {seq(i, i=1..n)}; T := sort([seq(seq([i,j,nops(G[i])+nops(G[j])-2*nops(G[i] intersect G[j])], j=i+1..n), i=1..n-1)], proc(a,b) option inline; evalb(a[3] < b[3]) end proc); while nops(B) > 1 do userinfo(3, 'Intersect', 'NoName', ` ideals :`, nops(B)); S := select(proc(a) evalb(a[3]=T[1][3]) end proc, T); S := sort(map(proc(a) [a[1],a[2], add(`if`(i[1]=a[1] or i[2]=a[1] or i[2]=a[2] or i[2]=a[2], i[3], 0), i=S), nops(G[a[1]] union G[a[2]])] end proc, S), proc(a,b) option inline; evalb(a[3] > b[3] or (a[3]=b[3] and a[4] > b[4])) end proc); i,j := op(S[1][1..2]); Jz := procname(PolynomialIdeal(G[i], 'variables'=vars, 'characteristic'=char), PolynomialIdeal(G[j], 'variables'=vars, 'characteristic'=char)); B := B minus {i,j}; n := n+1; G[n] := Generators(Jz) minus gensBoth; T := remove(proc(a) evalb(a[1]=i or a[2]=i or a[1]=j or a[2]=j) end proc, T); T := sort([op(T), seq([i, n, nops(G[i])+nops(G[n]) - 2*nops(G[i] intersect G[n])], i=B)], proc(a,b) option inline; evalb(a[3] < b[3]) end proc); B := B union {n}; end do; PolynomialIdeal(gensBoth, G[n], 'characteristic'=char, 'variables'=vars); else args end if end proc; # remove solutions of f from J, ie: compute J:f^infinity # optionally compute an exponent s (not minimal) with J:f^s = J:f^infinity Saturate := proc(J::PolynomialIdeal, f::Or(radalgfun,list(radalgfun),set(radalgfun)), s::name) option system, remember; local Z, vars, char, tord, F, G, C, i, j; vars := IdealInfo:-Variables(J); F := `if`(type(f, 'Or'('list','set')), convert(f,'set'), {f}); if nargs > 2 then char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J); F := `if`(type(f,'Or'('list','set')), convert(f,'list'), [f]); G := GroebnerBasis(J,tord); G, C := GroebnerBasis([1-Z*mul(i,i = F), op(G)], 'prod'('plex'(Z), tord), 'method'='extended', 'zap_pairs'=[{seq(i,i=2..nops(G)+1)}], char); s := max(1, seq(seq(degree(j,Z),j = i[2 .. -1]),i=C)); G := remove(has, G, Z); PolynomialIdeal(G, 'variables'=vars, 'characteristic'=char, 'known_groebner_bases'=[tord = G]) else EliminationIdeal(PolynomialIdeal(J, 1-Z*mul(i,i=F), 'variables'=(vars union {Z})), vars); end if; end proc; Contract := proc(J::PolynomialIdeal, X::set(name), $) option system, remember; local U, gb, tord, f, i; U := X intersect IdealInfo:-Variables(J); gb := GroebnerBasis(J, 'tord', U); f := lcm(seq(LeadingCoefficient(i, tord), i=gb)); Saturate(PolynomialIdeal(gb,'variables'=X), f); end proc; # simplifies an ideal using a Groebner basis. # you can specify a set of variables or the termorder to use as the second argument Simplify := proc(J::PolynomialIdeal) option system, remember; local ideals, allargs, gens, char, vars, sub, inq, tord, KGB, G, Z, i; allargs := [seq(args[i], i=1..nargs)]; ideals, allargs := selectremove(type, allargs, 'PolynomialIdeal'); if nops(ideals) > 1 then return seq(procname(i,op(allargs)), i=ideals); end if; char := IdealInfo:-Characteristic(J); vars, allargs := selectremove(type, allargs, 'Or'('list'('name'),'set'('name'))); vars := `if`(nops(vars) = 0, IdealInfo:-Variables(J), convert(vars[-1],set)); sub, allargs := selectremove(type, allargs, 'Or'('list'('equation'),'set'('equation'))); sub := convert(map(op, sub), 'set') union convert(select(type, allargs, 'equation'), 'set'); inq, allargs := selectremove(type, allargs, 'Or'('list'('`<>`'),'set'('`<>`'))); inq := convert(map(op, inq), 'set') union convert(select(type, allargs, '`<>`'), 'set'); tord, allargs := selectremove(type, allargs, 'function'); gens := IdealInfo:-Generators(J); if nops(allargs)>0 then error "invalid argument(s)", op(allargs) end if: if nops(sub)=0 and nops(inq)=0 then tord := `if`(nops(tord)=0, 'tord', tord[-1]); G := GroebnerBasis(J, tord, vars); vars := indets(tord, 'name'); KGB := select(proc(a) evalb(indets(a,'name')=vars) end proc, IdealInfo:-KnownGroebnerBases(J)); PolynomialIdeal(G, select(type, [op(J)], 'equation'), 'variables'=vars, 'known_groebner_bases'=[tord=G, seq(i=GroebnerBasis(J,i), i=KGB)]); else gens := gens union {'Z'*mul(lhs(i)-rhs(i), i=inq)-1}; ## the mul evaluates to 1 for empty sequence gens := subs(sub, gens); vars := vars intersect indets(gens,'name'); tord := `if`(nops(tord)=0, 'tdeg'(SuggestVariableOrder(gens,vars)), tord[-1]); vars := indets(tord, 'name'); G := GroebnerBasis(gens, 'prod'(plex('Z'), tord), char); G := remove(has, G, 'Z'); PolynomialIdeal(G, select(type, [op(J)], 'equation'), 'variables'=vars, 'characteristic'=char, 'known_groebner_bases'=[tord=G]); end if; end proc; # homogenize an ideal or a polynomial # maps onto lists and sets Homogenize := proc(f, h::name, vars1) local vars, fh, tord, i; if type(f, 'Or'('list','set')) then map(procname, f, args[2..-1]) elif type(f, 'PolynomialIdeal') then if nargs > 2 and type(vars1, 'set'('name')) then vars := vars1 minus IdealInfo:-NonVariables(f) else vars := IdealInfo:-Variables(f) end if; fh := map(procname, GroebnerBasis(f, 'tord', 'tdeg', vars), args[2..-1]); tord := 'tdeg'(op(tord),h); PolynomialIdeal(fh, 'variables'={op(tord)}, 'characteristic'=IdealInfo:-Characteristic(f), 'known_groebner_bases'=[tord=fh]) else if nargs > 2 and type(vars1, 'set'('name')) then vars := vars1 else vars := indets(f,'name') end if; vars := vars minus indets(indets(f, 'Or'('RootOf','radical','nonreal')), 'name'); subs({seq(i=i/h, i=vars)}, expand(h^degree(f,vars)*f)); end if; end proc: ## ## Algorithms for Polynomial Ideals # these are subroutines called by other programs in the system # they typically operate on a certain type of ideal, or on a Groebner basis # I put them in a separate module mostly to organize the code ideal_algorithms := module() export mis, isfinite, qr_dim, monomial_basis, vectorize, quotient1, quotient1gb, FACTORS, zisprime, zisprimary, zisradical, zradical, zprimedec, zprimarydec; # this module houses programs for computing maximal independent sets mis := module() local dimrec, leadterms, vars, lt, allsets, allindepsets; export maxindepsets; # recursive function for computing independent sets # checks all variables to see what ones can be added # then calls itself recursively on the larger sets dimrec := proc(k, oldvars) local i, newvars; for i from k to nops(vars) do newvars := oldvars union {vars[i]}; if nops(select(type, leadterms, 'polynom'('numeric',newvars))) = 0 then allsets := dimrec(i+1, newvars); end if; end do; for i in allsets do if oldvars intersect i = oldvars then return(allsets); end if; end do; allsets union {oldvars}; end proc; # compute all independent sets of variables allindepsets := proc(gb, tord) if member(1, gb) then return({}); end if; lt := evaluate_termorder(tord); leadterms := {op(map(proc(a) lt(a)[2] end proc, gb))}; vars := indets(tord, 'name'); allsets := {}; dimrec(1, {}); allsets; end proc; # compute all maximal independent sets of variables maxindepsets := proc(gb, tord) option system, remember; local i, s, m; s := allindepsets(gb, tord); m := max(seq(nops(i),i=s)); [op(select(proc(a) evalb(nops(a)=m) end proc, s))]; end proc; end module; # tests whether an ideal generates a finite dimensional vector space isfinite := proc(gb, tord) local lt; lt := evaluate_termorder(tord); evalb(indets(select(type, map(proc(a) lt(a)[2] end proc, gb), 'Or'('name','name'^'posint'))) = indets(tord, 'name') or member(1,gb)); end proc; # computes the dimension of the quotient ring as a vector space qr_dim := proc(gb, tord) local M, N, m, nextmon; if not isfinite(gb,tord) then return infinity; elif member(1,gb) then return 0; end if; nextmon := evaluate_nextmon(plex(op(indets(tord,'name')))); M := map(LeadingMonomial, gb, tord); N := 0; m := 1; while m <> FAIL do N := N+1; m := nextmon(m, M); end do; N; end proc: # computes a monomial basis for the quotient ring as a vector space # returns an unsorted set. you probably want to sort it yourself anyway. # optional arguments # # degree_bound=posint : don't consider monomials above this degree # stop_proc=procedure : halt algorithm if current monomial satisfies this condition # here is how you should sort :) # M := [op(sort(`+`(op(M)), [op(indets(tord,name))], 'plex'))]; monomial_basis := proc(gb, tord) local B, M, m, nextmon, opts, degreebound, stopproc; opts := table(['degree_bound'=infinity, 'stop_proc'=proc(m) option inline; false end proc, op(select(type, [args], 'equation'))]); degreebound := opts['degree_bound']; stopproc := eval(opts['stop_proc']); if nargs=2 and not isfinite(gb,tord) then return FAIL; elif member(1,gb) then return {}; end if; nextmon := evaluate_nextmon(plex(op(indets(tord,'name')))); B := {op(map(LeadingMonomial, gb, tord))}; M := {}; m := 1; while m <> FAIL and not stopproc(m)=FAIL do if degree(m) > degreebound then B := B union {m} else M := M union {m}; end if; m := nextmon(m, B); end do; M; end proc: # convert a polynomial to a vector of coefficients with insane speed # you must give the basis as a list # truncates stuff not in the basis vectorize := proc(f, B::list) local T, C, M, i; C := [coeffs(collect(f, indets(B, 'name'), 'distributed'), indets(B, 'name'), 'M')]; M := [M]; T := table('sparse', [seq(M[i]=C[i], i=1..nops(C))]); [seq(T[i], i=B)] end proc; # compute the quotient J:f quotient1 := proc(G, f, tord, char::nonnegint) local M, i; M := [seq([i,0],i=G), seq([0,i],i=G), [f,1]]; M := ModuleGB(M, tord, 'POT', 'method'='bb', char); M := [seq(`if`(i[1]=0, i[2], NULL), i=M)]; end proc; # compute the quotient J:f where a Groebner basis for J is given quotient1gb := proc(G, f, tord::function, char::nonnegint) local M, i, r; M := [seq([i,0],i=G), seq([0,i],i=G), [f,1]]; M := ModuleGB(M, tord, 'POT', 'method'='bb', 'zap_pairs'=[{seq(i,i=1..nops(G))}, {seq(i,i=nops(G)+1..2*nops(G))}], char); M := [seq(`if`(i[1]=0, i[2], NULL), i=M)]; end proc; # wrapper for polynomial factorization # tries to compensate for deficiencies in `mod/Factors` # input: polynomial, algebraic extensions, characteristic FACTORS := proc(f, K::set, char::nonnegint) local P, K1, dummy; if char=0 then # note: assumes RootOfs are independent (ie: prime decomposition) # so if you use this, you need to factor the system completely # can't be helped, since radfield might fail return factors(f,K) end if; K1 := K union indets(f, {'RootOf','radical'}) union `if`(has(f,I), {RootOf(dummy^2+1)}, {}); P := (Primfield(K1 union {'dummy'}) mod char)[2]; K1 := {op(select(type, map(rhs, P), 'RootOf'))}; if nops(K1) > 1 then error "multiple algebraic extensions could not be reduced to a single extension" else Factors(subs(P,f), op(K1)) mod char; end if end proc: # test if zero-dimensional ideal is prime # algorithm : compute a plex Groebner basis with x1 < x2 < ... < xn # let fi be the polynomial with leading monomial xi^k # factor fi using the roots of {f1,...,f(i-1)}, setting xj = RootOf(fj) for j < i # if it splits, or if it is a power, then ideal is not prime zisprime := proc(J::PolynomialIdeal, K1) local char, vars, tord, G, K, v, F, f; # convert 2nd argument into a set of algebraic extensions K := `if`(nargs=1, {}, `if`(type(K1, 'Or'('list','set')), {op(K1)}, {K1})); K := subs(radfield(K)[1], K); # note: is this a BUG over a finite field ? # ans : no vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); G := GroebnerBasis(J, tord); F := {}; for v in vars do f := select(proc(a) evalb(indets(LeadingMonomial(a, tord),'name')={v}) end proc, G); if nops(f) = 0 then next; else f := FACTORS(subs(F, f[1]), K, char)[2]; f := select(has, f, vars); end if; if nops(f) > 1 or f[1][2] > 1 then # ideal split, or ideal is not radical return false; elif degree(f[1][1], v) > 1 then # add algebraic extension K := K union {RootOf(f[1][1], v)}; end if; # add solution for variable v F := F union {v=RootOf(f[1][1], v)}; end do; true; end proc; # test if zero-dimensional ideal is primary # essentially tests if the radical is prime zisprimary := proc(J::PolynomialIdeal, K1) local char, vars, tord, G, K, v, F, R, f; # convert 2nd argument into a set of algebraic extensions K := `if`(nargs=1, {}, `if`(type(K1, 'Or'('list','set')), {op(K1)}, {K1})); K := subs(radfield(K)[1], K); vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); G := GroebnerBasis(J, tord); F := {}; R := {}; for v in vars do f := select(proc(a) evalb(indets(LeadingMonomial(a, tord))={v}) end proc, G); if nops(f) = 0 then next; else f := FACTORS(subs(F, f[1]), K, char)[2]; f := select(has, f, vars); end if; if nops(f) > 1 then # ideal split return false; elif degree(f[1][1], v) > 1 then # add algebraic extension K := K union {RootOf(f[1][1], v)}; # add substitution to undo the extension R := R union {RootOf(f[1][1], v)=v}; end if; if f[1][2] > 1 then # ideal is not radical add f[1][1] and recompute gb G := GroebnerBasis([op(G), subs(R, f[1][1])], tord, 'method'='bb', char) end if; # add solution for variable v F := F union {v=RootOf(f[1][1], v)}; end do; true; end proc; # test if a zero-dimensional ideal is radical # algorithm : compute a primary decomposition # verify that each factor is radical as you go zisradical := proc(J::PolynomialIdeal) local char, vars, tord, G, C, ISRADICAL, v, ans; vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); # this subroutine does all the work # it is called repeatedly, once per variable, on each component # it returns false if the component is not radical and # factors components when possible ISRADICAL := proc(a) local G, K, F, R, f, i; G, K, F, R := op(a); f := select(proc(a) evalb(indets(LeadingMonomial(a, tord),'name')={v}) end proc, G); if nops(f) = 0 then return a; end if; f := FACTORS(subs(F, f[1]), K, char)[2]; if max(seq(i[2],i=f)) > 1 then return false; elif nops(f) = 1 then [ G, `if`(degree(f[1][1],v) > 1, K union {RootOf(f[1][1], v)}, K), `if`(has(f[1][1],v), F union {v=RootOf(f[1][1], v)}, F), `if`(degree(f[1][1],v) > 1, R union {RootOf(f[1][1], v)=v}, R) ] else seq([ GroebnerBasis([op(G), subs(R, i[1])], tord, 'method'='bb', char), `if`(degree(i[1],v) > 1, K union {RootOf(i[1], v)}, K), `if`(has(i[1],v), F union {v=RootOf(i[1], v)}, F), `if`(degree(i[1],v) > 1, R union {RootOf(i[1], v)=v}, R) ], i=f) end if; end proc: # C is the set of components # elements are of the form [gb, alg_extensions, variable -> rootof subs, rootof -> variable subs] C := {[GroebnerBasis(J, tord), {}, {}, {}]}; ans := {}; for v in vars while not member(false, ans) do ans, C := selectremove(type, map(ISRADICAL, C), 'boolean') end do; not member(false, ans); end proc; # compute the radical of a zero-dimensional ideal # like prime decomposition - but we keep track of # polynomials which must be added to the original ideal zradical := proc(J::PolynomialIdeal) local char, vars, tord, G, C, RADICAL, v, vseq; vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); RADICAL := proc(a) local G, K, F, R, S, P, f, i, j; G, K, F, R, S, P := op(a); f := select(proc(a) evalb(indets(LeadingMonomial(a, tord),'name')={v}) end proc, G); if nops(f) = 0 then return a; end if; f := FACTORS(subs(F, f[1]), K, char)[2]; if nops(f) = 1 and f[1][2] = 1 then [ G, `if`(degree(f[1][1],v) > 1, K union {RootOf(f[1][1], v)}, K), `if`(has(f[1][1],v), F union {v=RootOf(f[1][1], v)}, F), `if`(degree(f[1][1],v) > 1, R union {RootOf(f[1][1], v)=v}, R), S, P ] else seq([ GroebnerBasis([op(G), subs(R, i[1])], tord, 'method'='bb', char), `if`(degree(i[1],v) > 1, K union {RootOf(i[1], v)}, K), `if`(has(i[1],v), F union {v=RootOf(i[1], v)}, F), `if`(degree(i[1],v) > 1, R union {RootOf(i[1], v)=v}, R), S union {v = mul(j[1], j={op(f)} minus {i})}, `if`(i[2] > 1, P union {subs(R, i[1]*subs(S union {v = mul(j[1], j={op(f)} minus {i})}, mul(j, j=subs(vseq, v))))}, P) ], i=f) end if; end proc: G := GroebnerBasis(J, tord); # for each f in G with leading monomial v^i # return a subtitution v = indets(f, 'name') vseq := map(proc(a) local v; v := indets(LeadingMonomial(a, tord), 'name'); if nops(v)=1 then v := op(v); v = indets(a, 'name'); else NULL end if; end proc, G); # C is the set of components, elements are of the form # [gb, alg extensions, variable -> rootof subs, rootof -> variable subs, intersection polys, new polys] # when a square-free part has to be added, the intersection polys are used # to compute a generator which is in the intersection of all the components C := {[G, {}, {}, {}, {}, {}]}; for v in vars do C := map(RADICAL, C) end do; C := map(op, {seq(expand(i[-1]), i=C)}); # get a plex basis for the radical as well # this is generally not very expensive, since we start with a plex basis for J `if`(nops(C)=0, J, PolynomialIdeal(J, C, 'known_groebner_bases'=[tord=GroebnerBasis([op(C),op(G)], tord, 'method'='bb', char)])) end proc; # compute the prime decomposition of a zero-dimensional ideal # algorithm : compute a plex Groebner basis with x1 < x2 < ... < xn # let fi be the polynomial with leading monomial xi^k # factor fi using the roots of {f1,...,f(i-1)}, setting xj = RootOf(fj) for j < i # remove powers, and split the ideal if possible, and recompute the Groebner basis if needed # note to self : this is as good as it gets, stop messing with it zprimedec := proc(J::PolynomialIdeal, K1) local char, vars, tord, K, C, ZPRIMEDEC, v, i, j, rofsF, rofsR; # convert 2nd argument into a set of algebraic extensions K := `if`(nargs=1, {}, `if`(type(K1, 'Or'('list','set')), {op(K1)}, {K1})); rofsF, rofsR := op(radfield(K)[1..2]); K := subs(rofsF, K); vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); ZPRIMEDEC := proc(a) local G, K, F, R, S, f, i; G, K, F, R, S := op(a); f := select(proc(a) evalb(indets(LeadingMonomial(a, tord),'name')={v}) end proc, G); if nops(f) = 0 then return a; else f := FACTORS(subs(F, f[1]), K, char)[2]; f := select(has, f, vars); end if; if nops(f) = 1 and f[1][2] = 1 then [ G, `if`(degree(f[1][1], v) > 1, K union {RootOf(f[1][1], v)}, K), `if`(has(f[1][1],v), F union {v=RootOf(f[1][1], v)}, F), `if`(degree(f[1][1], v) > 1, R union {RootOf(f[1][1], v)=v}, R), S ] else seq([ GroebnerBasis([op(G), subs(R, i[1])], tord, 'method'='bb', char), `if`(degree(i[1], v) > 1, K union {RootOf(i[1], v)}, K), `if`(has(i[1],v), F union {v=RootOf(i[1], v)}, F), `if`(degree(i[1], v) > 1, R union {RootOf(i[1], v)=v}, R), S union {subs(R, i[1])} ], i=f); end if; end proc; # elements of C are of the following form # [groebnerbasis, alg_extensions, variable->rootof subs, rootof->variable subs, new_polynomials_added] C := {[GroebnerBasis(J, tord), K, {}, {}, {}]}; for v in vars do C := map(ZPRIMEDEC, C); end do; seq(PolynomialIdeal(subs(rofsR, i[-1]), J, 'known_groebner_bases'=[tord=subs(rofsR,i[1])]), i=C); end proc; # compute a primary decomposition for a zero-dimensional ideal # algorithm : just like zprimedec, only we don't remove powers zprimarydec := proc(J::PolynomialIdeal, K1) local char, vars, tord, K, C, ZPRIMARYDEC, v, i, j, rofsF, rofsR; # convert 2nd argument into a set of algebraic extensions K := `if`(nargs=1, {}, `if`(type(K1, 'Or'('list','set')), {op(K1)}, {K1})); rofsF, rofsR := op(radfield(K)[1..2]); K := subs(rofsF, K); vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); tord := IdealInfo:-DefaultMonomialOrder(J, 'plex'); vars := ListTools:-Reverse([op(tord)]); ZPRIMARYDEC := proc(a) local G, K, F, R, S, f, i; G, K, F, R, S := op(a); f := select(proc(a) evalb(indets(LeadingMonomial(a, tord),'name')={v}) end proc, G); if nops(f) = 0 then return a; else f := FACTORS(subs(F, f[1]), K, char)[2]; f := select(has, f, vars); end if; if nops(f) = 1 then [ G, `if`(degree(f[1][1], v) > 1, K union {RootOf(f[1][1], v)}, K), `if`(has(f[1][1],v), F union {v=RootOf(f[1][1], v)}, F), `if`(degree(f[1][1], v) > 1, R union {RootOf(f[1][1], v)=v}, R), S ] else seq([ GroebnerBasis([op(G), subs(R, i[1]^i[2])], tord, 'method'='bb', char), `if`(degree(i[1], v) > 1, K union {RootOf(i[1], v)}, K), `if`(has(i[1],v), F union {v=RootOf(i[1], v)}, F), `if`(degree(i[1], v) > 1, R union {RootOf(i[1], v)=v}, R), S union {subs(R, i[1]^i[2])} ], i=f); end if; end proc; # elements of C are of the following form # [groebnerbasis, alg_extensions, variable->rootof subs, rootof->variable subs, new_polynomials_added] C := {[GroebnerBasis(J, tord), K, {}, {}, {}]}; for v in vars do C := map(ZPRIMARYDEC, C); end do; seq(PolynomialIdeal(subs(rofsR, i[-1]), J, 'known_groebner_bases'=[tord=subs(rofsR,i[1])]), i=C); end proc; end module; # compute a maximal independent set of variables MaximalIndependentSet := proc(J::PolynomialIdeal, X::set(name):=IdealInfo:-Variables(J), $) option system, remember; local G, tord, extra_vars, vars, S, vord, val, i; extra_vars := X minus indets(IdealInfo:-Generators(J),'name'); vars := X minus extra_vars; G := GroebnerBasis(J, 'tord', 'tdeg', vars); if ideal_algorithms:-isfinite(G, tord) then return(extra_vars); end if; vord := [SuggestVariableOrder(J, indets(tord, 'name'))]; val := table([seq(vord[i]=2^i, i=1..nops(vord))]); S := ideal_algorithms:-mis:-maxindepsets(G, tord); # ZeroDimensionalDecomposition will need to eliminate all # but these variables so choose the cheapest set possible. # IT IS NOT BACKWARDS, THIS IS RIGHT S := sort(map(a->[a,add(val[i],i=a)], S), (a,b)->evalb(a[2]>b[2])); S[1][1] union extra_vars; end proc; # compute the Hilbert dimension HilbertDimension := proc(J::PolynomialIdeal,X::set(name),$) option system, remember; `if`(IsZeroDimensional(args), 0, nops(MaximalIndependentSet(args))) end proc; # test only if ideal is zero dimensional IsZeroDimensional := proc(J::PolynomialIdeal,X::set(name),$) option system, remember; local gb, tord; gb := GroebnerBasis(J, 'tord', args[2..-1]); ideal_algorithms:-isfinite(gb, tord); end proc; # count the number of solutions over the algebraic closure NumberOfSolutions := proc(J::{PolynomialIdeal, list}, tord1::function, $) local G, tord; if type(J, PolynomialIdeal) then tord := `if`(nargs=1,'tord',tord1); G := GroebnerBasis(J,tord) else G, tord := J, tord1 end if; if ideal_algorithms:-isfinite(G,tord) then ideal_algorithms:-qr_dim(G,tord) else infinity end if; end proc: # compute the set of monomials not reducible by anything in the ideal # optional arguments # # degree_bound=posint : don't consider monomials above this degree # stop_proc=procedure : halt algorithm if current monomial satisfies this condition MonomialBasis := proc(J::{PolynomialIdeal, list}, tord::function) option system, remember; local G, M, lt; G := `if`(type(J,PolynomialIdeal), GroebnerBasis(J, tord), J); M := ideal_algorithms:-monomial_basis(G, tord, args[3..-1]); if M = FAIL then FAIL else lt := evaluate_termorder(tord); sort([op(M)], (a,b)->evalb(lt(a+b)[2]=a)); end if; end proc: # check for proper ideal (not whole ring) IsProper := proc(J::PolynomialIdeal, X::set(name), $) option system, remember; local tord; not member(1, GroebnerBasis(J, 'tord', args[2..-1])); end proc; # test for membership in an ideal IdealMembership := proc(f::Or(radalgfun,list(radalgfun),set(radalgfun)), J::PolynomialIdeal, $) option system, remember; local F, i; if type(f, 'Or'('list','set')) then F := f; else F := {f}; end if; if nops(IdealInfo:-Variables(J) intersect indets(indets(F, 'Or'('RootOf','radical')), 'name')) > 0 then error "the expression to be tested is not polynomial in the ideal variables" end if; for i in F do if NormalForm(i, J, IdealInfo:-DefaultMonomialOrder(J)) <> 0 then return(false); end if; end do; true; end proc; # test a sequence of left to right containments IdealContainment := proc(J::PolynomialIdeal, K::PolynomialIdeal) option system, remember; #local i, tord; if nargs > 2 then procname(args[2..-1]) and procname(J,K); elif IdealInfo:-Characteristic(J) <> IdealInfo:-Characteristic(K) then error "the characteristics are different" elif IdealInfo:-Generators(K) = {} then evalb(IdealInfo:-Generators(J) = {}); elif nops(IdealInfo:-Variables(J) intersect IdealInfo:-NonVariables(K)) > 0 then error "the generators of the first ideal are not polynomial in the variables of the second ideal" else IdealMembership(IdealInfo:-Generators(J), K); end if; end proc; # shortcut for IdealMembership `in` := proc(f, J) if type(f, 'radalgfun') and type(J, 'PolynomialIdeal') then IdealMembership(f, J) else :-`in`(args) end if; end proc; # shortcut for IdealContainment `subset` := proc(S, J) if type(S, 'Or'('set'('radalgfun'), 'PolynomialIdeal')) and type(J, 'PolynomialIdeal') then IdealContainment(PolynomialIdeal(S), J) else :-`subset`(args) end if; end proc; # test for membership in the radical by computing a Groebner basis RadicalMembership := proc(f::radalgfun, J::PolynomialIdeal, $) option system, remember; local vars, char, Z, GB, tord; vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); GB := GroebnerBasis(IdealInfo:-Generators(J) union {1-Z*f}, 'tord', vars union {Z}, char); evalb(GB=[1]); end proc; Add := proc() option system, remember; local user_options, ideals, polys, vars, pars, char, i, unknown; user_options, ideals := selectremove(type, [args], 'equation'); ideals, polys := selectremove(type, ideals, 'PolynomialIdeal'); polys, unknown := selectremove(type, polys, Or(radalgfun,list(radalgfun),set(radalgfun))); if nops(unknown)>0 then error "invalid argument(s)", op(unknown) end if: char := convert(map(IdealInfo:-Characteristic, ideals), 'set'); if nops(char) > 1 then error "ideals with different characteristics can not be added" end if; char := [op(char),0][1]; ideals := [op(ideals), op(map(PolynomialIdeal, polys))]; vars := `union`(seq(IdealInfo:-Variables(i), i=ideals)); pars := `union`(seq(IdealInfo:-Parameters(i), i=ideals)); if nops(vars intersect pars) > 0 then error "the variables of one ideal are parameters of another" end if; PolynomialIdeal(op(ideals), `if`(nops(vars)=0, NULL, 'variables'=vars), 'characteristic'=char); end proc; Multiply := proc() option system, remember; local user_options, ideals, polys, vars, pars, char, gens, i, j, k, unknown; user_options, ideals := selectremove(type, [args], 'equation'); ideals, polys := selectremove(type, ideals, 'PolynomialIdeal'); polys, unknown := selectremove(type, polys, Or(radalgfun,list(radalgfun),set(radalgfun))); if nops(unknown)>0 then error "invalid argument(s)", op(unknown) end if: char := convert(map(IdealInfo:-Characteristic, ideals), 'set'); if nops(char) > 1 then error "ideals with different characteristics can not be multiplied" end if; char := [op(char),0][1]; ideals := [op(ideals), op(map(PolynomialIdeal, polys))]; vars := `union`(seq(IdealInfo:-Variables(i), i=ideals)); pars := `union`(seq(IdealInfo:-Parameters(i), i=ideals)); if nops(vars intersect pars) > 0 then error "the variables of one ideal are parameters of another" end if; ideals := [op(map(IdealInfo:-Generators, ideals))]; gens := {1}; for i in ideals do gens := {seq(seq(j*k, j=gens), k=i)} end do; PolynomialIdeal(gens, `if`(nops(vars)=0, NULL, 'variables'=vars), 'characteristic'=char); end proc; # Compute an ideal quotient # note: Ideal properties are inherited from J Quotient := proc(J::Or(PolynomialIdeal,radalgfun,list(radalgfun),set(radalgfun)), K::Or(PolynomialIdeal,radalgfun,list(radalgfun),set(radalgfun))) option system, remember; local char, vars, pars, F, L, G, tord, i; if type(J, 'PolynomialIdeal') and type(K, 'PolynomialIdeal') then char := convert(map(IdealInfo:-Characteristic, [J,K]), 'set'); if nops(char) > 1 then error "ideals with different characteristics can not be divided" end if; char := [op(char),0][1]; vars := IdealInfo:-Variables(J) union IdealInfo:-Variables(K); pars := IdealInfo:-Parameters(J) union IdealInfo:-Parameters(K); if nops(vars intersect pars) > 0 then error "the variables of one ideal are parameters of another" end if; F := [op(Generators(K))] elif type(J, 'PolynomialIdeal') then char := IdealInfo:-Characteristic(J); vars := IdealInfo:-Variables(J); F := `if`(type(K, 'Or'('list', 'set')), [op(K)], [K]); elif type(K, 'PolynomialIdeal') then char := IdealInfo:-Characteristic(K); vars := IdealInfo:-Variables(K); return procname(PolynomialIdeal(J, 'variables'=vars, 'characteristic'=char), K, args[3..-1]) else return procname(PolynomialIdeal(J), PolynomialIdeal(K), args[3..-1]) end if; if nops(Generators(J))=0 then return J end if; # Groebner basis is known, use a graded basis if possible if nops(IdealInfo:-KnownGroebnerBases(J, vars)) > 0 then tord := IdealInfo:-DefaultMonomialOrder(J, {'tdeg', 'grlex'}); if not member(tord, IdealInfo:-KnownGroebnerBases(J)) then tord := IdealInfo:-DefaultMonomialOrder(J) end if; F := remove(member, F, Generators(J)); G := GroebnerBasis(J, tord); F := remove(member, map(NormalForm, F, G, tord, char), [1]); if nops(F)=0 then return PolynomialIdeal(1, 'characteristic'=char, 'variables'=vars) end if; L := [seq(ideal_algorithms:-quotient1gb(G, i, tord, char), i=F)]; else G := Generators(J); F := remove(member, F, G); if nops(F)=0 then return PolynomialIdeal(1, 'characteristic'=char, 'variables'=vars) end if; tord := 'tdeg'(SuggestVariableOrder(G,vars)); L := [seq(ideal_algorithms:-quotient1(G, i, tord, char), i=F)]; end if; Intersect(seq(PolynomialIdeal(i, 'characteristic'=char, 'variables'=vars, 'known_groebner_bases'=[tord=i]), i=L)) end proc; ZeroDimensionalDecomposition := proc(J::PolynomialIdeal, $) option system, remember; local U, X, tordUX, tordX, V, G, f, s, i, j, S, L; U := MaximalIndependentSet(J); X := IdealInfo:-Variables(J) minus U; if nops(U) = 0 then return(J); elif nops(indets(J, 'Or'('RootOf','radical','nonreal'))) > 0 then # this is a bug i := substitute_rootofs(J); S := [procname(PolynomialIdeal(subs(i[1], J), i[3], 'variables'=(IdealInfo:-Variables(J) union indets(i[3]))))]; S := map(Simplify, subs(i[2], S)); # not safe: assumes Groebner bases are not known return(op(S)); end if; tordUX := select(is_elim_order, IdealInfo:-KnownGroebnerBases(J, X union U), U); tordUX := select(type,tordUX,'lexdeg'('anything','anything')); if nops(tordUX) > 1 then tordUX := tordUX[1]; tordX := 'tdeg'(op(op(1,tordUX))); else V := [SuggestVariableOrder(J, U union X)]; tordUX := 'lexdeg'(select(member, V, X), select(member, V, U)); tordX := 'tdeg'(op(select(member, V, X))); end if; G := GroebnerBasis(J, tordUX); f := lcm(seq(LeadingCoefficient(i, tordX), i=G)); if f = 1 then return(Simplify(J,X)); else f := expand(f); S := {PolynomialIdeal(Saturate(J,f,'s'), 'variables'=X)}; S := S union {procname(PolynomialIdeal(J,f^s))}; end if; # nuke redundant components - this is expensive and absolutely neccesary L := sort([op(S)], length); for i in L do for j in S minus {i} do if IdealContainment(j,i) then S := S minus {i}; break; end if; end do; end do; if nops(S) = 1 then op(S); else op(remove(a->member(1,IdealInfo:-Generators(a)), S)); end if; end proc; # split an ideal into components of different dimension EquidimensionalDecomposition := proc(J::PolynomialIdeal, $) option system, remember; local vars, zdd, d, S, ans, i; vars := IdealInfo:-Variables(J); zdd := map(Contract, [ZeroDimensionalDecomposition(J)], vars); zdd := map(a->[a, HilbertDimension(a,vars)], zdd); zdd := sort(zdd, (a,b)->evalb(a[2] < b[2])); ans := NULL; while nops(zdd) > 0 do d := zdd[1][2]; S, zdd := selectremove(a->evalb(a[2]=d), zdd); ans := ans, Intersect(seq(i[1], i=S)); end do; ans; end proc; # test whether an ideal is radical IsRadical := proc(J::PolynomialIdeal, $) option system, remember; `if`(IsZeroDimensional(J), ideal_algorithms:-zisradical(J), IdealContainment(Radical(J), J)); end proc; # test whether an ideal is maximal (= prime and zero-dimensional) IsMaximal := proc(J::PolynomialIdeal, k::Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))):={}, $) option system, remember; IsZeroDimensional(J) and IsPrime(J,k) end proc; # test whether an ideal is prime IsPrime := proc(J::PolynomialIdeal, k::Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))):={}, $) option system, remember; local char, S; char := IdealInfo:-Characteristic(J); if char <> 0 and (not isprime(char)) then error "this test requires prime characteristic or characteristic zero" elif char <> 0 and nops(IdealInfo:-Parameters(J)) > 0 then error "non-perfect coefficient fields are not supported" end if; if IsZeroDimensional(J) then return IsProper(J) and ideal_algorithms:-zisprime(args); elif char = 0 then S := {ZeroDimensionalDecomposition(J)}; if nops(S) > 1 then # this is a total nightmare, we have to remove embedded components return false; else return IsProper(J) and ideal_algorithms:-zisprime(S[1],k); end if; else error "positive-dimensional ideals in positive characteristic are not supported" end if; end proc; # test whether an ideal is primary IsPrimary := proc(J::PolynomialIdeal, k::Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))):={}, $) option system, remember; local char, S; char := IdealInfo:-Characteristic(J); if char <> 0 and not isprime(char) then error "this test requires prime characteristic or characteristic zero" elif char <> 0 and nops(IdealInfo:-Parameters(J)) > 0 then error "non-perfect coefficient fields are not supported" end if; if IsZeroDimensional(J) then return IsProper(J) and ideal_algorithms:-zisprimary(args); elif char = 0 then S := {ZeroDimensionalDecomposition(J)}; if nops(S) > 1 then return false; else return IsProper(J) and ideal_algorithms:-zisprimary(S[1],k); end if; else error "positive-dimensional ideals in positive characteristic are not supported" end if; end proc; # compute the radical of an ideal Radical := proc(J::PolynomialIdeal, $) option system, remember; local vars, char, R, i; vars := IdealInfo:-Variables(J); char := IdealInfo:-Characteristic(J); if char <> 0 and nops(IdealInfo:-Parameters(J)) > 0 then error "non-perfect coefficient fields are not supported" end if; if IsZeroDimensional(J) then ideal_algorithms:-zradical(J) elif char = 0 then R := {ZeroDimensionalDecomposition(J)}; R := map(ideal_algorithms:-zradical, R); Intersect(seq(Contract(i,vars), i=R)); else error "positive-dimensional ideals in positive characteristic are not supported" end if; end proc; # compute a prime decomposition of the radical of an ideal PrimeDecomposition := proc(J::PolynomialIdeal, k::Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))):={}, $) option system, remember; local char, P; char := IdealInfo:-Characteristic(J); if char <> 0 and (not isprime(char)) then error "this test requires prime characteristic or characteristic zero" elif char <> 0 and nops(IdealInfo:-Parameters(J)) > 0 then error "non-perfect coefficient fields are not supported" end if; if IsZeroDimensional(J) then P := {ideal_algorithms:-zprimedec(args)}; elif char = 0 then P := {ZeroDimensionalDecomposition(J)}; P := map(ideal_algorithms:-zprimedec, P, k); P := select(IsProper, map(Contract, P, IdealInfo:-Variables(J))); else error "positive-dimensional ideals in positive characteristic are not supported" end if; op(P); end proc; # compute a primary decomposition of an ideal PrimaryDecomposition := proc(J::PolynomialIdeal, k::Or(algext,radext,set(Or(algext,radext)), list(Or(algext,radext))):={}, $) option system, remember; local char, P; char := IdealInfo:-Characteristic(J); if char <> 0 and (not isprime(char)) then error "this test requires prime characteristic or characteristic zero" elif char <> 0 and nops(IdealInfo:-Parameters(J)) > 0 then error "non-perfect coefficient fields are not supported" end if; if IsZeroDimensional(J) then return ideal_algorithms:-zprimarydec(args); elif char = 0 then P := {ZeroDimensionalDecomposition(J)}; P := map(ideal_algorithms:-zprimarydec, P, k); P := select(IsProper, map(Contract, P, IdealInfo:-Variables(J))); return `if`(nops(P)=1, J, op(P)); else error "positive-dimensional ideals in positive characteristic are not supported" end if; end proc; # auxilliary procedure for vanishing_ideal PointSort := proc(pointlist,ordervec) local myf; myf := proc(a,b) local j,myval; for j in ordervec while a[j]=b[j] do od: myval := [a[j],b[j]]: evalb(myval = sort(myval)); end proc: return sort(pointlist,myf); end proc: ## compute a vanishing ideal for a list of points ## constructs a monic Groebner basis vanishing_ideal := proc(points, xvars::list(name), torder::function, char) local i, j, k, l, n, m, j0, pos, small, gp0, cvarplace, pointseq, SortedPoints, GG, Peqns, P, g0, g, T, TT, currentvar, OrderVars, f, lm, lmplace; # Remove repeated points pointseq := `if`(char=0, {op(points)}, {op(points)} mod char): if map(nops, pointseq) <> {nops(xvars)} then error "some points have incorrect dimension" end if; n:=nops(pointseq): # n is number of points m:=nops(xvars): # m is number of variables userinfo(2, VanishingIdeal, "number of points: ", n); f := `+`(op(xvars)): OrderVars := NULL: for i to m do lm := LeadingMonomial(f,torder): member(lm, xvars, 'lmplace'): OrderVars := lmplace, OrderVars; f := f - lm: od: OrderVars:= [OrderVars]: SortedPoints := PointSort(pointseq,OrderVars): GG:=[1]: for i to n do j0:=0: P:=SortedPoints[i]: Peqns:= [seq(xvars[j]=P[j],j=1..m)]: for j from 1 to nops(GG) do if char=0 and eval(GG[j],Peqns) <> 0 then j0:=j; break: elif char > 0 and Eval(GG[j],Peqns) mod char <> 0 then j0:=j; break: fi: od: g0:= GG[j0]: # workaround : `mod/Eval` can not handle rational functions # gp0:=`if`(char=0, eval(1/g0,Peqns), 1/(Eval(g0,Peqns) mod char) mod char): gp0:=`if`(char=0, eval(1/g0,Peqns), 1/(eval(g0,Peqns) mod char) mod char): GG := subsop(j0=NULL,GG): # remove g0 from GG # update polynomials in GG if char=0 then GG := [op(1..j0-1,GG), seq(l-eval(l,Peqns)*gp0*g0, l=GG[j0..-1])]; else # workaround : `mod/Eval` can not handle rational functions # GG := [op(1..j0-1,GG), seq(l-modp((Eval(l,Peqns) mod char)*gp0, char)*g0 mod char, l=GG[j0..-1])]; GG := [op(1..j0-1,GG), seq(l-eval(l,Peqns)*gp0*g0 mod char, l=GG[j0..-1])]; end if; T:= LeadingMonomial(g0,torder): for j to m do cvarplace := OrderVars[j]: currentvar := xvars[cvarplace]: TT := currentvar * T: # check to see if LT is divisible by LT of any element in GG for k to nops(GG) while not divide(TT, LeadingMonomial(GG[k],torder)) do od: # if not div. by any poly in GG, then insert NF(g) into GG if k > nops(GG) then g := expand((currentvar - P[cvarplace]) * g0); g := NormalForm(g,GG,torder,char): # find the correct position for the new polynomial pos:=1: small:=0: while (small = 0) and (pos <= nops(GG)) do if TT <> LeadingMonomial(TT + GG[pos], torder) then small:=1: pos:=pos-1: fi: pos:=pos+1: od: GG:=[op( 1..pos-1, GG ), g, op( pos..-1, GG )]: fi: od: od: GG := evala(Normal(GG), {'expanded','independent'}); `if`(char=0, GG, GG mod char); end proc; # construct the vanishing ideal - returns an ideal VanishingIdeal := proc(S::Or(set(list),listlist), X::list(name), T) local G, char, vars, tord; char := `if`(type(args[-1], 'posint'), args[-1], 0); if char > 0 and not isprime(char) then error "characteristic must be zero or prime" end if; tord := `if`(nargs > 2 and type(args[3],'function'), args[3], 'plex'(op(X))); if nargs>2 and not type(args[3],'Or'('nonnegint','function')) then error "expects its 3rd argument, T, to be a valid monomial order, but received %1", T end if: G := vanishing_ideal(S, X, tord, char); vars := {op(X)}; if char=0 then G := map(primpart, G, vars); else G := map(`mod/Primpart`, map(numer@`mod/Normal`, G, char), vars, char) end if; PolynomialIdeal(G, 'variables'=vars, 'characteristic'=char, 'known_groebner_bases'=[tord=G]); end proc; ## provides arithmetic operators for ideals ## note : `in` and `subset` are provided by PolynomialIdeals Operators := module() option package; export Simplify, `+`, `*`, `^`; Simplify := proc(J::PolynomialIdeal) local tord, vars; vars := IdealInfo:-Variables(J); vars := sort(convert(vars, 'list'), 'lexorder'); tord := 'tdeg'(op(vars)); PolynomialIdeals:-Simplify(args, tord) end proc; `+` := overload([ proc(J::PolynomialIdeal, K::Or(PolynomialIdeal, radalgfun)) option overload; Simplify(Add(J, K)) end proc, proc(J::radalgfun, K::PolynomialIdeal) option overload; Simplify(Add(J, K)) end proc ]); `*` := overload([ proc(J::Or(radalgfun, PolynomialIdeal), K::PolynomialIdeal) option overload; Simplify(Multiply(J, K)) end proc, proc(J::Or(radalgfun, PolynomialIdeal), K::PolynomialIdeal^negint) option overload; Simplify(Quotient(J, op(1,K)^(-op(2,K)))) end proc, proc(J::PolynomialIdeal, K::radalgfun) option overload; local n,d,V,tord; n,d := numer(K),denom(K); V := IdealInfo:-Variables(J); if nops(indets(d) intersect V)=0 then Simplify(Multiply(J,K)) elif nops(indets(n) intersect V)=0 then Simplify(Quotient(J, op(1,K)^(-op(2,K)))) else Simplify(Quotient(Multiply(J,n), d)) end if end proc ]); `^` := overload([ proc(J1::PolynomialIdeal, k1::nonnegint) local S, J, k, r; option overload; k := k1; S := Simplify(J1); J := Simplify(Add(J1,1)); while k > 0 do r := irem(k,2,'k'); if r=1 then J := Simplify(Multiply(J,S)) end if; S := Simplify(Multiply(S,S)); end do; J end proc ]); end module: # ideal printing procedure `print/POLYNOMIALIDEAL` := proc() local gens, dummy; gens := remove(type, [args], 'equation'); `if`(nops(gens)=0, subs(dummy=0, ':-`<,>`(dummy)'), subs(dummy=op(gens), ':-`<,>`(dummy)')); end proc: # ideal type checker `type/PolynomialIdeal` := 'specfunc'('anything','POLYNOMIALIDEAL'): # convert to PolynomialIdeal `convert/PolynomialIdeal` := eval(PolynomialIdeal,1): # mod command for PolynomialIdeal # changes ring characteristic `mod/POLYNOMIALIDEAL` := proc() PolynomialIdeal('POLYNOMIALIDEAL'(args[1..-2]), 'characteristic'=args[-1]) end proc: # expands the generators of an ideal `expand/POLYNOMIALIDEAL` := proc() local gens, opts, char, kgb, i; char := 0; kgb := table([]); opts, gens := selectremove(type, [args], 'equation'); hasoption(select(type, opts, 'equation'), 'characteristic'='nonnegint', 'char'); hasoption(select(type, opts, 'equation'), 'known_groebner_bases'='table', 'kgb', 'opts'); kgb := op(op(kgb)); if char=0 then PolynomialIdeal(map(expand, gens), op(opts), 'known_groebner_bases'=kgb) else PolynomialIdeal(map(expand, gens) mod char, op(opts), 'known_groebner_bases'=kgb) end if; end proc: # latex export for ideals `latex/POLYNOMIALIDEAL` := proc() local gens; gens := remove(type, [args], 'equation'); LaTeX:-CountCharacters(2); `\\langle `, `latex/print`(op(gens)), `\\rangle `; end proc: end module: protect('PolynomialIdeals'): #savelib('PolynomialIdeals', '`print/POLYNOMIALIDEAL`', '`type/PolynomialIdeal`', '`convert/PolynomialIdeal`', '`mod/POLYNOMIALIDEAL`', '`expand/POLYNOMIALIDEAL`', '`latex/POLYNOMIALIDEAL`'):