--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Tue Mar 21 12:22:18 1995
--* Received: from yktvmv-ob.watson.ibm.com by watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA19269; Tue, 21 Mar 1995 12:22:18 -0500
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 7513; Tue, 21 Mar 95 12:21:19 EST
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.DOOLEY.NOTE.VAGENT2.0381.Mar.21.12:21:18.-0500>
--*           for asbugs@watson; Tue, 21 Mar 95 12:21:19 -0500
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 0375; Tue, 21 Mar 1995 12:21:18 EST
--* Received: from watson.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with TCP; Tue, 21 Mar 95 12:21:16 EST
--* Received: by watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA35136; Tue, 21 Mar 1995 12:21:43 -0500
--* Date: Tue, 21 Mar 1995 12:21:43 -0500
--* From: dooley@watson.ibm.com (Sam Dooley)
--* Message-Id: <9503211721.AA35136@watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: [2] Runtime bug in interoperable closure initialization.

--@ Fixed  by: <Who> <Date>
--@ Tested by: <Name of new or existing file in test directory>
--@ Summary:   <Description of real problem and the fix>

-- Command line: asharp -O -Fasy -Faso -Flsp -laxiom -Mno-AS_W_WillObsolete -DAxiom hilbert.as
-- Version: r1-0-6-2
-- Original bug file name: hilbert.as

#include "axiom.as"
#pile

-- 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

import from NonNegativeInteger
import from SingleInteger

Monomial : OrderedSet with
      totalDegree: % -> SI
      divides?: (%, %) -> B
      homogLess: (%, %) -> B
      quo: (%, %) -> %
      quo: (%, SI) -> %
      *: (%, %) -> %
      varMonom: (i:SI,n:SI, deg:SI) -> %
      variables: L % -> L SI
      apply: (%, SI) -> SI
      #: % -> SI
      monom: Tuple SI -> %
   == Array(SingleInteger) add
      Rep ==> Array(SingleInteger)
      import from Rep

      monom(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
	( m1 < m2)

      (m:%) quo (v:SI):% == --remove vth variable
#if 0
           --!! This version works in r1-0-6-2.
           per [((if i=v then 0 else (rep m).i) for i:SI in 1..#rep m)]
#else
           --!! This version fails when quo is inlined into decompose.
           --!! The EnvInfo slot mutates the closure for copy from an A#
           --!! closure to an Axiom closure and then tries to force it.
           m2:= copy rep m
           m2.v := 0
           per m2
#endif

      (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)]

      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 ==
	empty? I => nil
	n:SI:=# rep first I
	ans : L SI := nil
        v:SI:=0
        while (v:=v+1)<=n repeat
--	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:Monom := first lm
	  lm := rest lm
	  if m <= m1 then
	     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 := [m]
	mergeDiv(ris, risd)

      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]

      sort(I:L Monom, v:SI):L Monom ==
	sort((a:Monom,b:Monom):B+->(a.v < b.v), I)

      sortHomogRemDup(l:L Monom):L Monom ==
	l:=sort(homogLess, l)
	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
          local m:Monom
	  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)::Integer::SI)
	[size, idlist, deglist]


      var(n:SI) : SparseUnivariatePolynomial Integer ==
          monomial(1$Integer, n::Integer::NonNegativeInteger)


      Hilbert(I:L Monom):POLY ==
	empty? I => 1 -- no non-zero generators = 0 ideal
	empty? rest I => 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)::Integer::SI < Jbest.size repeat
	   JJ := decompose(I, v)
	   JJ.size < Jbest.size => Jbest := JJ
	import from L L Monom
	import from L SI
	Jold:List Monom := first(Jbest.ideals)
	dold:SI := first(Jbest.degs)
	f:SparseUnivariatePolynomial Integer:=var(dold)*Hilbert(Jold)
	for J:List Monom in rest Jbest.ideals for d:SI 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]
         i:SI:=0
	 [varMonom(i,n,1) while {free i; (i:=i+1)<=n}]

      varMonomsPower(n:SI, deg:SI):L Monom ==
	 n = 1 => [ monom(deg)]
	 ans : L Monom := nil
--	 for j in 0..deg repeat
         j:SI:=-1
	 while (j:=j+1)<=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:Array SI := new(#m + 1, i)
--	for k in 1..#m repeat nm.k :=m.k
        k:SI:=0
	while (k:=k+1)<=#m repeat nm.k :=m.k
	nm pretend Monom

      (i1:L Monom) * (i2:L Monom):L Monom ==
          import from HilbertFunctionPackage
	  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)



