document{quote programs,
"The file `programs' contains the following routines:
   Ainv, HilbCoeff, Hilbinfo, lastDegree, Polyfit, redno, ReesAlg, Reg,
   Reg1, RRclosure." 
}
document{quote Ainv,
"Ainv(R)
    --gives the last non-zero degree of the zero-th local cohomology of the
    --ring R with respect to the maximal ideal of R.  Gives the value of -1
    --if R has positive depth."
}
Ainv(Ring):= (R) -> (
I:=ideal(0_R);
M:=ideal(vars R);
H:=I:M;
lastDegree H
);
document{quote AssGr,
"AssGr(Ideal)-- computes the associated graded ring of the ideal."
}
AssGr(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x=quote x;
u=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
--Next, we get the ideal I in terms of the variables x_1,...,x_m:
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
--Now we create the ideal which will define the Rees algebra of I:
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F=ideal(J6);
S:=LL[x_1..x_m,u_1..u_n, Degrees=>drop(L3,1)];
F=substitute(F,S)+substitute(J,S);
G=S/F;
G 
);
document{quote depth,
"depth(Ring)-- computes the depth of a ring"
}
depth(Ring):= (X) -> (
I:=mingens(ideal(vars X));
m:= numgens(ideal(I));
d:=-1;
if (m>0) then (
L:=join(apply(i=0..m-1,i->(degree I_i)));
L=sort L;
r:=0;
B:=ideal(0_X);
E:=ideal(0_X);
U:=ideal(0_X);
while (B==E) do (
B=ideal(0_X);
r=first L;
E=ideal(1_X);
while (dim module E>0) do (
if not(hilbertFunction(r,X)==0) then (U=ideal(random(r,X));E=(B:U));
r=r+1);
d=d+1;
X=X/U)) else (d=d+1);
d
);
document{quote HilbCoeff,
"HilbCoeff(Function)
            -- given a polynomial function p(n) of degree d, this routine 
            -- returns a List consisting of the Hilbert Coefficients 
            -- {e_0,...,e_d}, where  
            -- p(n)= e_0*binomial(n+d-1,d)-e_1*binomial(n+d-2,d-1)+..."
               }
