newPackage("MarkedFamiliesWarsaw", Version => "beta", Date => "Decembre 2022", Authors => { {Name => "Paolo Lella", Email => "paolo.lella@polimi.it", HomePage => "http://www.paololella.it/"} }, Headline => "Computing marked schemes over quasi-stable ideals and modules", PackageImports => {"gfanInterface","Truncations","StronglyStableIdeals","Polyhedra","Elimination"} ) export { -- generic marked sets "supportMarkedSet", "genericMarkedSet", "genericMarkedSetTangentSpace", -- generic marked bases "markedBasisStarReduction", "markedBasisAffineStarReduction", "markedFamily", -- Groebner strata "supportGroebnerSet", "genericGroebnerSet", "genericGroebnerSetTangentSpace", "groebnerStratum", "tangentSpaceGroebnerStratum", "dimTangentSpaceGroebnerStratum", -- Hilbert and Quot schemes "openSubsetHilbertScheme", "openSubsetQuotScheme", -- tangent spaces "tangentSpaceMarkedFamily", "dimTangentSpaceMarkedFamily", "tangentSpaceHilbertScheme", "dimTangentSpaceHilbertScheme", "tangentSpaceQuotScheme", "dimTangentSpaceQuotScheme", -- marked bases and marked syzygies of an ideal "markedBasis", "markedSyzygyMatrix", -- methods quasi-stable ideals "isQuasiStable", "stableSet", -- options "DisplayProgress", "StarReduction", "AffineStarReduction", "TruncationDegree", "GotzmannNumber", "Regularity" } ------------------------------------------------------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------------------------------------------------------ supportMarkedSet = method(TypicalValue=>Sequence, Options=>{Verify=>true}) supportMarkedSet (Matrix) := opts -> M -> ( R := ring M; targetDegrees := for d in (degrees M)#0 list -d#0; r := numRows M; RmodJ := {}; hts := {}; for i from 0 to r - 1 do ( J := trim ideal flatten entries M^{i}; if opts.Verify then ( if not isQuasiStable J then error "expected argument 1 to be a quasi-stable monomial module"; ); RmodJ = append(RmodJ,R/J); JstableSet := stableSet J; hts = hts | ( for ht in JstableSet list map(R^targetDegrees,R^{targetDegrees#i - first degree ht},matrix for j from 0 to r - 1 list if j == i then {ht} else {0_R}) ); ); tails := {}; for ht in hts do ( columnDegree := (degrees ht)#1#0#0; newTail := flatten for i from 0 to r - 1 list for t in rsort flatten entries basis(columnDegree + targetDegrees#i,RmodJ#i) list map(R^targetDegrees,R^{-columnDegree},matrix for j from 0 to r - 1 list if j == i then {sub(t,R)} else {0_R}); tails = append(tails,newTail); ); return (hts,tails); ) supportMarkedSet (Ideal) := opts -> I -> ( return supportMarkedSet(matrix {sort I_*}, Verify=>opts.Verify); ) supportGroebnerSet = method(TypicalValue=>Sequence, Options=>{Verify=>true, MonomialOrder=>null}) supportGroebnerSet (Ideal) := opts -> I -> ( R := ring I; if opts.MonomialOrder =!= null then R = newRing(R,MonomialOrder=>opts.MonomialOrder); J := sub(I,R); if opts.Verify then ( if not isQuasiStable J then error "expected argument 1 to be a quasi-stable monomial module"; ); RmodJ := R/J; hts := stableSet J; tails := {}; for ht in hts do ( newTail := select (for t in flatten entries basis(first degree ht,RmodJ) list sub(t,R), x -> x < ht); tails = append(tails,newTail); ); return (hts,tails); ) ------------------------------------------------------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------------------------------------------------------ genericMarkedSet = method(TypicalValue=>Sequence, Options=>{VariableBaseName=>"A", Verify=>true, DisplayProgress=>true}) genericMarkedSet (Module) := opts -> M -> ( return genericMarkedSet(gens M, VariableBaseName=>opts.VariableBaseName, Verify=>opts.Verify, DisplayProgress=>opts.DisplayProgress); ) genericMarkedSet (Ideal) := opts -> I -> ( (hts,mps) := genericMarkedSet(matrix {sort I_*}, VariableBaseName=>opts.VariableBaseName, Verify=>opts.Verify, DisplayProgress=>opts.DisplayProgress); return markedPolynomialList({for ht in hts list ht_(0,0),for mp in mps list mp_(0,0)}); ) genericMarkedSet (Matrix) := opts -> M -> ( (hts,tails) := supportMarkedSet(M, Verify=>opts.Verify); return genericMarkedSet(hts,tails, VariableBaseName=>opts.VariableBaseName, Verify=>false, DisplayProgress=>opts.DisplayProgress) ) genericMarkedSet (List,List) := opts -> (hts,tails) -> ( R := ring hts#0; if instance(hts#0,RingElement) then ( return genericMarkedSet(for ht in hts list map(R^{0},R^{-first degree ht},matrix {{ht}}), for tail in tails list for t in tail list map(R^{0},R^{-first degree t},matrix {{t}}), VariableBaseName=>opts.VariableBaseName, Verify=>opts.Verify); ); KK := coefficientRing R; n := numgens R; r := numRows hts#0; if opts.Verify then ( M := hts#0; for j from 1 to #hts-1 do M = M | hts#j; for i from 0 to r - 1 do ( gensJ := delete(0_R, flatten entries M^{i}); if #gensJ > 0 then ( J := trim ideal gensJ; if not isQuasiStable J then error "expected argument 1 to be a set of generator of a quasi-stable monomial module"; if sort gensJ != stableSet J then error "expected argument 1 to be a stable set of generator of a quasi-stable monomial module"; ); ); allTails := unique flatten tails; T := allTails#0; for j from 1 to #allTails-1 do T = T | allTails#j; if T%M != T then error "expected argument 2 to be a set of elements not contained in argument 1"; ); weights := optimalWeights(hts,tails); variableWeights := for i from 0 to n - 1 list weights#i; positionWeights := for i from n to n + r - 1 list weights#i; Rw := newRing(R,Degrees=>variableWeights); N := length flatten tails; P := KK[Variables=>N,VariableBaseName=>opts.VariableBaseName]; PR := P[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x -> 0]; k := 0; markedElements := {}; coefficientsWeights := {}; for i from 0 to #hts-1 do ( newElement := sub(hts#i,PR); p := position(flatten entries hts#i, x -> x != 0_R); for t in tails#i do ( q := position(flatten entries t, x -> x != 0_R); coefficientsWeights = append(coefficientsWeights,first degree sub(hts#i_(p,0),Rw) + positionWeights#p - first degree sub(t_(q,0),Rw) - positionWeights#q); newElement = newElement + P_k*sub(t,PR); k=k+1; ); markedElements = append(markedElements,newElement); ); P = newRing(P,MonomialOrder=>{Weights=>coefficientsWeights,RevLex},Local=>true,Degrees=>coefficientsWeights); sortedCoefficients := rsort gens P; sortedCoefficientsWeights := for c in sortedCoefficients list first degree c; P = KK[sortedCoefficients,MonomialOrder=>{Weights=>sortedCoefficientsWeights,Lex},Degrees=>sortedCoefficientsWeights]; PR = P[gens PR,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x -> 0]; return (for ht in hts list sub(ht,PR), for me in markedElements list sub(me,PR)); ) genericMarkedSetTangentSpace = method(TypicalValue=>Sequence, Options=>{VariableBaseName=>"A", Verify=>true, DisplayProgress=>true}) genericMarkedSetTangentSpace (Module) := opts -> M -> ( return genericMarkedSetTangentSpace(gens M, VariableBaseName=>opts.VariableBaseName, Verify=>opts.Verify, DisplayProgress=>opts.DisplayProgress); ) genericMarkedSetTangentSpace (Ideal) := opts -> I -> ( (hts,mps) := genericMarkedSetTangentSpace(matrix {sort I_*}, VariableBaseName=>opts.VariableBaseName, Verify=>opts.Verify, DisplayProgress=>opts.DisplayProgress); return markedPolynomialList({for ht in hts list ht_(0,0),for mp in mps list mp_(0,0)}); ) genericMarkedSetTangentSpace (Matrix) := opts -> M -> ( (hts,tails) := supportMarkedSet(M, Verify=>opts.Verify); return genericMarkedSetTangentSpace(hts,tails, VariableBaseName=>opts.VariableBaseName, Verify=>false, DisplayProgress=>opts.DisplayProgress) ) genericMarkedSetTangentSpace (List,List) := opts -> (hts,tails) -> ( R := ring hts#0; if instance(hts#0,RingElement) then ( return genericMarkedSetTangentSpace(for ht in hts list map(R^{0},R^{-first degree ht},matrix {{ht}}), for tail in tails list for t in tail list map(R^{0},R^{-first degree t},matrix {{t}})); ); KK := coefficientRing R; n := numgens R; r := numRows hts#0; if opts.Verify then ( M := hts#0; for j from 1 to #hts-1 do M = M | hts#j; for i from 0 to r - 1 do ( gensJ := delete(0_R, flatten entries M^{i}); if #gensJ > 0 then ( J := trim ideal gensJ; if not isQuasiStable J then error "expected argument 1 to be a set of generator of a quasi-stable monomial module"; if sort gensJ != stableSet J then error "expected argument 1 to be a stable set of generator of a quasi-stable monomial module"; ); ); allTails := unique flatten tails; T := allTails#0; for j from 1 to #allTails-1 do T = T | allTails#j; if T%M != T then error "expected argument 2 to be a set of elements not contained in argument 1"; ); weights := optimalWeights(hts,tails); variableWeights := for i from 0 to n - 1 list weights#i; positionWeights := for i from n to n + r - 1 list weights#i; Rw := newRing(R,Degrees=>variableWeights); N := length flatten tails; P := KK[Variables=>N,VariableBaseName=>opts.VariableBaseName]; Pe := P["ee"]; Pe = Pe/(ideal(Pe_0^2)); PeR := Pe[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x -> 0]; k := 0; markedElements := {}; coefficientsWeights := {}; for i from 0 to #hts-1 do ( newElement := sub(hts#i,PeR); p := position(flatten entries hts#i, x -> x != 0_R); for t in tails#i do ( q := position(flatten entries t, x -> x != 0_R); coefficientsWeights = append(coefficientsWeights,first degree sub(hts#i_(p,0),Rw) + positionWeights#p - first degree sub(t_(q,0),Rw) - positionWeights#q); newElement = newElement + P_k*Pe_0*sub(t,PeR); k=k+1; ); markedElements = append(markedElements,newElement); ); P = newRing(P,MonomialOrder=>{Weights=>coefficientsWeights,RevLex},Local=>true,Degrees=>coefficientsWeights); sortedCoefficients := rsort gens P; sortedCoefficientsWeights := for c in sortedCoefficients list first degree c; P = KK[sortedCoefficients,MonomialOrder=>{Weights=>sortedCoefficientsWeights,Lex},Degrees=>sortedCoefficientsWeights]; Pe = P[Pe_0]; Pe = Pe/ideal(Pe_0^2); PeR = Pe[gens PeR,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x -> 0]; return (for ht in hts list sub(ht,PeR), for me in markedElements list sub(me,PeR)); ) genericGroebnerSet = method(TypicalValue=>Sequence, Options=>{VariableBaseName=>"A", Verify=>true, MonomialOrder=>null, DisplayProgress=>true}) genericGroebnerSet (Ideal) := opts -> I -> ( (hts,tails) := supportGroebnerSet(I, Verify=>opts.Verify, MonomialOrder=>opts.MonomialOrder); (hts',mps) := genericMarkedSet(hts,tails, VariableBaseName=>opts.VariableBaseName, Verify=>false, DisplayProgress=>opts.DisplayProgress); return markedPolynomialList({for ht in hts' list ht_(0,0),for mp in mps list mp_(0,0)}); ) genericGroebnerSetTangentSpace = method(TypicalValue=>Sequence, Options=>{VariableBaseName=>"T", Verify=>true, MonomialOrder=>null, DisplayProgress=>true}) genericGroebnerSetTangentSpace (Ideal) := opts -> I -> ( (hts,tails) := supportGroebnerSet(I, Verify=>opts.Verify, MonomialOrder=>opts.MonomialOrder); (hts',mps) := genericMarkedSetTangentSpace(hts,tails, VariableBaseName=>opts.VariableBaseName, Verify=>false, DisplayProgress=>opts.DisplayProgress); return markedPolynomialList({for ht in hts' list ht_(0,0),for mp in mps list mp_(0,0)}); ) markedBasisStarReduction = method(TypicalValue=>Sequence,Options=>{Verify=>true,Minimize=>false,Syzygies=>false,DisplayProgress=>false}) markedBasisStarReduction (MarkedPolynomialList) := opts -> MPL -> ( J := ideal MPL#0; R := ring J; n := numgens R - 1; r := #(MPL#0); if opts.Verify then ( if not isQuasiStable J then error "expected argument 1 to be a marked set of polynomials over a quasi-stable ideal"; if numgens J > 0 and sort J_* != stableSet J then error "expected argument 1 to be a marked set of polynomials over a stable set of monomials"; for i from 0 to r-1 do ( T := MPL#1#i - MPL#0#i; if T%J != T then error "expected argument 1 to be a marked set of polynomials over a quasi-stable ideal"; ); ); MF := {}; local T; CI := 51; if opts.DisplayProgress then ( stdio << "\t-- Computation of reduction of S-polynomials" << flush; T = cpuTime(); ); local targetDegrees; local syzygyMatrix; local newSyzygy; if opts.Syzygies then ( targetDegrees = for ht in MPL#0 list -first degree ht; syzygyMatrix = map(R^targetDegrees,R^{},0); ); for i from 0 to #(MPL#0) - 1 do ( ht := MPL#0#i; nonMultiplicativeVariables := sort select(gens R, x -> x > minVar ht); for v in nonMultiplicativeVariables do ( Spolynomial := v*MPL#1#i; if opts.Syzygies then newSyzygy = map(R^targetDegrees,R^{targetDegrees#i - 1},matrix for j from 0 to r - 1 list if j == i then {v} else {0_R}); SpolynomialTerms := terms Spolynomial; while (p := position(SpolynomialTerms, x -> (leadMonomial x)%J == 0)) =!= null do ( (pp,q) := starDecomposition(leadMonomial (SpolynomialTerms#p),MPL#0); Spolynomial = Spolynomial - (leadCoefficient (SpolynomialTerms#p))*q*MPL#1#pp; if opts.Syzygies then newSyzygy = newSyzygy + map(R^targetDegrees,R^{targetDegrees#i - 1},matrix for j from 0 to r - 1 list if j == pp then {- (leadCoefficient (SpolynomialTerms#p))*q} else {0_R}); SpolynomialTerms = terms Spolynomial; ); MF = MF | (for t in SpolynomialTerms list leadCoefficient t); if opts.Syzygies then syzygyMatrix = syzygyMatrix | newSyzygy; if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; P := coefficientRing R; if opts.Minimize then ( KK := coefficientRing P; lexP := KK[rsort gens P,MonomialOrder=>Lex,Degrees=>{numgens P:1}]; MF = sort unique for e in MF list sub(e,lexP); CI = 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; while (p := position (MF, x -> (first degree leadMonomial x) == 1)) =!= null do ( newElimEq := MF#p; nonEliminableCoeff = delete(sub(leadMonomial newElimEq,P),nonEliminableCoeff); newSubst := leadMonomial newElimEq => (leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq); substitutionsCoeff = append(substitutionsCoeff, sub(leadMonomial newElimEq,P) => sub((leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq),P)); MF = delete(0_lexP,sort for e in MF list sub(e,newSubst)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in MPL#1 list sub(mp,substitutionsCoeff); if opts.Syzygies then syzygyMatrix = sub(syzygyMatrix,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; MFideal := ideal()**newP; if #MF > 0 then MFideal = sub(ideal MF,newP); newPR := newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; if opts.Syzygies then ( return (markedPolynomialList({for ht in MPL#0 list sub(ht,newPR),for mp in newMB list sub(mp,newPR)}), MFideal, sub(syzygyMatrix,newPR)); ) else ( return (markedPolynomialList({for ht in MPL#0 list sub(ht,newPR),for mp in newMB list sub(mp,newPR)}), MFideal); ); ) else ( if opts.Syzygies then ( return (MPL,(ideal MF)**P, syzygyMatrix); ) else ( return (MPL,(ideal MF)**P); ); ); ) markedBasisStarReduction (List,List) := opts -> (hts,mes) -> ( MJ := hts#0; for j from 1 to #hts-1 do MJ = MJ | hts#j; R := ring MJ; n := numgens R - 1; r := #hts; if opts.Verify then ( for i from 0 to numRows MJ - 1 do ( J := trim ideal flatten entries MJ^{i}; if not isQuasiStable J then error "expected argument 1 to be a set of generator of a quasi-stable monomial module"; if numgens J > 0 and sort J_* != stableSet J then error "expected argument 1 to be a stable set of generator of a quasi-stable monomial module"; ); for i from 0 to r-1 do ( T := mes#i - hts#i; if T%MJ != T then error "expected argument 2 to be a set of marked elements with head terms described in argument 1"; ); ); MF := {}; local T; CI := 51; if opts.DisplayProgress then ( stdio << "\t-- Computation of reduction of S-polynomials" << flush; T = cpuTime(); ); local targetDegrees; local syzygyMatrix; local newSyzygy; if opts.Syzygies then ( targetDegrees = for ht in hts list -first degree ht; syzygyMatrix = map(R^targetDegrees,R^{},0); ); for i from 0 to r - 1 do ( ht := hts#i; nonMultiplicativeVariables := sort select(gens R, x -> x > minVar ht); for v in nonMultiplicativeVariables do ( Spolynomial := v*mes#i; if opts.Syzygies then newSyzygy = map(R^targetDegrees,R^{targetDegrees#i - 1},matrix for j from 0 to r - 1 list if j == i then {v} else {0_R}); SpolynomialMonomials := monomials Spolynomial; while (p := position(toList(0..numColumns SpolynomialMonomials - 1), x -> (SpolynomialMonomials_{x})%MJ == 0)) =!= null do ( (pp,q) := starDecomposition(SpolynomialMonomials_{p},hts); c := leadCoefficient (Spolynomial//SpolynomialMonomials_{p})_(0,0); Spolynomial = Spolynomial - c*q*mes#pp; if opts.Syzygies then newSyzygy = newSyzygy + map(R^targetDegrees,R^{targetDegrees#i - 1},matrix for j from 0 to r - 1 list if j == pp then {- c*q} else {0_R}); SpolynomialMonomials = monomials Spolynomial; ); MF = MF | (flatten for e in flatten entries Spolynomial list for t in terms e list leadCoefficient t); if opts.Syzygies then syzygyMatrix = syzygyMatrix | newSyzygy; if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; if opts.Minimize then ( P := coefficientRing R; KK := coefficientRing P; lexP := KK[rsort gens P,MonomialOrder=>Lex,Degrees=>{numgens P:1}]; MF = sort unique for e in MF list sub(e,lexP); CI = 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; while (p := position (MF, x -> (first degree leadMonomial x) == 1)) =!= null do ( newElimEq := MF#p; nonEliminableCoeff = delete(sub(leadMonomial newElimEq,P),nonEliminableCoeff); newSubst := leadMonomial newElimEq => (leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq); substitutionsCoeff = append(substitutionsCoeff, sub(leadMonomial newElimEq,P) => sub((leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq),P)); MF = delete(0_lexP,sort for e in MF list sub(e,newSubst)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in mes list sub(mp,substitutionsCoeff); if opts.Syzygies then syzygyMatrix = sub(syzygyMatrix,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; MFideal := ideal()**newP; if #MF > 0 then MFideal = sub(ideal MF,newP); newPR := newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; if opts.Syzygies then ( return (for ht in hts list sub(ht,newPR),for mp in newMB list sub(mp,newPR), MFideal, syzygyMatrix); ) else ( return (for ht in hts list sub(ht,newPR),for mp in newMB list sub(mp,newPR), MFideal); ); ) else ( if opts.Syzygies then ( return (hts,mes,ideal MF, syzygyMatrix); ) else ( return (hts,mes,ideal MF); ); ); ) markedBasisAffineStarReduction = method(TypicalValue => Sequence, Options => { Verify=>true, Minimize=>false, DisplayProgress=>false }) markedBasisAffineStarReduction (MarkedPolynomialList) := opts -> MPL -> ( J := ideal MPL#0; R := ring J; P := coefficientRing R; KK := coefficientRing P; n := numgens R - 1; Jsat := trim sub(J,R_n => 1); local T; CI := 51; if opts.Verify then ( if not isQuasiStable J then error "expected argument 1 to be a marked set of polynomials over a quasi-stable ideal"; if numgens J > 0 and sort J_* != stableSet J then error "expected argument 1 to be a marked set of polynomials over a stable set of monomials"; for i from 0 to #(MPL#0)-1 do ( T := MPL#1#i - MPL#0#i; if T%J != T then error "expected argument 1 to be a marked set of polynomials over a quasi-stable ideal"; ); ); superminimalGenPositions := positions(MPL#0, x -> member(sub(x,R_n=>1),stableSet Jsat)); htsSuperminimal := for p in superminimalGenPositions list sub(MPL#0#p,R_n=>1); mpSuperminimal := for p in superminimalGenPositions list sub(MPL#1#p,R_n=>1); superminimalCoefficients := rsort flatten for i from 0 to #mpSuperminimal-1 list for t in terms (mpSuperminimal#i - htsSuperminimal#i) list leadCoefficient t; superminimalCoefficientsDegrees := for g in superminimalCoefficients list first degree g; nonSuperminimalCoefficients := rsort select (gens P, x -> not member(x,superminimalCoefficients)); nonSuperminimalCoefficientsDegrees := for g in nonSuperminimalCoefficients list first degree g; newP := KK[nonSuperminimalCoefficients | superminimalCoefficients,MonomialOrder=>{Lex=>#nonSuperminimalCoefficientsDegrees,Weights=>superminimalCoefficientsDegrees,Lex},Degrees=>(nonSuperminimalCoefficientsDegrees|superminimalCoefficientsDegrees)]; newPR := newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; htsSuperminimal = for ht in htsSuperminimal list sub(ht,newPR); mpSuperminimal = for mp in mpSuperminimal list sub(mp,newPR); hts := for ht in MPL#0 list sub(ht,newPR); mps := for mp in MPL#1 list sub(mp,newPR); Jsat = sub(Jsat,newPR); superminimalCoefficients = for g in superminimalCoefficients list sub(g,newP); if opts.DisplayProgress then ( stdio << "\t-- Computation of reduction of S-polynomials" << flush; T = cpuTime(); ); MF := {}; for i from 0 to #htsSuperminimal - 1 do ( ht := htsSuperminimal#i; nonMultiplicativeVariables := sort select(gens newPR, x -> x > minVar ht); for v in nonMultiplicativeVariables do ( Spolynomial := v*mpSuperminimal#i; SpolynomialTerms := terms Spolynomial; while (p := position(SpolynomialTerms, x -> (leadMonomial x)%Jsat == 0)) =!= null do ( (pp,q) := starDecomposition(leadMonomial (SpolynomialTerms#p),htsSuperminimal); Spolynomial = Spolynomial - (leadCoefficient (SpolynomialTerms#p))*q*mpSuperminimal#pp; SpolynomialTerms = terms Spolynomial; ); MF = MF | (for t in SpolynomialTerms list leadCoefficient t); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; if opts.DisplayProgress then stdio << "\t-- Computation of normal forms of non-superminimal generators" << flush; CI = 68; MFelimination := {}; for i from 0 to #hts-1 do ( if not member(i,superminimalGenPositions) then ( normalForm := sub(hts#i,newPR_n => 1); normalFormTerms := terms normalForm; while (p := position(normalFormTerms, x -> (leadMonomial x)%Jsat == 0)) =!= null do ( (pp,q) := starDecomposition(leadMonomial (normalFormTerms#p),htsSuperminimal); normalForm = normalForm - (leadCoefficient (normalFormTerms#p))*q*mpSuperminimal#pp; normalFormTerms = terms normalForm; ); for t in terms (sub(mps#i - hts#i, newPR_n => 1) + normalForm) do ( if member(sub(leadMonomial (newEq := leadCoefficient t),P),nonSuperminimalCoefficients) then ( MFelimination = append(MFelimination, newEq); ) else ( MF = append(MF, newEq); ); ); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; if opts.Minimize then ( if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); substitutionsCoeff := {}; for e in MFelimination do ( substitutionsCoeff = append(substitutionsCoeff,leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); lexP := KK[superminimalCoefficients,MonomialOrder=>Lex,Degrees=>{#superminimalCoefficients:1}]; MF = sort unique for e in MF list sub(e,lexP); CI = 45; nonEliminableCoeff := superminimalCoefficients; while (p := position (MF, x -> (first degree leadMonomial x) == 1)) =!= null do ( newElimEq := MF#p; nonEliminableCoeff = delete(sub(leadMonomial newElimEq,newP),superminimalCoefficients); newSubst := leadMonomial newElimEq => (leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq); substitutionsCoeff = append(substitutionsCoeff, sub(leadMonomial newElimEq,newP) => sub((leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq),newP)); MF = delete(0_lexP,sort for e in MF list sub(e,newSubst)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); mps = for mp in mps list sub(mp,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP = KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees,Lex},Degrees=>nonEliminableCoeffDegrees]; MFideal := ideal()**newP; if #MF > 0 then MFideal = sub(ideal MF,newP); newPR = newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return (markedPolynomialList({for ht in hts list sub(ht,newPR),for mp in mps list sub(mp,newPR)}), MFideal); ) else ( return (markedPolynomialList({hts,mps}),ideal (MFelimination | MF)); ); ) markedBasisAffineStarReduction (List,List) := opts -> (hts,mes) -> ( MJ := hts#0; for j from 1 to #hts-1 do MJ = MJ | hts#j; R := ring MJ; P := coefficientRing R; KK := coefficientRing P; n := numgens R - 1; r := numRows MJ; local T; CI := 51; targetDegrees := for d in (degrees MJ)#0 list -d#0; MJsat := map(R^targetDegrees,R^{},0); htsSuperminimal := {}; mesSuperminimal := {}; superminimalGenPositions := {}; for i from 0 to r - 1 do ( J := trim ideal flatten entries MJ^{i}; Jsat := trim sub(J, R_n => 1); if opts.Verify then ( if not isQuasiStable J then error "expected argument 1 to be a set of generator of a quasi-stable monomial module"; if numgens J > 0 and sort J_* != stableSet J then error "expected argument 1 to be a stable set of generator of a quasi-stable monomial module"; ); for g in stableSet Jsat do ( newSatGen := map(R^targetDegrees,R^{targetDegrees#i-first degree g},matrix for j from 0 to r -1 list if j == i then {g} else {0_R}); htsSuperminimal = append(htsSuperminimal,newSatGen); p := position(hts, x-> (sub(x,R_n => 1)//newSatGen)_(0,0) == 1); superminimalGenPositions = append(superminimalGenPositions,p); mesSuperminimal = append(mesSuperminimal,sub(mes#p,R_n => 1)); MJsat = MJsat | newSatGen; ); ); if opts.Verify then ( for i from 0 to r-1 do ( T := mes#i - hts#i; if T%MJ != T then error "expected argument 2 to be a set of marked elements with head terms described in argument 1"; ); ); superminimalCoefficients := rsort flatten flatten for i from 0 to #mesSuperminimal-1 list for e in flatten entries (mesSuperminimal#i - htsSuperminimal#i) list for t in terms e list leadCoefficient t; superminimalCoefficientsDegrees := for g in superminimalCoefficients list first degree g; nonSuperminimalCoefficients := rsort select (gens P, x -> not member(x,superminimalCoefficients)); nonSuperminimalCoefficientsDegrees := for g in nonSuperminimalCoefficients list first degree g; newP := KK[nonSuperminimalCoefficients | superminimalCoefficients,MonomialOrder=>{Lex=>#nonSuperminimalCoefficientsDegrees,Weights=>superminimalCoefficientsDegrees,Lex},Degrees=>(nonSuperminimalCoefficientsDegrees|superminimalCoefficientsDegrees)]; newPR := newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; htsSuperminimal = for ht in htsSuperminimal list sub(ht,newPR); mesSuperminimal = for me in mesSuperminimal list sub(me,newPR); allHts := for ht in hts list sub(ht,newPR); allMes := for me in mes list sub(me,newPR); MJsat = sub(MJsat,newPR); superminimalCoefficients = for g in superminimalCoefficients list sub(g,newP); if opts.DisplayProgress then ( stdio << "\t-- Computation of reduction of S-polynomials" << flush; T = cpuTime(); ); MF := {}; for i from 0 to #htsSuperminimal - 1 do ( ht := htsSuperminimal#i; nonMultiplicativeVariables := sort select(gens newPR, x -> x > minVar ht); for v in nonMultiplicativeVariables do ( Spolynomial := v*mesSuperminimal#i; SpolynomialMonomials := monomials Spolynomial; while (p := position(toList(0..numColumns SpolynomialMonomials - 1), x -> (SpolynomialMonomials_{x})%MJsat == 0)) =!= null do ( (pp,q) := starDecomposition(SpolynomialMonomials_{p},htsSuperminimal); c := leadCoefficient (Spolynomial//SpolynomialMonomials_{p})_(0,0); Spolynomial = Spolynomial - c*q*mesSuperminimal#pp; SpolynomialMonomials = monomials Spolynomial; ); MF = MF | (flatten for e in flatten entries Spolynomial list for t in terms e list leadCoefficient t); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; if opts.DisplayProgress then stdio << "\t-- Computation of normal forms of non-superminimal generators" << flush; CI = 68; MFelimination := {}; for i from 0 to #allHts - 1 do ( if not member(i,superminimalGenPositions) then ( normalForm := sub(allHts#i,newPR_n => 1); normalFormMonomials := monomials normalForm; while (p := position(toList(0..numColumns normalFormMonomials - 1), x -> (normalFormMonomials_{x})%MJsat == 0)) =!= null do ( (pp,q) := starDecomposition(normalFormMonomials_{p},htsSuperminimal); c := leadCoefficient (normalForm//normalFormMonomials_{p})_(0,0); normalForm = normalForm - c*q*mesSuperminimal#pp; normalFormMonomials = monomials normalForm; ); for e in flatten entries (sub(allMes#i - allHts#i, newPR_n => 1) + normalForm) do ( for t in terms e do ( if member(sub(leadMonomial (newEq := leadCoefficient t),P),nonSuperminimalCoefficients) then ( MFelimination = append(MFelimination, newEq); ) else ( MF = append(MF, newEq); ); ); ); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); ); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; if opts.Minimize then ( if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); substitutionsCoeff := {}; for e in MFelimination do ( substitutionsCoeff = append(substitutionsCoeff,leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); lexP := KK[superminimalCoefficients,MonomialOrder=>Lex,Degrees=>{#superminimalCoefficients:1}]; MF = sort unique for e in MF list sub(e,lexP); CI = 45; nonEliminableCoeff := superminimalCoefficients; while (p := position (MF, x -> (first degree leadMonomial x) == 1)) =!= null do ( newElimEq := MF#p; nonEliminableCoeff = delete(sub(leadMonomial newElimEq,newP),superminimalCoefficients); newSubst := leadMonomial newElimEq => (leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq); substitutionsCoeff = append(substitutionsCoeff, sub(leadMonomial newElimEq,newP) => sub((leadCoefficient newElimEq)^(-1)*(leadTerm newElimEq - newElimEq),newP)); MF = delete(0_lexP,sort for e in MF list sub(e,newSubst)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); allMes = for me in allMes list sub(me,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP = KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees,Lex},Degrees=>nonEliminableCoeffDegrees]; MFideal := ideal()**newP; if #MF > 0 then MFideal = sub(ideal MF,newP); newPR = newP[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return (for ht in allHts list sub(ht,newPR),for me in allMes list sub(me,newPR), MFideal); ) else ( return (allHts,allMes,ideal (MFelimination | MF)); ); ) markedFamily = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "A", Algorithm => StarReduction, Syzygies => false, Minimize=>false, DisplayProgress=>false }) markedFamily (Module) := opts -> M -> ( (hts,tails) := genericMarkedSet(M, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); if opts.Algorithm === StarReduction then ( return markedBasisStarReduction (hts, tails, Verify=>false, Minimize => opts.Minimize, Syzygies => opts.Syzygies, DisplayProgress => opts.DisplayProgress); ) else if opts.Algorithm === AffineStarReduction then ( return markedBasisAffineStarReduction (hts, tails, Verify=>false, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( error "unknown Algorithm option"; ); ) markedFamily (Ideal) := opts -> J -> ( MPL := genericMarkedSet( J, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); if opts.Algorithm === StarReduction then ( return markedBasisStarReduction (MPL, Verify=>false, Minimize => opts.Minimize, Syzygies => opts.Syzygies, DisplayProgress => opts.DisplayProgress); ) else if opts.Algorithm === AffineStarReduction then ( return markedBasisAffineStarReduction (MPL, Verify=>false, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( error "unknown Algorithm option"; ); ) tangentSpaceMarkedFamily = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "T", Minimize=>false, DisplayProgress=>false }) tangentSpaceMarkedFamily (Module) := opts -> M -> ( (hts,mes) := genericMarkedSetTangentSpace(M, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); local mf; (hts,mes,mf) = markedBasisStarReduction (hts, mes, Verify=>false, Minimize => false, Syzygies => false, DisplayProgress => opts.DisplayProgress); R := ring M; Pe := coefficientRing ring hts#0; P := coefficientRing Pe; if numgens mf == 0 then ( mf = ideal()**P; ) else ( mf = ideal gens gb ideal for e in mf_* list leadCoefficient e; ); if opts.Minimize then ( KK := coefficientRing P; local T; CI := 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; for e in mf_* do ( nonEliminableCoeff = delete(leadMonomial e,nonEliminableCoeff); substitutionsCoeff = append(substitutionsCoeff, leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in mes list sub(mp,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; newPe := newP[Pe_0]; newPe = newPe/ideal(newPe_0^2); newPeR := newPe[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return (for ht in hts list sub(ht,newPeR),for mp in newMB list sub(mp,newPeR)); ); return (hts,mes,mf); ) tangentSpaceMarkedFamily (Ideal) := opts -> J -> ( MPL := genericMarkedSetTangentSpace( J, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); local mf; (MPL,mf) = markedBasisStarReduction (MPL, Verify=>false, Minimize => false, Syzygies => false, DisplayProgress => opts.DisplayProgress); R := ring J; Pe := coefficientRing ring MPL#0#0; P := coefficientRing Pe; if numgens mf == 0 then ( mf = ideal()**P; ) else ( mf = ideal gens gb ideal for e in mf_* list leadCoefficient e; ); if opts.Minimize then ( KK := coefficientRing P; hts := MPL#0; mes := MPL#1; local T; CI := 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; for e in sort mf_* do ( nonEliminableCoeff = delete(leadMonomial e,nonEliminableCoeff); substitutionsCoeff = append(substitutionsCoeff, leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in mes list sub(mp,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; newPe := newP[Pe_0]; newPe = newPe/ideal(newPe_0^2); newPeR := newPe[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return markedPolynomialList({for ht in hts list sub(ht,newPeR),for mp in newMB list sub(mp,newPeR)}); ); return (MPL,mf); ) tangentSpaceMarkedFamily (Ideal,Ideal) := opts -> (J,I) -> ( MPL0 := markedBasis(J,I); MPL := genericMarkedSetTangentSpace( J, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); PeR := ring MPL#0#0; mesDef := for i from 0 to #(MPL#0)-1 list sub(MPL0#1#i,PeR) + MPL#1#i - MPL#0#i; MPL = markedPolynomialList({MPL#0,mesDef}); local mf; (MPL,mf) = markedBasisStarReduction (MPL, Verify=>false, Minimize => false, Syzygies => false, DisplayProgress => opts.DisplayProgress); mf = ideal gens gb ideal for e in mf_* list leadCoefficient e; if opts.Minimize then ( R := ring J; Pe := coefficientRing ring MPL#0#0; P := ring mf; KK := coefficientRing P; hts := MPL#0; mes := MPL#1; local T; CI := 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; for e in sort mf_* do ( nonEliminableCoeff = delete(leadMonomial e,nonEliminableCoeff); substitutionsCoeff = append(substitutionsCoeff, leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in mes list sub(mp,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; newPe := newP[Pe_0]; newPe = newPe/ideal(newPe_0^2); newPeR := newPe[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return markedPolynomialList({for ht in hts list sub(ht,newPeR),for mp in newMB list sub(mp,newPeR)}); ); return (MPL,mf); ) dimTangentSpaceMarkedFamily = method(TypicalValue => ZZ, Options => { Verify => true, DisplayProgress=>false }) dimTangentSpaceMarkedFamily (Module) := opts -> M -> ( (hts,tails,mf) := tangentSpaceMarkedFamily (M, Verify => opts.Verify, DisplayProgress => opts.DisplayProgress); return (numgens ring mf) - (numgens mf); ) dimTangentSpaceMarkedFamily (Ideal) := opts -> J -> ( (MPL,mf) := tangentSpaceMarkedFamily (J, Verify => opts.Verify, DisplayProgress => opts.DisplayProgress); return (numgens ring mf) - (numgens mf); ) dimTangentSpaceMarkedFamily (Ideal,Ideal) := opts -> (J,I) -> ( (MPL,mf) := tangentSpaceMarkedFamily (J, I, Verify => opts.Verify, DisplayProgress => opts.DisplayProgress); return (numgens ring mf) - (numgens mf); ) groebnerStratum = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "A", MonomialOrder => null, Algorithm => StarReduction, Syzygies => false, Minimize=>false, DisplayProgress=>false }) groebnerStratum (Ideal) := opts -> J -> ( MPL := genericGroebnerSet(J, Verify=>opts.Verify, VariableBaseName => opts.VariableBaseName, MonomialOrder => opts.MonomialOrder, DisplayProgress => opts.DisplayProgress); if opts.Algorithm === StarReduction then ( return markedBasisStarReduction (MPL, Verify=>false, Minimize => opts.Minimize, Syzygies => opts.Syzygies, DisplayProgress => opts.DisplayProgress); ) else if opts.Algorithm === AffineStarReduction then ( return markedBasisAffineStarReduction (MPL, Verify=>false, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( error "unknown Algorithm option"; ); ) tangentSpaceGroebnerStratum = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "T", MonomialOrder => null, Minimize=>false, DisplayProgress=>false }) tangentSpaceGroebnerStratum (Ideal) := opts -> J -> ( MPL := genericGroebnerSetTangentSpace( J, Verify=>opts.Verify, MonomialOrder => opts.MonomialOrder, VariableBaseName => opts.VariableBaseName, DisplayProgress => opts.DisplayProgress); local mf; (MPL,mf) = markedBasisStarReduction (MPL, Verify=>false, Minimize => false, Syzygies => false, DisplayProgress => opts.DisplayProgress); mf = ideal gens gb ideal for e in mf_* list leadCoefficient e; if opts.Minimize then ( R := ring J; Pe := coefficientRing ring MPL#0#0; P := ring mf; KK := coefficientRing P; hts := MPL#0; mes := MPL#1; local T; CI := 45; if opts.DisplayProgress then ( stdio << "\t-- Elimination of redundant parameters" << flush; T = cpuTime(); ); nonEliminableCoeff := gens P; substitutionsCoeff := {}; for e in sort mf_* do ( nonEliminableCoeff = delete(leadMonomial e,nonEliminableCoeff); substitutionsCoeff = append(substitutionsCoeff, leadMonomial e => (leadCoefficient e)^(-1)*(leadTerm e - e)); if opts.DisplayProgress then ( stdio << "." << flush; CI = CI+1; if CI%printWidth == 0 then ( stdio << "\n\t-- " << flush; CI = 0;); ); ); newMB := for mp in mes list sub(mp,substitutionsCoeff); if opts.DisplayProgress then stdio << " completed [" << toString(cpuTime()-T) << " s]\n" << flush; nonEliminableCoeffDegrees := for g in nonEliminableCoeff list first degree g; newP := KK[nonEliminableCoeff,MonomialOrder=>{Weights=>nonEliminableCoeffDegrees},Degrees=>nonEliminableCoeffDegrees]; newPe := newP[Pe_0]; newPe = newPe/ideal(newPe_0^2); newPeR := newPe[gens R,MonomialOrder=>(options R).MonomialOrder,Join=>false,DegreeMap => x-> 0]; return markedPolynomialList({for ht in hts list sub(ht,newPeR),for mp in newMB list sub(mp,newPeR)}); ); return (MPL,mf); ) dimTangentSpaceGroebnerStratum = method(TypicalValue => ZZ, Options => { Verify => true, MonomialOrder => null, DisplayProgress=>false }) dimTangentSpaceGroebnerStratum (Ideal) := opts -> J -> ( (MPL,mf) := tangentSpaceGroebnerStratum (J, Verify => opts.Verify, MonomialOrder => opts.MonomialOrder, DisplayProgress => opts.DisplayProgress); return (numgens ring mf) - (numgens mf); ) tangentSpaceHilbertScheme = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "T", TruncationDegree => null, Minimize => false, DisplayProgress=>false }) tangentSpaceHilbertScheme (Ideal) := opts -> J -> ( R := ring J; n := numgens R; lastButOne := n - 2; local r; local d; local s; if (not isMonomialIdeal J) or (not isQuasiStable J) then ( I := J; inI := ideal leadTerm I; genericChangeOfCoordinates := gens R; while not isQuasiStable inI do ( genericChangeOfCoordinates = flatten entries random(R^{0},R^{n:-1}); changeOfCoordinates := map(R,R,genericChangeOfCoordinates); I = changeOfCoordinates I; inI = ideal leadTerm I; ); if I != J then ( stdio << "\t-- warning: change of coordinates needed: "; for i from 0 to n-1 do stdio << "\n\t-- \t" << toString(R_i) << " => " << toString(genericChangeOfCoordinates#i); stdio << flush; ); if opts.TruncationDegree === GotzmannNumber then ( r = gotzmannNumber hilbertPolynomial inI; inIr := truncate(r,inI); Ir := truncate(r,I); return tangentSpaceMarkedFamily( inIr, Ir, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( stableSetInI := stableSet inI; if opts.TruncationDegree === Regularity then ( d = max for m in stableSetInI list first degree m; inId := truncate(d,inI); Id := truncate(d,I); return tangentSpaceMarkedFamily( inId, Id, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( s = max for m in select(stableSetInI, m -> degree(R_lastButOne,m) > 0) list first degree m; inIs := truncate(s-1,inI); Is := truncate(s-1,I); return tangentSpaceMarkedFamily( inIs, Is, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ); ); ) else ( if opts.TruncationDegree === GotzmannNumber then ( r = gotzmannNumber hilbertPolynomial J; Jr := truncate(r,J); return tangentSpaceMarkedFamily( Jr, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( stableSetJ := stableSet J; if opts.TruncationDegree === Regularity then ( d = max for m in stableSetJ list first degree m; Jd := truncate(d,J); return tangentSpaceMarkedFamily( Jd, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( s = max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Js := truncate(s-1,J); return tangentSpaceMarkedFamily( Js, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ); ); ); ) dimTangentSpaceHilbertScheme = method(TypicalValue => ZZ, Options => { Verify => true, DisplayProgress=>false }) dimTangentSpaceHilbertScheme (Ideal) := opts -> J -> ( R := ring J; n := numgens R; lastButOne := n - 2; local s; if (not isMonomialIdeal J) or (not isQuasiStable J) then ( I := J; inI := ideal leadTerm I; genericChangeOfCoordinates := gens R; while not isQuasiStable inI do ( genericChangeOfCoordinates = flatten entries random(R^{0},R^{n:-1}); changeOfCoordinates := map(R,R,genericChangeOfCoordinates); I = changeOfCoordinates I; inI = ideal leadTerm I; ); stableSetInI := stableSet inI; s = max for m in select(stableSetInI, m -> degree(R_lastButOne,m) > 0) list first degree m; inIs := truncate(s-1,inI); Is := truncate(s-1,I); return dimTangentSpaceMarkedFamily( inIs, Is, Verify => false, DisplayProgress => opts.DisplayProgress); ) else ( stableSetJ := stableSet J; s = max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Js := truncate(s-1,J); return dimTangentSpaceMarkedFamily( Js, Verify => false, DisplayProgress => opts.DisplayProgress); ); ) openSubsetHilbertScheme = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "A", Algorithm => StarReduction, TruncationDegree => null, Syzygies => false, Minimize => false, DisplayProgress=>false }) openSubsetHilbertScheme (Ideal) := opts -> J -> ( if (not isMonomialIdeal J) or (not isQuasiStable J) then error "expected argument 1 to be a quasi-stable monomial ideal"; R := ring J; if opts.TruncationDegree === GotzmannNumber then ( r := gotzmannNumber hilbertPolynomial J; Jr := truncate(r,J); return markedFamily( Jr, Verify => false, VariableBaseName => opts.VariableBaseName, Algorithm => opts.Algorithm, Syzygies => opts.Syzygies, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( stableSetJ := stableSet J; if opts.TruncationDegree === Regularity then ( d := max for m in stableSetJ list first degree m; Jd := truncate(d,J); return markedFamily( Jd, Verify => false, VariableBaseName => opts.VariableBaseName, Algorithm => opts.Algorithm, Syzygies => opts.Syzygies, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) else ( lastButOne := numgens R - 2; s := max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Js := truncate(max(s-1,0),J); return markedFamily( Js, Verify => false, VariableBaseName => opts.VariableBaseName, Algorithm => opts.Algorithm, Syzygies => opts.Syzygies, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ); ); ) openSubsetQuotScheme = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "A", Algorithm => StarReduction, TruncationDegree => null, Syzygies => false, Minimize => false, DisplayProgress=>false }) openSubsetQuotScheme (Module) := opts -> M -> ( R := ring M; gensM := gens M; m := numRows gensM; Mtrunc := (module ideal ())**R; for i from 0 to m - 1 do ( iRow := gensM^{i}; baseDeg := first (degrees iRow)#0#0; J := trim ideal flatten entries iRow; if opts.Verify then ( if (not isMonomialIdeal J) or (not isQuasiStable J) then error "expected argument 1 to be a quasi-stable monomial module"; ); if opts.TruncationDegree === GotzmannNumber then ( r := gotzmannNumber hilbertPolynomial J; Mtrunc = Mtrunc ++ truncate(r+baseDeg, image iRow); ) else ( stableSetJ := stableSet J; if opts.TruncationDegree === Regularity then ( d := max for m in stableSetJ list first degree m; Mtrunc = Mtrunc ++ truncate(d + baseDeg, image iRow); ) else ( lastButOne := numgens R - 2; s := max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Mtrunc = Mtrunc ++ truncate(s - 1 + baseDeg, image iRow); ); ); ); return markedFamily(Mtrunc, Verify => false, VariableBaseName => opts.VariableBaseName, Algorithm => opts.Algorithm, Syzygies => opts.Syzygies, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) tangentSpaceQuotScheme = method(TypicalValue => Sequence, Options => { Verify => true, VariableBaseName => "T", TruncationDegree => null, Minimize => false, DisplayProgress=>false }) tangentSpaceQuotScheme (Module) := opts -> M -> ( R := ring M; gensM := gens M; m := numRows gensM; Mtrunc := (module ideal ())**R; for i from 0 to m - 1 do ( iRow := gensM^{i}; baseDeg := first (degrees iRow)#0#0; J := trim ideal flatten entries iRow; if opts.Verify then ( if (not isMonomialIdeal J) or (not isQuasiStable J) then error "expected argument 1 to be a quasi-stable monomial module"; ); if opts.TruncationDegree === GotzmannNumber then ( r := gotzmannNumber hilbertPolynomial J; Mtrunc = Mtrunc ++ truncate(r+baseDeg, image iRow); ) else ( stableSetJ := stableSet J; if opts.TruncationDegree === Regularity then ( d := max for m in stableSetJ list first degree m; Mtrunc = Mtrunc ++ truncate(d + baseDeg, image iRow); ) else ( lastButOne := numgens R - 2; s := max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Mtrunc = Mtrunc ++ truncate(s - 1 + baseDeg, image iRow); ); ); ); return tangentSpaceMarkedFamily(Mtrunc, Verify => false, VariableBaseName => opts.VariableBaseName, Minimize => opts.Minimize, DisplayProgress => opts.DisplayProgress); ) dimTangentSpaceQuotScheme = method(TypicalValue => ZZ, Options => { Verify => true, DisplayProgress=>false }) dimTangentSpaceQuotScheme (Module) := opts -> M -> ( R := ring M; gensM := gens M; m := numRows gensM; Mtrunc := (module ideal ())**R; for i from 0 to m - 1 do ( iRow := gensM^{i}; baseDeg := first (degrees iRow)#0#0; J := trim ideal flatten entries iRow; if opts.Verify then ( if (not isMonomialIdeal J) or (not isQuasiStable J) then error "expected argument 1 to be a quasi-stable monomial module"; ); stableSetJ := stableSet J; lastButOne := numgens R - 2; s := max for m in select(stableSetJ, m -> degree(R_lastButOne,m) > 0) list first degree m; Mtrunc = Mtrunc ++ truncate(s - 1 + baseDeg, image iRow); ); return dimTangentSpaceMarkedFamily(Mtrunc, Verify => false, DisplayProgress => opts.DisplayProgress); ) isQuasiStable = method(TypicalValue=>Boolean) isQuasiStable (Ideal) := I -> ( if not isMonomialIdeal I then return false; M := max for g in I_* list first degree g; for g in I_* do ( xMin := minVar g; gg := g//xMin; for x in select(gens ring I, y -> y > xMin) do ( d := 1; while d <= M and (x^d*gg)%I != 0 do d = d+1; if d == M+1 then return false; ); ); return true; ) stableSet = method(TypicalValue=>List) stableSet (Ideal) := J -> ( if not isQuasiStable J then error "expected argument 1 to be a quasi-stable ideal"; R := ring J; stabilized := false; G := sort J_*; g := 0; while g < #G do ( xMin := minVar G#g; for y in gens R when y > xMin do ( yg := y*G#g; if (starDecomposition(yg,G) === null) then G = append(G,yg) ); g = g+1; ); return sort G; ) markedBasis = method(TypicalValue=>MarkedPolynomialList) markedBasis (Ideal,Ideal) := (J,I) -> ( if not isQuasiStable J then error "expected argument 1 to be a quasi-stable ideal"; R := ring J; KK := coefficientRing(R); hts := sort stableSet J; D := max for ht in hts list first degree ht; markedPolynomials := {}; for d from 1 to D do ( htsd := sort select(hts, x -> first degree x == d); if #htsd > 0 then ( Rd := flatten entries basis (d,R); Jd := flatten entries super basis(d,J); Id := flatten entries super basis(d,I); if #Jd != #Id then error "expected argument 2 to have a marked basis over argument 1"; M := matrix for g in Id list for x in Rd list sub(g//x,KK); P := for g in Jd list position(Rd, x -> x == g); MP := M_P; if det MP == 0 then error "expected argument 2 to have a marked basis over argument 1"; M = MP^(-1)*M; for ht in htsd do ( p := position(Jd, x -> x == ht); htRow := flatten entries M^{p}; markedPolynomials = append(markedPolynomials,sum(for i from 0 to #htRow-1 list htRow#i*Rd#i)); ); ); ); return markedPolynomialList({hts,markedPolynomials}); ) markedSyzygyMatrix = method(TypicalValue=>Matrix) markedSyzygyMatrix (MarkedPolynomialList) := MPL -> ( hts := MPL#0; J := trim ideal hts; R := ring J; n := numgens R - 1; mp := MPL#1; syzMatrixSupport := {}; for j from 0 to #hts - 1 do ( ht := hts#j; nonMultiplicativaVar := sort select(gens R, x -> x > minVar ht); for v in nonMultiplicativaVar do ( newSyzygy := new MutableList from for g in hts list 0_R; Spolynomial := v*mp#j; newSyzygy#j = v; while Spolynomial != 0_R do ( SpolynomialTerms := terms Spolynomial; p := position(SpolynomialTerms, x -> (leadMonomial x)%J == 0); (pp,q) := starDecomposition(leadMonomial (SpolynomialTerms#p),hts); --pp := position(hts, x -> x == g); Spolynomial = Spolynomial - (leadCoefficient (SpolynomialTerms#p))*q*mp#pp; newSyzygy#pp = newSyzygy#pp - (leadCoefficient (SpolynomialTerms#p))*q; ); syzMatrixSupport = append(syzMatrixSupport, toList newSyzygy); ); ); return matrix entries transpose matrix syzMatrixSupport; ) -------------------------------------------------------------------------------- ---------- UNEXPORTED METHODS ---------- -------------------------------------------------------------------------------- monomial = method(TypicalValue=>RingElement) monomial (Matrix) := M -> ( -- assume M a monomial term return first select(flatten entries M, x -> x != 0); ) minVar = method(TypicalValue=>RingElement) minVar (RingElement) := (m) -> ( -- Assume m a monomial R := ring m; j := numgens R - 1; while j >= 0 and degree(R_j,m) == 0 do j = j-1; if j >= 0 then return R_j; ) minVar (Matrix) := (M) -> ( -- Assume M a monomial term R := ring M; return minVar monomial M; ) maxVar = method(TypicalValue=>RingElement) maxVar (RingElement) := (m) -> ( -- Assume m a monomial R := ring m; j := 0; while j < numgens R and degree(R_j,m) == 0 do j = j+1; if j < numgens R then return R_j; ) maxVar (Matrix) := (M) -> ( -- Assume M a monomial term R := ring M; return maxVar monomial M; ) starDecomposition = method(TypicalValue=>Sequence) starDecomposition (RingElement,List) := (m,G) -> ( -- Assume m a monomial for i from 0 to #G-1 do ( --if m == G#i then return (i,1); if (q := m//(G#i)) != 0 and ( q == 1 or minVar(G#i) >= maxVar(q) ) then return (i,q); ); return null; ) starDecomposition (Matrix,List) := (m,G) -> ( -- Assume m a monomial term for i from 0 to #G-1 do ( if (q := m//(G#i)) != 0 and ( q_(0,0) == 1 or minVar(G#i) >= maxVar(q) ) then return (i,q_(0,0)); ); return null; ) optimalWeights = method(TypicalValue=>List) optimalWeights (List,List) := (hts,tails) -> ( R := ring hts#0; n := numgens R; inequalities := {}; if instance(hts#0,RingElement) then ( inequalities = for i from 0 to n-2 list for j from 0 to n-1 list if j == i then 1 else if j == i+1 then -1 else 0; inequalities = append(inequalities, (for i from 0 to n-2 list 0) | {1}); for j from 0 to #hts-1 do for t in tails#j do inequalities = unique append(inequalities, for i from 0 to n-1 list degree(R_i,hts#j) - degree(R_i,t)); ) else ( r := numRows hts#0; inequalities = for i from 0 to n-2 list for j from 0 to n+r-1 list if j == i then 1 else if j == i+1 then -1 else 0; inequalities = append(inequalities, (for i from 0 to n-2 list 0) | {1} | (for i from 1 to r list 0)); for j from 0 to #hts-1 do ( p := position(flatten entries hts#j, x -> x != 0_R); for t in tails#j do ( q := position(flatten entries t, x -> x != 0_R); newInequality := for i from 0 to n-1 list degree(R_i,hts#j_(p,0)) - degree(R_i,t_(q,0)); newInequality = newInequality | (for i from 0 to r-1 list if i == p and i != q then 1 else if i == q and i != p then -1 else 0); inequalities = unique append(inequalities, newInequality); ); ); ); segmentCone := coneFromHData matrix inequalities; extremalRays := entries transpose rays segmentCone; interiorPoint := sum extremalRays; return interiorPoint; ) monomialTerms = method(TypicalValue=>List) monomialTerms (Matrix) := M -> ( (T,C) := coefficients )