1+ #####
2+ with(ListTools):
3+ with(Groebner):
4+
5+
6+
7+
8+
9+ #####
10+ Binomials := module()
11+ description "Some functions for Binomial Ideals";
12+ option package;
13+ #with(ListTools):
14+ #with(Groebner):
15+ export isBinomial, isBinomialGroebner, coefficientsMatrix, monomialCoefficient;
16+
17+
18+
19+
20+
21+ #####Main Functions:
22+
23+ monomialCoefficient := proc(f::polynom,variablesf::list, m::monomial)
24+ # This procedure takes a polynomial f and it returns the coefficient of the monomial m in with respect to the variables variablesf
25+ local monomialsf, TermsList, termsf, coefficientsf, fd, kk, indexm, positionm;
26+ if m = 1 then return (tcoeff(f,variablesf)); end if;
27+ fd := collect(f,variablesf,'distributed');
28+ TermsList := e->`if`(type(e,`+`), convert(e,list), [e]); #taken from Internet
29+ termsf := TermsList(fd);
30+ termsf := collect(termsf,variablesf);
31+ coefficientsf := seq(coeffs(termsf[kk],variablesf),kk=1..(nops(termsf)));
32+ monomialsf := termsf/~[coefficientsf];
33+ positionm := member(m,monomialsf,'indexm');
34+ if (positionm = true) then coeffs(termsf[indexm],variablesf); else 0; end if;
35+ end proc;
36+
37+ coefficientsMatrix := proc(F::list,variablesF::list)
38+ #Input: a list of polynomials
39+ #Output: a matrix C of coefficients and a list m of monomials such that
40+ local vectorF, Fd, TermsList, termsF, kkk, coefficientsF, monomialsF, kk, nrElemsTermsF, nrElemsF;
41+ Fd := collect(F,variablesF,distributed);
42+ TermsList := e->`if`(type(e,`+`), convert(e,list), [e]); #taken from Internet
43+ nrElemsF := nops(F);
44+ termsF := [seq(TermsList(Fd[kkk]),kkk=1..nrElemsF)];
45+ termsF := op~(termsF);
46+ termsF := collect(termsF,variablesF);
47+ nrElemsTermsF := nops(termsF);
48+ coefficientsF := seq(coeffs(termsF[kkk],variablesF),kkk=1..nrElemsTermsF);
49+ monomialsF := termsF/~coefficientsF;convert(%,set);
50+ monomialsF := convert(%,list);
51+ coefficientsF :=
52+ [seq(
53+ [seq(
54+ monomialCoefficient(F[kkk],variablesF,monomialsF[kk])
55+ ,kk=1..nops(monomialsF))]
56+ ,kkk=1..nops(F))];
57+ return (Matrix(coefficientsF), convert(monomialsF,Vector));
58+ end proc;
59+
60+ isBinomial := proc(F::list,variablesF::list)
61+ description "tests binomiality";
62+ # --Checks binomiality of an ideal (here binomial means no monomials are in the ideal)
63+ local isB, kk,A,M,rowDimA,maxNrRowElementsA,minNrRowElementsA;
64+ (A,M) := coefficientsMatrix (F,variablesF);
65+ A := LinearAlgebra[ReducedRowEchelonForm](A);
66+ rowDimA := LinearAlgebra[RowDimension](A);
67+ if (rowDimA = 1) then maxNrRowElementsA := ArrayNumElems(Array(A),NonZero);
68+ else
69+ maxNrRowElementsA := max(seq(ArrayNumElems(Array(A[kk]),NonZero),kk=1..rowDimA));
70+ minNrRowElementsA := min(seq(ArrayNumElems(Array(A[kk]),NonZero),kk=1..rowDimA));
71+ end if;
72+ if maxNrRowElementsA < 2 or minNrRowElementsA = 1 then return (false); end if;
73+ if maxNrRowElementsA = 2 and minNrRowElementsA = 2 then return (true); end if;
74+ isB := isBinomialGroebnerFree(F,variablesF)[1];
75+ if isB then return isB; end if;
76+ if (F = Homogenize(F,ttt,variablesF)) then return isB; end if;
77+ if not isB then
78+ print("Gröbner free method failed.");
79+ end if;
80+ return 0;
81+ end proc;
82+
83+ isBinomialGroebner := proc(F::list,variablesF::list)
84+ #test binomiality using a Groebner Basis
85+ local B, numTermsB, maxNumTermsB, minNumTermsB;
86+ B := Basis(F,tdeg(op(variablesF)));
87+ numTermsB := nops~(B);
88+ maxNumTermsB := max(numTermsB);
89+ minNumTermsB := min(numTermsB);
90+ if maxNumTermsB = 2 and minNumTermsB = 2 then true else false end if
91+ end proc;
92+
93+
94+
95+
96+
97+ #####Local Functions:
98+
99+ local isBinomialGroebnerFree;
100+ isBinomialGroebnerFree := proc(F::list,variablesF::list)
101+ description "test for Binomiality with Conradi-Kahle algorithm";
102+ #Detecting Binomiality using Groebner Free method.
103+ #Implements Algorithm 3.3 in [CK15] and part of Recipe 4.5 in [CK15]
104+ local needToDehom, FF, variablesFF, ttt, isB, bins;
105+ needToDehom := false;
106+ FF := F;
107+ variablesFF := variablesF;
108+ if not (FF = Homogenize(FF,ttt,variablesF)) then
109+ # homogenize the generators of the ideal as in Recipe 4.5 in [CK15]
110+ # The last variable is the homogenization variable
111+ FF := Homogenize(FF,ttt,variablesF);
112+ needToDehom := true;
113+ variablesFF := [op(variablesF),ttt];
114+ end if;
115+ (isB, bins) := CKbasisRecursion (FF,variablesFF);
116+ if not isB then return (false, {}); end if;
117+ if needToDehom then
118+ bins := subs(ttt = 1, bins);
119+ end if;
120+ return (true,bins);
121+ end proc;
122+
123+ local minDegreeElements;
124+ minDegreeElements := proc(F::list,variablesF::list)
125+ local kkkk, degreesList, minDegree, elementsMinDegree, indicesMinDegree;
126+ degreesList := [seq(degree(F[kkkk],variablesF),kkkk=1..nops(F))];
127+ minDegree := min(degreesList);
128+ indicesMinDegree := SearchAll(minDegree, degreesList);
129+ indicesMinDegree := [SearchAll(minDegree, degreesList)];
130+ elementsMinDegree := [seq(F[kkkk], kkkk in indicesMinDegree)];
131+ return elementsMinDegree;
132+ end proc;
133+
134+ local CKbasisRecursion;
135+ CKbasisRecursion := proc(F::list,variablesF::list)
136+ # Taken and adapted to Maple from the Macaulay2 package Binomials, by Thomas Kahle
137+ #Input: A set of standard-graded polynomials,
138+ #Output: (Boolean, a matrix of binomials)
139+ #Algorithm 3.3. in [CK15]
140+ local FF, Fmin, A, M, setFF, setFmin, diffFFmin, maxNrRowElementsA, k, B, setAM, GB, minNrRowElementsA;
141+ if nops (F) = 0 then return (true, Matrix(0)); end if;
142+ B:={};
143+ Fmin:={};
144+ FF := F; # MakeUnique(F);
145+ while (nops(FF) > 0) do
146+ Fmin := minDegreeElements(FF,variablesF);
147+ (A,M) := coefficientsMatrix (Fmin,variablesF);
148+ A := LinearAlgebra[ReducedRowEchelonForm](A);
149+ local rowDimA; rowDimA := LinearAlgebra[RowDimension](A);
150+ if (rowDimA = 1) then maxNrRowElementsA := ArrayNumElems(Array(A),NonZero);
151+ else
152+ maxNrRowElementsA := max(seq(ArrayNumElems(Array(A[k]),NonZero),k=1..rowDimA));
153+ minNrRowElementsA := min(seq(ArrayNumElems(Array(A[k]),NonZero),k=1..rowDimA));
154+ end if;
155+ if maxNrRowElementsA > 2 or minNrRowElementsA = 1 then return (false, {}); end if;
156+ setAM := convert(op~([entries(LinearAlgebra[Multiply](A,M))]),set);
157+ B := B union setAM;
158+ B := remove(member, B, [0]);
159+ setFF := convert(FF,set);
160+ setFmin := convert(Fmin,set);
161+ diffFFmin := setFF minus setFmin;
162+ FF := (convert(FF,set)) minus (convert(Fmin,set));
163+ GB := Basis(Fmin,tdeg(op(variablesF)));
164+ convert([seq( NormalForm(FF[k],GB,tdeg(op(variablesF))), k=1..(nops(FF)) )],set); FF:= convert(%,list);
165+ FF := remove(member, FF, [0]);
166+ if nops(FF) = 0 then return (true, B); end if; ##this avoids the next step in case F is empty
167+ #Check again if the reduction made things zero
168+ if nops(FF) = 0 then return [true, B]; end if; ##this avoids the next step in case F is empty
169+ end do;
170+ # # -- We should never reach this:
171+ print ("You have found a bug in the binomials package. Please report.");
172+ return (0,Matrix[0]);
173+ end proc;
174+
175+
176+
177+
178+
179+ end module; # Binomials
180+
181+
182+
183+
184+
185+ #####Bibliography:
186+ #[CK15] Conradi, Kahle. Detecting Binomiality, 2015.
0 commit comments