HilbCoeff(Function):= (p) -> (
x:=quote x;
R:=QQ[x];
g:=p x;
L:=degree g;
d:=L#0;
i:=d;
J:={};
while i>-1 do (
e:=substitute(g,x=>0_QQ);
e1:=(-1)^i*e;
J=join({e1},elements(J));
g=g-substitute(g,x=>x-1);
i=i-1;);
J 
);
document{quote Hilbinfo,
"Hilbinfo(I,n)
             -- this command returns a List consisting of the Hilbert
             -- polynomial of the ideal I, the postulation number of I, 
             -- and a List {e_0,..,e_d} of the Hilbert coefficients
             -- of I.  This routine uses n as an upper bound on the
             -- postulation number.

 Hilbinfo(I) -- computes Hilbinfo(I,n) where n=Reg I.  Since Reg I is
             -- an upper bound on the postulation number, this is 
             -- guaranteed to give you the correct Hilbert polynomial.
             -- However, Reg I may take a very long time for ideals
             -- with generators in several different degrees."}
Hilbinfo=method();
Hilbinfo(Ideal,ZZ):= (I,r) -> (
if dim I>0 then error "ideal is not zero-dimensional";
S:=ring I;
d:=dim S;
f:=n->degree I^n;
L:=elements(apply(i=1..d+1,i->(r+i,f(r+i))));
p:=Polyfit L;
v:=r+1;
c:=(p(v)-f(v))^2;
h:=f(v);
if c>0 then error "oops, there's a mistake";
while (c==0) do (
v=v-1;
if (v<0) then (h=0) else (h=f(v));
c=p(v)-h;);
C:=HilbCoeff p;
HP=p;
{HP,v,C}
);
Hilbinfo(Ideal):=(I) -> (
r:=Reg I;
Hilbinfo(I,r)
);
document{quote lastDegree,
"lastDegree(M)
          --for a non-negatively graded module, ring or ideal M of 
          --finite length, this gives the largest integer n such that
          --M_n is non-zero.  The value -1 is returned if the module
          --is zero." 
}
lastDegree= method();
lastDegree(Module):= (M) -> (
-- We first make sure that M has finite length:
if dim M>0 then error "module is not of finite length";
l:=degree M;
n:=0;
while l>0 do (
l=l-hilbertFunction(n,M);
n=n+1;
      );
n-1
);
lastDegree(Ideal):= (I) -> (
M:=module I;
lastDegree M
);
lastDegree(Ring):= (R) -> (
I:=ideal(vars R);
if (I==ideal(0_R)) then (0) else (lastDegree I)
);
document{quote Polyfit, 
"Polyfit(List)
            -- given a list of d points in Q x Q, Polyfit returns the
            -- polynomial function of degree d-1 which passes through
            -- those points."
 }
Polyfit(List):= (L) -> (
d:=#L;
j:=0;
i:=0;
x:=quote x;
X:=QQ[x];
g:=0_X;
h:=1_X;
while j<d do (
i=0;
h=1_X;
while i<d do (
if not(i==j) then (
h=(1_QQ/(L#j#0-L#i#0))*h*(x-L#i#0););
i=i+1;);
g=g+h*L#j#1;
j=j+1;);
f:=m->substitute(g,x=>m*1_QQ);
f
);
document{quote redno,
"redno(Ideal)
             -- this routine gives the reduction number of a generic
             -- reduction of an ideal"
 }
redno(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x:=quote x;
u:=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F:=ideal(J6);
S:=LL[u_1..u_n];
F=substitute(F,S);
F=ideal(mingens F);
T:=S/F;
U:=ideal(0_T);
E:=ideal(0_T);
while (dim T)>0 do (
U=ideal(random(1,T));
E=(ideal(0_T):U);
while (dim module E>0) do (U=ideal(random(1,T));E=(ideal(0_T):U););
T=T/U;);
r:=0;
A:=ideal(vars T);
while not(A==ideal(0_T)) do (r=r+1;A=ideal(vars T)*A;);
r
);
document{quote ReesAlg,
"ReesAlg(Ideal)-- produces the Rees Ring R[It] of the ideal I"
}
ReesAlg(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x=quote x;
u=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
--Next, we get the ideal I in terms of the variables x_1,...,x_m:
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
--Now we create the ideal which will define the Rees algebra of I:
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F=ideal(J6);
S:=LL[x_1..x_m,u_1..u_n, Degrees=>drop(L3,1)];
F=substitute(F,S);
T=S/F;
T 
);
document{quote Reg,
" Reg(Ring)-- returns the Castelnuovo-Mumford regularity of R.
 Reg(Ideal)-- gives Reg(G) where G is the associated graded ring (or the
              Rees ring) of the ideal.  Note:  the ideal must be zero- 
              dimensional.  Also, Reg(I) may run very slowly if I has 
              several minimal generators of different degrees."
}
Reg= method();
Reg(Ring):= (R) -> (
a:= Ainv R;
while dim R>0 do (
d:=1;
while d>0 do (
I:=ideal(random(1,R));
F:=(ideal(0_R):I);
d=dim module F;);
R=R/I;
b:=Ainv R;
if b>a then a=b;);
a
);
Reg(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x:=quote x;
u=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
--Next, we get the ideal I in terms of the variables x_1,...,x_m:
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
--Now we create the ideal which will define the Rees algebra of I:
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F:=ideal(J6);
S:=LL[x_1..x_m,u_1..u_n, Degrees=>drop(L3,1)];
F=substitute(F,S)+substitute(J,S);
T:=S/F;
A:=ideal(apply(i=1..m,i->(x_i)*(1_T)));
T1:=ideal(apply(i=1..n,i->(u_i)*(1_T)));
a:=0;
R:=sort L1;
r:=0;
q:=0;
p:=0;
d:=0;
b:=0;
B:=ideal(0_T);
E:=ideal(0_T);
U:=ideal(0_T);
g:=quote g;
while (dim A)>0 do (
r=first R;
q=r;
p=1;
d=0;
R=drop(R,1);
c1=matrix{{apply(i=1..m,i->(x_i)*(0_T))}};
c2=matrix{{apply(i=1..n,i->(u_i)*(1_T))}};
c3=matrix{{c1,c2}};
g=map(T,T,c3);
U=g(ideal(random(r,T)));
E=(ideal(0_T):U);
while (dim module E>0) do (p=p+1;r=p*q;U=g(ideal(random(r,T)));
E=(ideal(0_T):U));
T1=ideal(apply(i=1..n,i->(u_i)*(1_T)));
B=ideal(0_T);
b=a+d;
E=intersect(E,T1^b);
while not(E==B) do (
b=b+1;
E=intersect(E,T1^(b+1)));
T=T/U;
A=ideal(apply(i=1..m,i->(x_i)*(1_T)));
if (b-d)>a then a=b-d;
d=d+p-1);
T1=ideal(apply(i=1..n,i->(u_i)*(1_T)));
b=a+d;
B=ideal(0_T);
while not(T1^(b+1)==B) do (
b=b+1;);
if b>a+d then a=b-d;
a
);
document{quote Reg1,
"Reg1(Ideal)-- gives Reg(S) where S=R[It]/mR[It]"
}
Reg1(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x:=quote x;
u:=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F:=ideal(J6);
S:=LL[u_1..u_n];
F=substitute(F,S);
F=ideal(mingens F);
T:=S/F;
Reg(T)
);
document{quote RRclosure,
"RRclosure(I)-- yields the Ratliff-Rush closure of I."
}
RRclosure(Ideal):= (I) -> (
X:= ring I;
m:= numgens X;
LL:= coefficientRing X;
-- We next find the number of generators of I:
I=ideal(mingens I);
n:=numgens I;
W:=gens I;
L1:=join(apply(i=0..n-1,i->(degree W_i)+{1}));
L2:=join(apply(i=0..m-1,i->(degree (vars X)_i)));
L3:=join({1},L2,L1);
-- Now we create a polynomial ring over LL in m+n+1 variables:
t:=quote t;
x:=quote x;
u:=quote u;
R1:=LL[t,x_1..x_m,u_1..u_n, MonomialOrder=> Eliminate 1,Degrees=>L3];
--Next, we get the ideal I in terms of the variables x_1,...,x_m:
f:=map(R1, X, elements(x_1..x_m));
J:=f I;
--Now we create the ideal which will define the Rees algebra of I:
J2:=gens J;
J3:=matrix{elements(u_1..u_n)};
J4:=J3-t*J2;
J5:=ideal(J4);
J6:=selectInSubring(1, gens gb J5);
F:=J+ideal(J6);
S:=LL[x_1..x_m,u_1..u_n, Degrees=>drop(L3,1)];
F=substitute(F,S);
H:=saturate(F);
--Now map back to the original ring, sending the extra variables to 0:
c1:=matrix{apply(elements(1..n), i->0_X)};
c2:=matrix{{vars X, c1}};
g:=map(X,S,c2);
B:=g H;
--The Ratliff-Rush closure of I is:
ideal(mingens B)
);