open-axiom-test 1.4.1+svn~2626-2ubuntu2 / usr / lib / open-axiom / input / hilbert.as

-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are
-- met:
--
--     - Redistributions of source code must retain the above copyright
--       notice, this list of conditions and the following disclaimer.
--
--     - Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in
--       the documentation and/or other materials provided with the
--       distribution.
--
--     - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--       names of its contributors may be used to endorse or promote products
--       derived from this software without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
-- IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-- TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
-- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
-- OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


)compile hilbert.as
mon1 := monom(4,0,0,0)
mon2:= monom(3,3,0,0)
mon3 := monom(3,2,1,0)
mon4 := monom(3,1,2,0)
mon5 := monom(0,2,0,1)
mon6 := monom(0,1,0,5)
l := [mon1, mon2, mon3, mon4, mon5, mon6]
Hilbert l
idA := varMonomsPower(6,5);
#idA
Hilbert idA
idB := varMonomsPower(6,6);
#idB
Hilbert idB
idC := varMonomsPower(12,3);
#idC
Hilbert idC
idD:=[monom(2,0,0,0),monom(1,1,0,0),monom(1,0,1,0),monom(1,0,0,1),_
 monom(0,3,0,0),monom(0,2,1,0)]^4;
#idD
Hilbert idD

#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
import from Segment SI
import from Integer

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

      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)