Artifact f7d7cd1f7576d4183e0bf3134631947a31f54b92948ad25e8a1c62b320c752db:
- Executable file
r37/packages/groebner/hilberts.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: 13835) [annotate] [blame] [check-ins using] [more...]
module hilberts; % Hilbert series of a set of Monomials. % Author: Joachim Hollman, Royal Institute for Technology, % Stockholm, Sweden % email: <joachim@nada.kth.se> 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); list('times,c,list('expt,v,e)); % -------------- safety check -------------- smacro procedure var!-p(m); and(m,atom(m),not(numberp(m))); smacro procedure check!-expt(m); and(eqcar(m,'expt),var!-p(cadr(m)),numberp(caddr(m))); smacro procedure check!-single!-var(m); if var!-p(m) then t else check!-expt(m); smacro procedure check!-mon(m); if check!-single!-var(m) then t else if eqcar(m,'times) then check!-times(cdr(m)) else nil; smacro procedure check!-args(monl,var); and(listp monl,eqcar(monl,'list),var!-p(var), check!-monl(monl)); symbolic procedure make!-vector (n,pat); begin scalar v; v:=mkvect n; for i:=1:n do putv(v,i,pat); return v; end; % -------------- monomials -------------- smacro procedure alloc!-mon(n); make!-vector(n,0); smacro procedure get!-nth!-exp(mon,n); getv(mon,n); smacro procedure set!-nth!-exp(mon,n,d); putv(mon,n,d); smacro procedure get!-tdeg(mon); getv(mon,0); smacro procedure set!-tdeg(mon,d); putv(mon,0,d); % -------------- ideals -------------- smacro procedure the!-empty!-ideal(); list(nil,nil); smacro procedure get!-next!-mon(ideal); <<x:=caadr ideal; if cdadr ideal then ideal:=list(car ideal,cdadr ideal) else ideal:=the!-empty!-ideal(); x>>; smacro procedure not!-empty!-ideal(ideal); cadr ideal; smacro procedure first!-mon(ideal); caadr ideal; smacro procedure append!-ideals(ideal1,ideal2); list(car ideal2,append(cadr ideal1,cadr ideal2)); symbolic procedure insert!-var(var,ideal); % inserts variable var as last generator of ideal begin scalar last; last:=list(make!-one!-var!-mon(var)); return(list(last,append(cadr ideal,last))); end; symbolic procedure add!-to!-ideal(mon,ideal); % add mon as generator to the ideal begin scalar last; last := list(mon); if ideal = the!-empty!-ideal() 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 check!-args(monl,var) then rerror(groebnr2,18, "Usage: hilb({monomial_1,...,monomial_n} [,variable])"); % return(aeval % list('QUOTIENT, % coefflist2prefix(NPol(gltb2arrideal(monl)), var), % list('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 check!-monl(monl); begin scalar flag,tmp; flag:=t; monl:=gltb!-fix(monl); while monl and flag do << tmp:=car monl; flag:=check!-mon(tmp); monl:=cdr monl; >>; return(flag); end; symbolic procedure check!-times(m); begin scalar flag,tmp; flag:=t; while m and flag do << tmp:= car m; flag:=check!-single!-var(tmp); m:=cdr m; >>; return(flag); end; symbolic procedure coefflist2prefix(cl,var); begin scalar i,poly; i:=0; for each c in cl do << poly:= cons(term(c,var,i),poly); i:=i+1; >>; return (cons('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(list(m),varlist) else if eqcar(m,'expt) then varlist:=union(list(cadr(m)),varlist) else varlist:=union(indets(cdr(m)),varlist); return(varlist); end; symbolic procedure build!-assoc(l); % Given a list of indeterminates (x1 x2 ...xn) we produce % an a-list of the form ((x1.1) (x2.2)... (xn.n)) begin scalar i; i:=0; return(for each var in l collect progn(i:=i+1,cons(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 cons(cons(m,1),monlist) else if eqcar(m,'expt) then cons(cons(cadr(m),caddr(m)),monlist) else (for each x in cdr(m) collect mons!-aux(x)) . monlist; return(monlist); end; symbolic procedure mons!-aux(m); if atom(m) then cons(m,1) else if eqcar(m,'expt) then cons(cadr(m),caddr(m)); 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 build!-assoc() begin scalar tdeg,mon; mon:=alloc!-mon(nvars!*); tdeg:=0; if listp m then for each varnodotexp in m do << set!-nth!-exp(mon,car varnodotexp, cdr varnodotexp); tdeg:=tdeg+cdr varnodotexp; >> else << set!-nth!-exp(mon,car m, cdr m); tdeg:=tdeg+cdr m; >>; set!-tdeg(mon,tdeg); return(mon); end; symbolic procedure gltb!-fix(l); % sometimes GLTB has the form (list (list...)) % instead of (list ...) if and(listp cadr l,caadr(l) = 'list) then cadr l else l; symbolic procedure ge(m1,m2); if get!-tdeg(m1) >= get!-tdeg(m2) then t else nil; symbolic procedure get!-end!-ptr(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(gltb!-fix(xgltb)); nvars!* := length(l); l:= sublis(build!-assoc(l),mons(gltb!-fix(xgltb))); l:=for each m in l collect lmon2arrmon(m); l:=sort(l,'ge); return(list(get!-end!-ptr(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:=next!-var(ideal); if not v then return(base!-case!-pol(ideal)); si:=split!-ideal(ideal,v); return(shift!-add(npol(car si),npol(cadr si))); end; symbolic procedure divides!-by!-var(var,mon); begin scalar div; if get!-nth!-exp(mon,var) = 0 then return(nil); div:=alloc!-mon(nvars!*); for i:=1 : nvars!* do set!-nth!-exp(div,i,get!-nth!-exp(mon,i)); set!-nth!-exp(div,var,get!-nth!-exp(mon,var)-1); set!-tdeg(div,get!-tdeg(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:=alloc!-mon(nvars!*); set!-tdeg(m,d:=get!-tdeg(m2)-get!-tdeg(m1)); while d >= 0 and i <= nvars!* do << set!-nth!-exp(m,i,d:=get!-nth!-exp(m2,i) - get!-nth!-exp(m1,i)); i:= i+1; >>; if d < 0 then return (nil) else return(m); end; symbolic procedure shift!-add(p1,p2); % p1+z*p2 % p1 and p2 are polynomials (nonempty coefficient lists) begin scalar p,pptr; pptr:=p:=cons(car p1,nil); p1:=cdr p1; while p1 and p2 do << rplacd(pptr,cons(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 rem!-mult(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... inew := the!-empty!-ideal(); while not!-empty!-ideal(ipp1) and not!-empty!-ideal(ipp2) do begin if get!-tdeg(first!-mon(ipp2)) < get!-tdeg(first!-mon (ipp1)) then << fmon:=get!-next!-mon(ipp1); isearch:=ipp2; >> else << fmon:=get!-next!-mon(ipp2); isearch:=ipp1; >>; primeflag:=T; while primeflag and not!-empty!-ideal(isearch) do if divides(get!-next!-mon(isearch),fmon) then primeflag:=nil; if primeflag then add!-to!-ideal(fmon,inew); end; if not!-empty!-ideal(ipp1) then return(append!-ideals(inew,ipp1)) else return(append!-ideals(inew,ipp2)); end; symbolic procedure next!-var(ideal); % extracts a variable in the ideal suitable for division begin scalar x,m,var; repeat << m:=get!-next!-mon(ideal); var:=get!-var!-if!-not!-single(m); >> until var or ideal = the!-empty!-ideal(); return(var); end; symbolic procedure get!-var!-if!-not!-single(mon); % returns nil if the monomial is in a single variable, % otherwise the index of the second variable of the monomial begin scalar i,foundvarflag,exp; i := 0; foundvarflag:=nil; while not foundvarflag do << i:=i+1; exp:=get!-nth!-exp(mon,i); if exp > 0 then foundvarflag:=T; >>; foundvarflag:=nil; while i < nvars!* and not foundvarflag do << i:=i+1; exp:=get!-nth!-exp(mon,i); if exp > 0 then foundvarflag:=T; >>; if foundvarflag then return(i) else return(nil); end; symbolic procedure make!-one!-var!-mon(vindex); % returns the monomial consisting of the single variable vindex begin scalar mon; mon := alloc!-mon(nvars!*); for i := 1:nvars!* do set!-nth!-exp(mon,i,0); set!-nth!-exp(mon,vindex,1); set!-tdeg(mon,1); return(mon); end; symbolic procedure split!-ideal(ideal,var); % splits the ideal into two simpler ideals begin scalar div,ideal1,ideal2,m,x; ideal1 := the!-empty!-ideal(); ideal2 := the!-empty!-ideal(); while not!-empty!-ideal(ideal) do << m:=get!-next!-mon(ideal); if div:=divides!-by!-var(var,m) then add!-to!-ideal(div,ideal2) else add!-to!-ideal(m,ideal1); >>; ideal2:=rem!-mult(ideal1,ideal2); ideal1:=insert!-var(var,ideal1); return(list(ideal1,ideal2)); end; symbolic procedure base!-case!-pol(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,tdeg; tdeg:=0; for each mon in cadr ideal do tdeg:=tdeg + get!-tdeg(mon); p:=make!-vector(tdeg,0); putv(p,0,1); degsofar:=0; for each mon in cadr ideal do << e:=get!-tdeg(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;