module hilberts;% Hilbert series of a set of Monomials .
% Author : Joachim Hollman,Royal Institute for Technology,Stockholm,Sweden
% email : < joachim@nada.kth.se >
% Improvement : Herbert Melenk,ZIB Berlin,Takustr 9,email : < melenk@zib.de >
comment
A very brief " description " of the method used.
M=k[x,y,z]/(x^2*y,x*z^2,y^2)
x.
0 --> ker(x.) --> M --> M --> M/x --> 0
M/x = k[x,y,z]/(x^2*y,x*z^2,y^2,x) = k[x,y,z]/(x,y^2)
ker(x.) =((x) +(x^2*y,x*z^2,y^2))/(x^2*y,x*z^2,y^2) =
=(x,y^2)/(x^2*y,x*z^2,y^2)
Hilb(ker(x.)) = Hilb - Hilb
(x,y^2) (x^2*y,x*z^2,y^2)
= 1/(1-t)^3 - Hilb -
k[x,y,z]/(x,y^2)
-(1/(1-t)^3 - Hilb
k[x,y,z]/(x^2*y,x*z^2,y^2)
= Hilb -Hilb
M k[x,y,z]/(x,y^2)
If you only keep the numerator in Hilb = N(t)/(1-t)^3
M
then you get
(1-t)N(t) = N(t) - t(N(t) - N(t) )
I I+(x) I Ann(x) + I
i.e.
N(t) = N(t) + t N(t) (*)
I I+(x) Ann(x) + I
Where
I =(x^2*y,x*z^2,y^2)
I +(x) =(x,y^2)
I + Ann(x) =(x*y,z^2,y^2)
N(t) is the numerator polynomial in Hilb
I k[x,y,z]/I
Equation(*)is what we use to compute the numerator polynomial,i.e.
we " divide out " one variable at a time until we reach a base case.
( One is not limited to single variables but I don't know of any good
strategy for selecting a monomial.)
Usage : hilb({ monomial_1,...,monomial_n } [,variable ]);
fluid '(nvars!*);
% ************** MACROS ETC. **************
smacro procedure term(c,v,e);{ ' times,c,{ ' expt,v,e } };
% -------------- safety check --------------
smacro procedure varp m;and(m,atom(m), not(numberp(m)));
smacro procedure checkexpt m;
eqcar(m,' expt)and varp(cadr(m)) and numberp(caddr(m));
smacro procedure checksinglevar m;
if varp(m)then t else checkexpt(m);
smacro procedure checkmon m;
if checksinglevar(m)then t
else if eqcar(m,' times)then checktimes(cdr(m)) else nil;
smacro procedure checkargs(monl,var);
listp monl and eqcar(monl,' list)and
varp(var)and checkmonl(monl);
symbolic procedure makevector(n,pat);
begin scalar v;v:=mkvect n;
for i:=1:n do putv(v,i,pat);return v end;
% -------------- monomials --------------
smacro procedure allocmon n;makevector(n,0);
smacro procedure getnthexp(mon,n);getv(mon,n);
smacro procedure setnthexp(mon,n,d);putv(mon,n,d);
smacro procedure gettdeg mon;getv(mon,0);
smacro procedure settdeg(mon,d);putv(mon,0, d);
% -------------- ideals --------------
smacro procedure theemptyideal();{ nil,nil };
smacro procedure getnextmon ideal;
<< x:=caadr ideal;
if cdadr ideal then ideal:={ car ideal,cdadr ideal }
else ideal:=theemptyideal();x >>;
smacro procedure notemptyideal ideal;cadr ideal;
smacro procedure firstmon ideal;caadr ideal;
smacro procedure appendideals(ideal1,ideal2);
{ car ideal2,append(cadr ideal1,cadr ideal2)};
symbolic procedure insertvar(var,ideal);
% Inserts variable var as last generator of ideal
begin scalar last;last:={ makeonevarmon(var)};
return({ last,append(cadr ideal,last)})end;
symbolic procedure addtoideal(mon,ideal);
% Add mon as generator to the ideal
begin scalar last;last:={ mon };
if ideal = theemptyideal() then rplaca(cdr(ideal), last)
else rplacd(car(ideal), last);
rplaca(ideal,last)end;
% ************** END OF MACROS ETC. **************
% ************** INTERFACE TO ALGEBRAIC MODE **************
symbolic procedure hilbsereval u;
begin scalar l,monl,var;l:=length u;
if l < 1 or l > 2 then rerror(groebnr2,17,
"Usage: hilb({monomial_1,...,monomial_n} [,variable])")
else if l = 1 then
<< monl:=reval car u;var:=' x >> else
<< monl:= reval car u;var:=reval cadr u >>;
monl:= ' list . for each aa in(cdr monl)collect reval aa;
if not checkargs(monl,var)then rerror(groebnr2,18,
"Usage: hilb({monomial_1,...,monomial_n} [,variable])");
% return(aeval
% {'QUOTIENT,
% coefflist2prefix(NPol(gltb2arrideal(monl)), var),
% {'EXPT,list('PLUS,1,list('TIMES,-1,var)},
% nvars!*)});
return(aeval coefflist2prefix(npol(gltb2arrideal(monl)),var)) end;
% Define "hilb" to be the algebraic mode function
put(' hilb,' psopfn,' hilbsereval);
symbolic procedure checkmonl monl;
begin scalar flag,tmp;flag:=t;monl:=gltbfix(monl);
while monl and flag do
<< tmp:=car monl;
flag:= checkmon(tmp);monl:=cdr monl >>;
return flag end;
symbolic procedure checktimes m;
begin scalar flag,tmp;flag:=t;
while m and flag do
<< tmp:=car m;flag:=checksinglevar(tmp);
m:=cdr m >>;return flag end;
symbolic procedure coefflist2prefix(cl,var);
begin scalar poly;integer i;
for each c in cl do
<< poly:=term(c,var,i). poly;
i:=i + 1 >>;return ' plus . poly end;
symbolic procedure indets l;
% "Indets" returns a list containing all the
% indeterminates of l.
% L is supposed to have a form similar to the variable
% GLTB in the Groebner basis package.
%(LIST(EXPT Z 2)(EXPT X 2) Y)
begin scalar varlist;
for each m in l do
if m neq ' list then
if atom(m) then varlist:=union({ m },varlist)
else if eqcar(m,' expt)then varlist:=union({ cadr(m)},varlist)
else varlist:=union(indets(cdr(m)),varlist);
return varlist end;
symbolic procedure buildassoc l;
% Given a list of indeterminates(x1 x2 ...xn) we produce
% an a-list of the form(( x1 . 1)(x2 . 2)...(xn . n)).
begin integer i;
return(for each var in l collect progn(i:=i #+1,var . i)) end;
symbolic procedure mons l;
% Rewrite the leading monomials(i . e . GLTB).
% the result is a list of monomials of the form :
%(variable . exponent)or(( variable1 . exponent1)...
% (variablen . exponentn))
%
% mons('(LIST(EXPT Z 2)(EXPT X 2)(TIMES Y(EXPT X 3))));
%(((Y . 1)(X . 3))(X . 2)(Z . 2)).
begin scalar monlist;
for each m in l do
if m neq ' list then monlist:=
if atom(m)then(m . 1). monlist
else if eqcar(m,' expt)
then(cadr m . caddr m). monlist
else(for each x in cdr(m)collect monsaux(x)) . monlist;
return monlist end;
symbolic procedure monsaux m;
if eqcar(m,'expt)then cadr m . caddr m else m . 1;
symbolic procedure lmon2arrmon m;
% List-monomial to array-monomial
% a list-monomial has the form:(variable_number . exponent)
% or is a list with entries of this form.
% "variable_number" is the number associated with the variable,
% see buildassoc().
begin scalar mon;integer tdeg;mon:=allocmon nvars!*;
if listp m then
for each varnodotexp in m do
<< setnthexp(mon,car varnodotexp,cdr varnodotexp);
tdeg:=tdeg + cdr varnodotexp >>
else
<< setnthexp(mon,car m,cdr m);tdeg:=tdeg + cdr m >>;
settdeg(mon,tdeg);return mon end;
symbolic procedure gltbfix l;
% Sometimes GLTB has the form(list(list ...))
% instead of(list ...).
if listp cadr l and caadr(l)= ' list then cadr l else l;
symbolic procedure gege(m1,m2);
if gettdeg(m1)>= gettdeg(m2)then t else nil;
symbolic procedure getendptr l;
begin scalar ptr;while l do << ptr:=l;l:=cdr l >>;
return ptr end;
symbolic procedure gltb2arrideal xgltb;
% Convert the monomial ideal given by GLTB(in list form)
% to a list of vectors where each vector represents a monomial.
begin scalar l;l:=indets(gltbfix(xgltb));nvars!*:=length(l);
l:=sublis(buildassoc(l), mons(gltbfix(xgltb)));
l:=for each m in l collect lmon2arrmon(m);
l:=sort(l,' gege);
return { getendptr(l), l } end;
% ************** END OF INTERFACE TO ALGEBRAIC MODE **************
%************** PROCEDURES **************
symbolic procedure npol ideal;
% Recursively computes the numerator of the Hilbert series.
begin scalar v,si;v:=nextvar ideal;
if not v then return basecasepol ideal;
si:=splitideal(ideal,v);
return shiftadd(npol car si,npol cadr si)end;
symbolic procedure dividesbyvar(var,mon);
begin scalar div;if getnthexp(mon,var)= 0 then return nil;
div:=allocmon nvars!*;
for i:=1 : nvars!* do setnthexp(div,i,getnthexp(mon,i));
setnthexp(div,var, getnthexp(mon,var)- 1);
settdeg(div,gettdeg mon - 1);return div end;
symbolic procedure divides(m1,m2);
% Does m1 divide m2?
% m1 and m2 are monomials;
% result: either nil(when m1 does not divide m2)or m2 / m1.
begin scalar m,d,i;i:=1;m:=allocmon(nvars!*);
settdeg(m,d:=gettdeg(m2)- gettdeg(m1));
while d >= 0 and i <= nvars!* do
<< setnthexp(m,i,d:=getnthexp(m2,i)- getnthexp(m1,i));
i:= i+1 >>;
return if d < 0 then nil else m end;
symbolic procedure shiftadd(p1,p2);
% p1 + z * p2;
% p1 and p2 are polynomials(nonempty coefficient lists).
begin scalar p,pptr;pptr:=p:=car p1 . nil;
p1:=cdr p1;
while p1 and p2 do
<< rplacd(pptr,(car p1 + car p2). nil);
p1:=cdr p1;p2:=cdr p2;pptr:=cdr pptr >>;
if p1 then rplacd(pptr,p1)
else rplacd(pptr,p2);return p end;
symbolic procedure remmult(ipp1,ipp2);
% The union of two ideals with redundancy of generators eliminated.
begin scalar fmon,inew,isearch,primeflag,x;
% fix;x is used in the macro...
x:=nil;inew:=theemptyideal();
while notemptyideal(ipp1)and notemptyideal(ipp2)do
begin if gettdeg(firstmon(ipp2)) < gettdeg(firstmon(ipp1))
then << fmon:=getnextmon(ipp1);isearch:=ipp2 >>
else << fmon:=getnextmon(ipp2);isearch:=ipp1 >>;
primeflag:=t;
while primeflag and notemptyideal(isearch)do
if divides(getnextmon(isearch), fmon)then primeflag:=nil;
if primeflag then addtoideal(fmon,inew)end;
return if notemptyideal(ipp1)then appendideals(inew,ipp1)
else appendideals(inew,ipp2)end;
symbolic procedure nextvar ideal;
% Extracts a variable in the ideal suitable for division.
begin scalar m,var,x;x:=nil;
repeat
<< m:=getnextmon ideal;
var:=getvarifnotsingle m;
>> until var or ideal = theemptyideal();
return var end;
symbolic procedure getvarifnotsingle mon;
% Returns nil if the monomial is in a single variable,
% otherwise the index of the second variable of the monomial.
begin scalar foundvarflag,exp;integer i;
while not foundvarflag do
<< i:=i + 1;exp:=getnthexp(mon,i);
if exp > 0 then foundvarflag:=t >>;
foundvarflag:=nil;
while i < nvars!* and not foundvarflag do
<< i:=i + 1;exp:=getnthexp(mon,i);
if exp > 0 then foundvarflag:=t >>;
if foundvarflag then return i else return nil end;
symbolic procedure makeonevarmon vindex;
% Returns the monomial consisting of the single variable vindex.
begin scalar mon;mon:=allocmon nvars!*;
for i:=1 : nvars!* do setnthexp(mon,i,0);
setnthexp(mon,vindex,1);
settdeg(mon,1);return mon end;
symbolic procedure splitideal(ideal,var);
% Splits the ideal into two simpler ideals.
begin scalar div,ideal1,ideal2,m,x;x:=nil;
ideal1:=theemptyideal();ideal2:=theemptyideal();
while notemptyideal(ideal)do
<< m:=getnextmon(ideal);
if div:=dividesbyvar(var,m)then addtoideal(div,ideal2)
else addtoideal(m,ideal1)>>;
ideal2:=remmult(ideal1,ideal2);ideal1:=insertvar(var,ideal1);
return { ideal1,ideal2 } end;
symbolic procedure basecasepol ideal;
% In the base case every monomial is of the form Xi ^ ei;
% result : the numerator polynomial of the Hilbert series
% i.e.(1 - z ^ e1)*(1 - z ^ e2)* ...
begin scalar p,degsofar,e;integer tdeg;
for each mon in cadr ideal do tdeg:=tdeg + gettdeg mon;
p:=makevector(tdeg,0);putv(p,0,1);degsofar:=0;
for each mon in cadr ideal do
<< e:=gettdeg mon;
for j:= degsofar step -1 until 0 do
putv(p,j + e,getv(p,j+e)- getv(p,j));
degsofar:=degsofar + e >>;
return vector2list p end;
symbolic procedure vector2list v;
% Convert a vector v to a list. No type checking is done.
begin scalar u;
for i:=upbv v step -1 until 0 do u:=getv(v,i).u;
return u end;
endmodule;;end;