1module hilberts;% Hilbert series of a set of Monomials . 2 3% Author : Joachim Hollman,Royal Institute for Technology,Stockholm,Sweden 4% email : < joachim@nada.kth.se > 5% Improvement : Herbert Melenk,ZIB Berlin,Takustr 9,email : < melenk@zib.de > 6 7% Redistribution and use in source and binary forms, with or without 8% modification, are permitted provided that the following conditions are met: 9% 10% * Redistributions of source code must retain the relevant copyright 11% notice, this list of conditions and the following disclaimer. 12% * Redistributions in binary form must reproduce the above copyright 13% notice, this list of conditions and the following disclaimer in the 14% documentation and/or other materials provided with the distribution. 15% 16% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 17% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 18% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 19% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 20% CONTRIBUTORS 21% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27% POSSIBILITY OF SUCH DAMAGE. 28% 29 30 31COMMENT 32 33A very brief " description " of the method used. 34 35M=k[x,y,z]/(x^2*y,x*z^2,y^2) 36 x. 370 --> ker(x.) --> M --> M --> M/x --> 0 38 39M/x = k[x,y,z]/(x^2*y,x*z^2,y^2,x) = k[x,y,z]/(x,y^2) 40 41ker(x.) =((x) +(x^2*y,x*z^2,y^2))/(x^2*y,x*z^2,y^2) = 42 43 =(x,y^2)/(x^2*y,x*z^2,y^2) 44 45Hilb(ker(x.)) = Hilb - Hilb 46 (x,y^2) (x^2*y,x*z^2,y^2) 47 48 = 1/(1-t)^3 - Hilb - 49 k[x,y,z]/(x,y^2) 50 51 -(1/(1-t)^3 - Hilb 52 k[x,y,z]/(x^2*y,x*z^2,y^2) 53 54 = Hilb -Hilb 55 M k[x,y,z]/(x,y^2) 56 57If you only keep the numerator in Hilb = N(t)/(1-t)^3 58 M 59then you get 60 61(1-t)N(t) = N(t) - t(N(t) - N(t) ) 62 I I+(x) I Ann(x) + I 63 64i.e. 65 66 N(t) = N(t) + t N(t) (*) 67 I I+(x) Ann(x) + I 68 69Where 70 I =(x^2*y,x*z^2,y^2) 71 I +(x) =(x,y^2) 72 I + Ann(x) =(x*y,z^2,y^2) 73 N(t) is the numerator polynomial in Hilb 74 I k[x,y,z]/I 75 76Equation(*)is what we use to compute the numerator polynomial,i.e. 77we " divide out " one variable at a time until we reach a base case. 78( One is not limited to single variables but I don't know of any good 79strategy for selecting a monomial.) 80 81Usage : hilb({ monomial_1,...,monomial_n } [,variable ]); 82 83fluid '(nvars!*); 84 85% ************** MACROS ETC. ************** 86 87smacro procedure term(c,v,e);{ ' times,c,{ ' expt,v,e } }; 88 89% -------------- safety check -------------- 90 91smacro procedure varp m;and(m,atom(m), not(numberp(m))); 92 93smacro procedure checkexpt m; 94 eqcar(m,' expt)and varp(cadr(m)) and numberp(caddr(m)); 95 96smacro procedure checksinglevar m; 97 if varp(m)then t else checkexpt(m); 98 99smacro procedure checkmon m; 100 if checksinglevar(m)then t 101 else if eqcar(m,' times)then checktimes(cdr(m)) else nil; 102 103smacro procedure checkargs(monl,var); 104 listp monl and eqcar(monl,' list)and 105 varp(var)and checkmonl(monl); 106 107symbolic procedure makevector(n,pat); 108begin scalar v;v:=mkvect n; 109 for i:=1:n do putv(v,i,pat);return v end; 110 111% -------------- monomials -------------- 112 113smacro procedure allocmon n;makevector(n,0); 114 115smacro procedure getnthexp(mon,n);getv(mon,n); 116 117smacro procedure setnthexp(mon,n,d);putv(mon,n,d); 118 119smacro procedure gettdeg mon;getv(mon,0); 120 121smacro procedure settdeg(mon,d);putv(mon,0, d); 122 123% -------------- ideals -------------- 124 125smacro procedure theemptyideal();{ nil,nil }; 126 127smacro procedure getnextmon ideal; 128<< x:=caadr ideal; 129 if cdadr ideal then ideal:={ car ideal,cdadr ideal } 130 else ideal:=theemptyideal();x >>; 131 132smacro procedure notemptyideal ideal;cadr ideal; 133 134smacro procedure firstmon ideal;caadr ideal; 135 136smacro procedure appendideals(ideal1,ideal2); 137{ car ideal2,append(cadr ideal1,cadr ideal2)}; 138 139symbolic procedure insertvar(var,ideal); 140% Inserts variable var as last generator of ideal 141begin scalar last;last:={ makeonevarmon(var)}; 142 return({ last,append(cadr ideal,last)})end; 143 144symbolic procedure addtoideal(mon,ideal); 145% Add mon as generator to the ideal 146begin scalar last;last:={ mon }; 147 if ideal = theemptyideal() then rplaca(cdr(ideal), last) 148 else rplacd(car(ideal), last); 149 rplaca(ideal,last)end; 150 151% ************** END OF MACROS ETC. ************** 152 153% ************** INTERFACE TO ALGEBRAIC MODE ************** 154 155symbolic procedure hilbsereval u; 156begin scalar l,monl,var;l:=length u; 157 if l < 1 or l > 2 then rerror(groebnr2,17, 158 "Usage: hilb({monomial_1,...,monomial_n} [,variable])") 159 else if l = 1 then 160 << monl:=reval car u;var:=' x >> else 161 << monl:= reval car u;var:=reval cadr u >>; 162 monl:= ' list . for each aa in(cdr monl)collect reval aa; 163 if not checkargs(monl,var)then rerror(groebnr2,18, 164 "Usage: hilb({monomial_1,...,monomial_n} [,variable])"); 165% return(aeval 166% {'QUOTIENT, 167% coefflist2prefix(NPol(gltb2arrideal(monl)), var), 168% {'EXPT,list('PLUS,1,list('TIMES,-1,var)}, 169% nvars!*)}); 170 return(aeval coefflist2prefix(npol(gltb2arrideal(monl)),var)) end; 171 172% Define "hilb" to be the algebraic mode function 173put(' hilb,' psopfn,' hilbsereval); 174 175symbolic procedure checkmonl monl; 176begin scalar flag,tmp;flag:=t;monl:=gltbfix(monl); 177 while monl and flag do 178 << tmp:=car monl; 179 flag:= checkmon(tmp);monl:=cdr monl >>; 180 return flag end; 181 182symbolic procedure checktimes m; 183begin scalar flag,tmp;flag:=t; 184 while m and flag do 185 << tmp:=car m;flag:=checksinglevar(tmp); 186 m:=cdr m >>;return flag end; 187 188symbolic procedure coefflist2prefix(cl,var); 189begin scalar poly;integer i; 190 for each c in cl do 191 << poly:=term(c,var,i). poly; 192 i:=i + 1 >>;return ' plus . poly end; 193 194symbolic procedure indets l; 195% "Indets" returns a list containing all the 196% indeterminates of l. 197% L is supposed to have a form similar to the variable 198% GLTB in the Groebner basis package. 199%(LIST(EXPT Z 2)(EXPT X 2) Y) 200begin scalar varlist; 201 for each m in l do 202 if m neq ' list then 203 if atom(m) then varlist:=union({ m },varlist) 204 else if eqcar(m,' expt)then varlist:=union({ cadr(m)},varlist) 205 else varlist:=union(indets(cdr(m)),varlist); 206 return varlist end; 207 208symbolic procedure buildassoc l; 209% Given a list of indeterminates(x1 x2 ...xn) we produce 210% an a-list of the form(( x1 . 1)(x2 . 2)...(xn . n)). 211begin integer i; 212 return(for each var in l collect progn(i:=i #+1,var . i)) end; 213 214symbolic procedure mons l; 215% Rewrite the leading monomials(i . e . GLTB). 216% the result is a list of monomials of the form : 217%(variable . exponent)or(( variable1 . exponent1)... 218% (variablen . exponentn)) 219% 220% mons('(LIST(EXPT Z 2)(EXPT X 2)(TIMES Y(EXPT X 3)))); 221%(((Y . 1)(X . 3))(X . 2)(Z . 2)). 222begin scalar monlist; 223 for each m in l do 224 if m neq ' list then monlist:= 225 if atom(m)then(m . 1). monlist 226 else if eqcar(m,' expt) 227 then(cadr m . caddr m). monlist 228 else(for each x in cdr(m)collect monsaux(x)) . monlist; 229 return monlist end; 230 231symbolic procedure monsaux m; 232 if eqcar(m,'expt)then cadr m . caddr m else m . 1; 233 234symbolic procedure lmon2arrmon m; 235% List-monomial to array-monomial 236% a list-monomial has the form:(variable_number . exponent) 237% or is a list with entries of this form. 238% "variable_number" is the number associated with the variable, 239% see buildassoc(). 240begin scalar mon;integer tdeg;mon:=allocmon nvars!*; 241 if listp m then 242 for each varnodotexp in m do 243 << setnthexp(mon,car varnodotexp,cdr varnodotexp); 244 tdeg:=tdeg + cdr varnodotexp >> 245 else 246 << setnthexp(mon,car m,cdr m);tdeg:=tdeg + cdr m >>; 247 settdeg(mon,tdeg);return mon end; 248 249symbolic procedure gltbfix l; 250% Sometimes GLTB has the form(list(list ...)) 251% instead of(list ...). 252 if listp cadr l and caadr(l)= ' list then cadr l else l; 253 254symbolic procedure gege(m1,m2); 255 if gettdeg(m1)>= gettdeg(m2)then t else nil; 256 257symbolic procedure getendptr l; 258begin scalar ptr;while l do << ptr:=l;l:=cdr l >>; 259 return ptr end; 260 261symbolic procedure gltb2arrideal xgltb; 262% Convert the monomial ideal given by GLTB(in list form) 263% to a list of vectors where each vector represents a monomial. 264begin scalar l;l:=indets(gltbfix(xgltb));nvars!*:=length(l); 265 l:=sublis(buildassoc(l), mons(gltbfix(xgltb))); 266 l:=for each m in l collect lmon2arrmon(m); 267 l:=sort(l,' gege); 268 return { getendptr(l), l } end; 269 270% ************** END OF INTERFACE TO ALGEBRAIC MODE ************** 271 272%************** PROCEDURES ************** 273 274symbolic procedure npol ideal; 275% Recursively computes the numerator of the Hilbert series. 276begin scalar v,si;v:=nextvar ideal; 277 if not v then return basecasepol ideal; 278 si:=splitideal(ideal,v); 279 return shiftadd(npol car si,npol cadr si)end; 280 281symbolic procedure dividesbyvar(var,mon); 282begin scalar div;if getnthexp(mon,var)= 0 then return nil; 283 div:=allocmon nvars!*; 284 for i:=1 : nvars!* do setnthexp(div,i,getnthexp(mon,i)); 285 setnthexp(div,var, getnthexp(mon,var)- 1); 286 settdeg(div,gettdeg mon - 1);return div end; 287 288symbolic procedure divides(m1,m2); 289% Does m1 divide m2? 290% m1 and m2 are monomials; 291% result: either nil(when m1 does not divide m2)or m2 / m1. 292begin scalar m,d,i;i:=1;m:=allocmon(nvars!*); 293 settdeg(m,d:=gettdeg(m2)- gettdeg(m1)); 294 while d >= 0 and i <= nvars!* do 295 << setnthexp(m,i,d:=getnthexp(m2,i)- getnthexp(m1,i)); 296 i:= i+1 >>; 297 return if d < 0 then nil else m end; 298 299symbolic procedure shiftadd(p1,p2); 300% p1 + z * p2; 301% p1 and p2 are polynomials(nonempty coefficient lists). 302begin scalar p,pptr;pptr:=p:=car p1 . nil; 303 p1:=cdr p1; 304 while p1 and p2 do 305 << rplacd(pptr,(car p1 + car p2). nil); 306 p1:=cdr p1;p2:=cdr p2;pptr:=cdr pptr >>; 307 if p1 then rplacd(pptr,p1) 308 else rplacd(pptr,p2);return p end; 309 310symbolic procedure remmult(ipp1,ipp2); 311% The union of two ideals with redundancy of generators eliminated. 312begin scalar fmon,inew,isearch,primeflag,x; 313 % fix;x is used in the macro... 314 x:=nil;inew:=theemptyideal(); 315 while notemptyideal(ipp1)and notemptyideal(ipp2)do 316 begin if gettdeg(firstmon(ipp2)) < gettdeg(firstmon(ipp1)) 317 then << fmon:=getnextmon(ipp1);isearch:=ipp2 >> 318 else << fmon:=getnextmon(ipp2);isearch:=ipp1 >>; 319 primeflag:=t; 320 while primeflag and notemptyideal(isearch)do 321 if divides(getnextmon(isearch), fmon)then primeflag:=nil; 322 if primeflag then addtoideal(fmon,inew)end; 323 return if notemptyideal(ipp1)then appendideals(inew,ipp1) 324 else appendideals(inew,ipp2)end; 325 326symbolic procedure nextvar ideal; 327% Extracts a variable in the ideal suitable for division. 328begin scalar m,var,x;x:=nil; 329 repeat 330 << m:=getnextmon ideal; 331 var:=getvarifnotsingle m; 332 >> until var or ideal = theemptyideal(); 333 return var end; 334 335symbolic procedure getvarifnotsingle mon; 336% Returns nil if the monomial is in a single variable, 337% otherwise the index of the second variable of the monomial. 338begin scalar foundvarflag,exp;integer i; 339 while not foundvarflag do 340 << i:=i + 1;exp:=getnthexp(mon,i); 341 if exp > 0 then foundvarflag:=t >>; 342 foundvarflag:=nil; 343 while i < nvars!* and not foundvarflag do 344 << i:=i + 1;exp:=getnthexp(mon,i); 345 if exp > 0 then foundvarflag:=t >>; 346 if foundvarflag then return i else return nil end; 347 348symbolic procedure makeonevarmon vindex; 349% Returns the monomial consisting of the single variable vindex. 350begin scalar mon;mon:=allocmon nvars!*; 351 for i:=1 : nvars!* do setnthexp(mon,i,0); 352 setnthexp(mon,vindex,1); 353 settdeg(mon,1);return mon end; 354 355symbolic procedure splitideal(ideal,var); 356% Splits the ideal into two simpler ideals. 357begin scalar div,ideal1,ideal2,m,x;x:=nil; 358 ideal1:=theemptyideal();ideal2:=theemptyideal(); 359 while notemptyideal(ideal)do 360 << m:=getnextmon(ideal); 361 if div:=dividesbyvar(var,m)then addtoideal(div,ideal2) 362 else addtoideal(m,ideal1)>>; 363 ideal2:=remmult(ideal1,ideal2);ideal1:=insertvar(var,ideal1); 364 return { ideal1,ideal2 } end; 365 366symbolic procedure basecasepol ideal; 367% In the base case every monomial is of the form Xi ^ ei; 368% result : the numerator polynomial of the Hilbert series 369% i.e.(1 - z ^ e1)*(1 - z ^ e2)* ... 370begin scalar p,degsofar,e;integer tdeg; 371 for each mon in cadr ideal do tdeg:=tdeg + gettdeg mon; 372 p:=makevector(tdeg,0);putv(p,0,1);degsofar:=0; 373 for each mon in cadr ideal do 374 << e:=gettdeg mon; 375 for j:= degsofar step -1 until 0 do 376 putv(p,j + e,getv(p,j+e)- getv(p,j)); 377 degsofar:=degsofar + e >>; 378 return vector2list p end; 379 380symbolic procedure vector2list v; 381% Convert a vector v to a list. No type checking is done. 382begin scalar u; 383 for i:=upbv v step -1 until 0 do u:=getv(v,i).u; 384 return u end; 385 386endmodule;;end; 387