Artifact af8a9eebbdb6686dbbd0935fe938c4e765d264bc0785dcfb002fd72cab506bea:
- Executable file
r38/packages/groebner/hilbertp.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5940) [annotate] [blame] [check-ins using] [more...]
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;