module hilbertp; % Computing Hilbert Polynomial from the Hilbert series.
comment
Authors: H. Michael Moeller, Fernuniversitaet
Hagen, Germany
email: MA105@DHAFEU11.bitnet
H. Melenk, Konrad-Zuse-Zentrum
Berlin, Germany
email: melenk@sc.ZIB-Berlin.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 list(list('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('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('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 list(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;
list ('iplus2,ind2,list('itimes2,list('cadr,dims),
list('iplus2,ind1,-1)));
end;
macro procedure getRarray u;
begin scalar arry,inds;
arry := cadr(u);
inds := cddr u;
if length(inds) = 1 then
return list('getv,list('cdr,arry),car inds)
else
return list('getv,list('cdr,arry),
'functionIndex2 . list('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 list('putv,list('cdr,arry),car inds,val)
else
return list('putv,list('cdr,arry), 'functionIndex2 .
list('car,arry).car inds. cadr inds . NIL , val);
end;
procedure HilbertZeroDimp(nrall, n, rarray);
begin integer i, k, count, vicount;
count := 0;
i := 0;
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 list(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;