module hilbertp;% Computing Hilbert Polynomial from the Hilbert series.
comment
Authors: H. Michael Moeller, now Universitaet
Dortmund, Germany
email: moeller@math.uni-dortmund.de
H. Melenk, Konrad-Zuse-Zentrum
Berlin, Germany
email: melenk@zib.de
;
symbolic procedure newhilbi(bas,var,vars);
begin scalar baslt,n,u,grad,h,joa,a,ii,dim0,varx,dmode!*,!*modular;
% Extract leading terms .
baslt:= for each p in cdr bas collect
<< u:=hgspliteval list(p,vars);cadr u >>;
% Replace non atomic elements in the varlist by gensyms .
for each x in cdr vars do
if(pairp x)then
baslt:=cdr subeval {{'equal,x,gensym()},' list . baslt};
% Compute the Hilbertseries .
joa:=hilbsereval list(' list . baslt,var);
% Get the Hilbert polynomial .
grad:=deg(joa,var);
a:=for i:=0 : grad collect coeffn(joa,var,i);
n:= length cdr vars;
% dim0:=( for i:=1 : n product(var + i)) /(for i:=1 : n product i);
varx:=!*a2f var;
dim0:=1;
for i:=1 : n do dim0:=multf(addd(i,varx),dim0);
dim0:=multsq(dim0 ./ 1,1 ./(for i:=1 : n product i));
h:=multsq(car a ./ 1,dim0);
a:=cdr a;
ii:=0;
while a do
<< dim0:=multsq(dim0,addf(varx,numr simp(minus ii))
./ addf(varx,numr simp(n - ii)));
ii:=ii + 1;
if not(car a = 0)then h:=addsq(h,multsq(car a ./ 1,dim0));
a:=cdr a >>;
return mk!*sq h end;
symbolic procedure psnewhilbi u;
begin scalar zz,pl,vl;pl:=reval car u;
if cdr u then vl:=listeval(cadr u,nil);
zz:='list.groebnervars(cdr pl,vl);
return newhilbi(pl,'x,zz)end;
put('hilbertpolynomial,'psopfn,'psnewhilbi);
symbolic procedure hgspliteval pars;
% A variant of Gsplit from grinterf.red.
% Split a polynomial into leading monomial and reductum.
begin scalar vars,x,u,v,w,oldorder,!*factor,!*exp;
integer n,pcount!*;!*exp:=t;
n:=length pars;
u:=reval car pars;
v:=if n > 1 then reval cadr pars else nil;
u:={'list,u};
w:=for each j in groerevlist u
collect if eqexpr j then !*eqn2a j else j;
vars:=groebnervars(w,v);
if not vars then vdperr ' hilbertpolynomial;
oldorder:=vdpinit vars;
w:=a2vdp car w;
if vdpzero!? w then x:=w else
% <<x:=vdpfmon(vdplbc w,vdpevlmon w);
<<x:=vdpfmon('(1 . 1),vdpevlmon w);w:=vdpred w>>;
w:={' list,vdp2a x,vdp2a w};
setkorder oldorder;
return w end;
% Simple Array access method for one- and two-dimensional arrays .
% NO check against misusage is done !
% Usage: Rar:=makeRarray list dim1;Rar:=makeRarray {dim1,dim2};
% val:=getRarray(Rar,ind1);val:=getrarray(Rar,ind1,ind2);
% putRarray(Rar,ind1,val); PutRarray(Rar,in1,ind2,val);
% For two dimensional array access only !
macro procedure functionindex2 u;
begin scalar dims,ind1,ind2;
dims:=cadr u;ind1:=caddr u;ind2:=cadddr u;
return %%%%((ind1 #- 1) #* cadr dims) #+ ind2;
{'iplus2,ind2,{'itimes2,{'cadr,dims},
{'iplus2,ind1,-1}}} end;
macro procedure getrarray u;
begin scalar arry,inds;
arry:=cadr u;inds:=cddr u;
if length inds = 1 then
return {'getv,{'cdr,arry},car inds}
else return {'getv,{'cdr,arry},
'functionIndex2.{'car,arry}.inds} end;
symbolic procedure makerarray dims;
begin scalar u,n;
n:=for each i in dims product i;
u:=mkvect n;return dims . u end;
macro procedure putrarray u;
begin scalar arry,inds,val;
arry:=cadr u;
inds:=cddr u;
val:=nth(u,length u); % PSL: lastcar u;
if length inds = 2 then
return {'putv,{'cdr,arry},car inds,val}
else return {'putv,{'cdr,arry},'functionindex2 .
{' car,arry}.car inds.cadr inds.nil,val} end;
symbolic procedure hilbertzerodimp(nrall,n,rarray);
begin integer i,k,count,vicount;
while(( i:=i+1)<= nrall and count < n)do
begin vicount:=1;
for k:=1 : n do
if(getrarray(rarray,i,k)= 0)then vicount:=vicount + 1;
if vicount = n then count:=count + 1;
end;return count = n end;
symbolic procedure groezerodim!?(f,n);
begin scalar explist,a;integer r;
%explist:= list( vev(lt(f1)),...,vev(lt(fr)));
explist:= for each fi in f collect vdpevlmon fi;
r:= length f;
a:=makerarray {r,n};
for i:=1 step 1 until r do
for k:=1 step 1 until n do
putrarray(a,i,k,nth(nth(explist,i),k));
return hilbertzerodimp(r,n,a)end;
symbolic procedure gzerodimeval u;
begin scalar vl;
if cdr u then vl:=reval cadr u;return gzerodim1(reval car u,vl)end;
put('gzerodim!?,'psopfn,'gzerodimeval);
symbolic procedure gzerodim1(u,v);
begin scalar vars,w,oldorder;
w:=for each j in getrlist u
collect if eqexpr j then !*eqn2a j else j;
if null w then rerror(groebnr2,21,"empty list in hilbertpolynomial");
vars:=groebnervars(w,v);
oldorder:=vdpinit vars;
w:=for each j in w collect f2vdp numr simp j;
w:=groezerodim!?(w,length vars);
setkorder oldorder;
return if w then newhilbi(u,'x,'list.v)else nil end;
symbolic procedure gbtest g;
% Test,if the given set of polynomials is a Groebner basis .
% Only fast to compute plausilbility test .
begin scalar fredu,g1,r,s;
g:=vdplsort g;
% Make abbreviated version of g .
g1:= for each p in g collect
<<r:=vdpred p;
if vdpzero!? r then p else
vdpsum(vdpfmon(vdplbc p,vdpevlmon p),
vdpfmon(vdplbc r,vdpevlmon r))>>;
while g1 do
<<for each p in cdr g1 do
if not groebbuchcrit4t(vdpevlmon car g1,vdpevlmon p)then
<< s:=groebspolynom(car g1,p);
if not vdpzero!? s and
null groebsearchinlist(vdpevlmon s,cddr g1)
then rerror(groebnr2,22,
"****** Not a Groebner basis wrt current ordering")>>;
if groebsearchinlist(vdpevlmon car g1,cdr g1)then fredu:=t;
g1:=cdr g1>>;
if fredu then
<<terpri!* t;
prin2t "WARNING: system is not a fully reduced Groebner basis";
prin2t "with current term ordering">> end;
endmodule;; end;