--* From BMT%WATSON.vnet.ibm.com@yktvmv.watson.ibm.com Wed Apr 13 06:59:19 1994 --* Received: from yktvmv.watson.ibm.com by leonardo.watson.ibm.com (AIX 3.2/UCB 5.64/920123) --* id AA26628; Wed, 13 Apr 1994 06:59:19 -0400 --* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3) --* with BSMTP id 7613; Wed, 13 Apr 94 06:58:53 EDT --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id --* for asbugs@watson; Wed, 13 Apr 94 06:58:53 -0400 --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id 4351; Wed, 13 Apr 1994 06:58:52 EDT --* Received: from spadserv.watson.ibm.com by yktvmv.watson.ibm.com --* (IBM VM SMTP V2R3) with TCP; Wed, 13 Apr 94 06:58:52 EDT --* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524) --* id AA14733; Wed, 13 Apr 1994 07:00:12 -0400 --* Date: Wed, 13 Apr 1994 07:00:12 -0400 --* From: bmt@spadserv.watson.ibm.com --* X-External-Networks: yes --* Message-Id: <9404131100.AA14733@spadserv.watson.ibm.com> --* To: asbugs@watson.ibm.com --* Subject: [1] bogus complaint about no-meaning for 'm' already used in previous line, then compiler dies. [hilbert4.as][34.6] --@ Fixed by: SSD Thu Jun 23 22:58:14 EDT 1994 --@ Tested by: none --@ Summary: Bogus message no longer appears. --> testcomp -O -Q inline-all -l asdem --> testrun -O -Q inline-all -l asdem #include "axiom.as" -- This file computes hilbert functions for monomial ideals -- ref: "On the Computation of Hilbert-Poincare Series", Bigatti, Caboara, Robbiano, -- AAECC vol 2 #1 (1991) pp 21-33 macro Monom == Monomial L == List SI == SingleInteger B == Boolean POLY == SparseUnivariatePolynomial Integer Array == Vector export new: (SI, SI) -> Array SI export generator: Segment SI -> Generator SI export ..: (SI, SI) -> Segment SI export ..: SI -> Segment SI import from SingleInteger Monomial : OrderedSet with totalDegree: % -> SI divides?: (%, %) -> B homogLess: (%, %) -> B quo: (%, %) -> % quo: (%, SI) -> % *: (%, %) -> % pure?: % -> B varMonom: (i:SI,n:SI, deg:SI) -> % variables: L % -> L SI apply: (%, SI) -> SI #: % -> SI monomial: Tuple SI -> % == Array(SingleInteger) add Rep ==> Array(SingleInteger) import from Rep monomial(t:Tuple SI):% == per [ t ] totalDegree(m:%):SI == sum:SI := 0 for e in rep m repeat sum := sum + e sum divides?(m1:%, m2:%):B == for e1 in rep m1 for e2 in rep m2 repeat if e1 > e2 then return false true -- (m1:%) < (m2:%):B == -- for e1 in rep m1 for e2 in rep m2 repeat -- if e1 < e2 then return true -- if e1 > e2 then return false -- false -- (m1:%) > (m2:%):B == m2 < m1 homogLess(m1:%, m2:%):B == (d1:=totalDegree(m1)) < (d2:=totalDegree(m2)) => true d1 > d2 => false -- rep(m1) < rep(m2) true (m:%) quo (v:SI):% == --remove vth variable -- per [((if i=v then 0 else (rep m).i) for i in 1..#rep m)] m2:= copy rep m m2.v := 0 per m2 (m1:%) quo (m2:%):% == per [(max(a1-a2,0) for a1 in rep m1 for a2 in rep m2)] (m1:%) * (m2:%):% == per [(a1+a2 for a1 in rep m1 for a2 in rep m2)] pure?(m:%):B == varFlag := false for e in rep m repeat e ~= 0 => varFlag => return false varFlag := true true varMonom(i:SI,n:SI, deg:SI):% == -- per [((if j=i then deg else 0$SI) for j in 1..n)] m:Rep := new(n,0) m.i:=deg per m variables(I:L %) :L SI == #I < 1 => nil n:SI:=# rep first I ans : L SI := nil for v in 1..n repeat for m in I repeat (rep m).v ~= 0 => ans := cons(v, ans) break ans HilbertFunctionPackage: with Hilbert: L Monom -> POLY adjoin: (Monom, L Monom) -> L Monom == add adjoin(m:Monom, lm:L Monom):L Monom == empty?(lm) => cons(m, nil) ris1:L Monom:= empty() ris2:L Monom:= empty() while not empty? lm repeat m1 := first lm lm := rest lm m <= m1 => if not divides?(m,m1) then (ris1 := cons(m1, ris1)) iterate ris2 := cons(m1, ris2) if divides?(m1, m) then return concat!(reverse!(ris1), concat!(reverse! ris2, lm)) concat!(reverse!(ris1), cons(m, reverse! ris2)) reduce(lm:L Monom):L Monom == lm := sortHomogRemDup(lm) empty? lm => lm ris :L Monom := nil risd:L Monom := list first lm d := totalDegree first lm for m in rest lm repeat if totalDegree(m)=d then risd := cons(m, risd) else ris := mergeDiv(ris, risd) d := totalDegree m risd :=list m mergeDiv(ris, risd) merge(l1:L Monom, l2:L Monom):L Monom == #l1 > #l2 => merge(l2,l1) ris := l2 for m1 in l1 repeat ris := adjoin(m1, ris) ris mergeDiv( small:L Monom, big:L Monom): L Monom == ans : L Monom := small for m in big repeat if not contained?(m,small) then ans := cons(m, ans) ans contained?(m:Monom, id: L Monom) : B == for mm in id repeat divides?(mm, m) => return true false (I:L Monom) quo (m:Monom):L Monom == reduce [mm quo m for mm in I] pureSplit(l:List Monom):Record(pure:List Monom, mixed:List Monom) == pures := nil mixeds := nil for m in l repeat if pure? m then pures:=cons(m,pures) else mixeds := cons(m,mixeds) [pures, mixeds] sort(I:L Monom, v:SI):L Monom == sort((a:Monom,b:Monom):B+->(a.v < b.v), I)$ListSort(Monom) sortHomogRemDup(l:L Monom):L Monom == l:=sort(homogLess, l)$ListSort(Monom) empty? l => l ans:L Monom := list first l for m in rest l repeat if m ~= first(ans) then ans:=cons(m, ans) reverse! ans decompose(I:L Monom, v:SI):Record(size:SI, ideals:L L Monom, degs:L SI) == I := sort(I, v) idlist: L L Monom := nil deglist : L SI := nil size : SI := 0 J: L Monom := nil while not empty? I repeat d := first(I).v tj : L Monom := nil while not empty? I and d=(m:=first I).v repeat tj := cons(m quo v, tj) I := rest I J := mergeDiv(tj, J) idlist := cons(J, idlist) deglist := cons(d, deglist) size := size + #J [size, idlist, deglist] var(n:SI) : POLY == monomial(1$Integer, n::Integer) Hilbert(I:L Monom):POLY == #I = 0 => 1 -- no non-zero generators = 0 ideal #I = 1 => var(0) - var(totalDegree first I) lvar :L SI := variables I import from Record(size:SI, ideals:L L Monom, degs:L SI) Jbest := decompose(I, first lvar) for v in rest lvar while #I < Jbest.size repeat JJ := decompose(I, v) JJ.size < Jbest.size => Jbest := JJ import from L L Monom import from L SI Jold := first(Jbest.ideals) dold := first(Jbest.degs) f:=var(dold)*Hilbert(Jold) for J in rest Jbest.ideals for d in rest Jbest.degs repeat f := f + (var(d) - var(dold)) * Hilbert(J) dold := d var(0) - var(dold) + f MonomialIdealPackage: with varMonomsPower: (SI, SI) -> L Monom *: (L Monom, L Monom) -> L Monom ^: (L Monom, SI) -> L Monom == add varMonoms(n:SI):L Monom == [varMonom(i,n,1) for i in 1..n] varMonomsPower(n:SI, deg:SI):L Monom == n = 1 => [ monomial(deg)] ans : L Monom := nil for j in 0..deg repeat ans := concat(varMonomMult(j,varMonomsPower(n-1,deg-j)), ans) ans varMonomMult(i:SI, mlist : L Monom) : L Monom == [varMonomMult(i, m) for m in mlist] varMonomMult(i:SI, m:Monom) : Monom == nm:Rep := new(#m + 1, i) for k in 1..#m repeat nm.k :=m.k nm pretend Monom (i1:L Monom) * (i2:L Monom):L Monom == ans : L Monom := nil for m1 in i1 repeat for m2 in i2 repeat ans := adjoin(m1*m2, ans) ans (i:L Monom) ^ (n:SI) : L Monom == n = 1 => i odd? n => i * (i*i)^shift(n, -1) (i*i)^shift(n,-1)