-- May 15 2010:

-- Changes compared to version 5: handles degree alpha+1 better. Before the final output assumed gens were known in that degree
-- only when h_Z(alpha)=1 although the verbose version gave more info for intermediate steps using q and l.
-- Now it uses q and l both for the final output and for the intermediate steps.

-- Run as follows: findGrBettisVerbose({m1,...,mn}), where mi are the multiplicities
-- of a fat point subscheme. 

findGrBettisVerbose = (lll) -> (l:=prmt(lll);
<< endl << "Z = " << l << endl;
<< "The number of points here is: " << #l << endl;
<< "(When the number of points is more than 8 the result may be only conjectural.)" << endl << endl;
i:=0;
j:=0;
tmp1:=0;
tmp2:=0;
gens:=0;
nef:={};
ql:={};
ll:={};
syz:=0;
v:={};
aa:=findalpha(l);
sig:=findtau(l)+1;
<< "alpha = " << aa << endl << endl;
<< "tau = " << sig-1 << endl << endl;
v={aa-1,0,0,0};
MM:=matrix({v});
MMM:=matrix({v});
ddim:=0;
ddim0:=0;
ddim1:=0;
ddim2:=0;
ddim3:=0;

for tt from aa to sig+2 do (
<< "Zariski decomposition in degree " << tt << ":" << endl;
decomp({tt,l});
nef=nefpart({tt,l});
nef={nef#0+1,nef#1};
ql={};
ll={};
scan(#(nef#1), i->(
        if i==0 then ql=join(ql,{1+(nef#1)#i});
        if i > 0 then ql=join(ql,{(nef#1)#i})));
scan(#(nef#1), i->(
        if i==0 then ll=join(ll,{-1+(nef#1)#i});
        if i > 0 then ll=join(ll,{(nef#1)#i})));
-- << ql << " " << ll << endl;
ddim3=homcompdim(fundom({tt,l}));
tmp2=homcompdim(fundom({nef#0-1,ql}));
tmp1=homcompdim(fundom({nef#0-2,ll}));
j=homcompdim(nef);
<< "H-E_1 is Cremona equivalent to: " << fundom({nef#0-1,ql}) << endl;
<< "H-(L-E_1) is Cremona equivalent to: " << fundom({nef#0-2,ll}) << endl;
<< "h^0(F) = h^0(H) = " << i << endl;
if tt == aa then v={tt,ddim3,ddim3,0};
if tt == aa then MM=MM||matrix({v});
if tt == aa then MMM=MMM||matrix({v});
tmp2=tmp1+tmp2;
if 3*ddim3-j>tmp1 then tmp1=3*ddim3-j;
<< tmp1 << " <= dim ker mu_H <= " << tmp2 << " by q + l theorem (where q = h^0(H-E_1) and l = h^0(H-(L-E_1)))" << endl;
<< "h^0(H + L) = " << j << endl;
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=homcompdim(fundom({tt+1,l}));
<< "h^0(F + L) = " << ddim3 << endl;
tmp1=j-(3*ddim2-tmp1);
tmp2=j-(3*ddim2-tmp2);
if tmp1 < 0 then tmp1 = 0;
if tmp2 < 0 then tmp2 = 0;
<< tmp1 << " <= dim cok mu_H <= " << tmp2 << endl << endl << endl;
gens=tmp2+ddim3-j;
syz=gens-(ddim3-3*ddim2+3*ddim1-ddim0);
-- if ddim2 == 1 then v={tt+1, ddim3,ddim3-3,0};
-- if ddim2>1 and tmp1 == tmp2 then v={tt+1, ddim3,gens,syz};
-- if ddim2>1 and tmp1 < tmp2 then v={tt+1, ddim3,-1,-1};
if tmp1 == tmp2 then v={tt+1, ddim3,gens,syz};
if tmp1 < tmp2 then v={tt+1, ddim3,-1,-1};
MM=MM||matrix({v});
if tt == aa then MMM=MMM||matrix({v}));
   << "Graded Betti numbers for Z = " << l << endl << endl;
   << "Column 1 = deg t" << endl;
   << "Column 2 = hilbert function in degree t" << endl;
   << "Column 3 = # generators in degree t" << endl;
   << "Column 4 = # syzygies in degree t" << endl << endl;
if #l <= 8 then findresReticent(l) else (
   << "Note: an entry of -1 means that entry cannot be computed by this script in the usual way," << endl;
   << "but bounds (which may be exact) in that case are given above" << endl;
   << "(see bounds on dim cok mu_H plus (h^0(F + L)-h^0(H + L)) in each degree)." << endl;
   << "The only conjecture used here is the SHGH Conjecture giving the Hilbert " << endl;
   << "function of fat point schemes. No Grobner basis computation was used." << endl << endl;
<< MM << endl;
<< endl;
   << "We now redo the computation, assuming conjectures about nonreduced structures on P1 and" << endl;
   << "using Grobner bases to compute splittings of bundles on P1. This allows us to give" << endl;
   << "expected values for the graded Betti numbers in all degrees except possibly alpha+1." << endl << endl;
ddim=homcompdim(fundom({aa,l}));
ddim0=0;
ddim1=0;
ddim2=0;
ddim3=ddim;
tt=aa+1;
gens=ddim;
ddim=homcompdim(fundom({tt,l}));
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=ddim;
while (tt <= sig+1) do (
tt=tt+1;
ddim=homcompdim(fundom({tt,l}));
nef=nefpart({tt-1,l});
nef={nef#0+1,nef#1};
gens=decompDeluxe({tt-2,l},ddim,homcompdim(nef));
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=ddim;
syz=gens-(ddim3-3*ddim2+3*ddim1-ddim0);
v={tt,ddim,gens,syz};
MMM=MMM||matrix({v});
);
<< MMM << endl))


-- Run as follows: findGrBettis({m1,...,mn}), where mi are the multiplicities
-- of a fat point subscheme. 

findGrBettis = (lll) -> (l:=prmt(lll);
<< "Graded Betti numbers for Z = " << l << endl;
<< "The number of points here is: " << #l << endl;
<< "(When the number of points is more than 8 the result may be only conjectural.)" << endl << endl;
i:=0;
j:=0;
tmp1:=0;
tmp2:=0;
gens:=0;
nef:={};
ql:={};
ll:={};
syz:=0;
v:={};
aa:=findalpha(l);
sig:=findtau(l)+1;
v={aa-1,0,0,0};
MM:=matrix({v});
MMM:=matrix({v});
ddim:=0;
ddim0:=0;
ddim1:=0;
ddim2:=0;
ddim3:=0;

for tt from aa to sig+2 do (
decompnonverbose({tt,l});
nef=nefpart({tt,l});
nef={nef#0+1,nef#1};
ql={};
ll={};
scan(#(nef#1), i->(
        if i==0 then ql=join(ql,{1+(nef#1)#i});
        if i > 0 then ql=join(ql,{(nef#1)#i})));
scan(#(nef#1), i->(
        if i==0 then ll=join(ll,{-1+(nef#1)#i});
        if i > 0 then ll=join(ll,{(nef#1)#i})));
ddim3=homcompdim(fundom({tt,l}));
tmp2=homcompdim(fundom({nef#0-1,ql}));
tmp1=homcompdim(fundom({nef#0-2,ll}));
j=homcompdim(nef);
if tt == aa then v={tt,ddim3,ddim3,0};
if tt == aa then MM=MM||matrix({v});
if tt == aa then MMM=MMM||matrix({v});
tmp2=tmp1+tmp2;
if 3*ddim3-j>tmp1 then tmp1=3*ddim3-j;
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=homcompdim(fundom({tt+1,l}));
tmp1=j-(3*ddim2-tmp1);
tmp2=j-(3*ddim2-tmp2);
if tmp1 < 0 then tmp1 = 0;
if tmp2 < 0 then tmp2 = 0;
gens=tmp2+ddim3-j;
syz=gens-(ddim3-3*ddim2+3*ddim1-ddim0);
if tmp1 == tmp2 then v={tt+1, ddim3,gens,syz};
if tmp1 < tmp2 then v={tt+1, ddim3,-1,-1};
MM=MM||matrix({v});
if tt == aa then MMM=MMM||matrix({v}));
   << "Column 1 = deg t" << endl;
   << "Column 2 = hilbert function in degree t" << endl;
   << "Column 3 = # generators in degree t" << endl;
   << "Column 4 = # syzygies in degree t" << endl << endl;
if #l <= 8 then findresReticent(l) else (
   << "Note: The only conjecture used here is the SHGH Conjecture giving the Hilbert " << endl;
   << "function of fat point schemes. No Grobner basis computation was used." << endl;
   << "An entry of -1 means that entry cannot be computed by the method used." << endl << endl;
<< MM << endl;
<< endl;
   << "We now redo the computation, assuming conjectures about nonreduced structures on P1 and" << endl;
   << "using Grobner bases to compute splittings of bundles on P1. This allows us to give" << endl;
   << "expected values for the graded Betti numbers in all degrees except possibly alpha+1." << endl << endl;
ddim=homcompdim(fundom({aa,l}));
ddim0=0;
ddim1=0;
ddim2=0;
ddim3=ddim;
tt=aa+1;
gens=ddim;
ddim=homcompdim(fundom({tt,l}));
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=ddim;
while (tt <= sig+1) do (
tt=tt+1;
ddim=homcompdim(fundom({tt,l}));
nef=nefpart({tt-1,l});
nef={nef#0+1,nef#1};
gens=decompDeluxenonverbose({tt-2,l},ddim,homcompdim(nef));
ddim0=ddim1;
ddim1=ddim2;
ddim2=ddim3;
ddim3=ddim;
syz=gens-(ddim3-3*ddim2+3*ddim1-ddim0);
v={tt,ddim,gens,syz};
MMM=MMM||matrix({v});
);
<< MMM << endl))


-- March 20, 2004
-- Run bundledecomp(l), where l is an exceptional class,
-- to get the bundle decomp of cotangent bundle of P2 pulled back
-- to l.

R=ZZ/31991[x,y, MonomialSize=>16, MonomialOrder => ProductOrder {1,1}];

S=ZZ/31991[a,b,c, MonomialSize=>16, MonomialOrder => ProductOrder {1,1,1}];

-- input: list of form l = {d,{m1,...,mn}} which reduces by 
-- quadratic transformations to {1,{1,1,0,...,0}}. Thus there is a
-- birational morphism f: P1 -> P2 parametrizing the curve {d,{m1,...,mn}}.
-- output: returns the result of restrictionIdeal(l,0) 
-- and prints the degrees d1 and d2 of the splitting of the bundle f*\Omega_P2.
-- Note that d1+d2=3d. The script uses the fact
-- that the splitting is known if d-m1 <= m1 + 1 to avoid doing any
-- computation in this case.

bundledecomp = (l) -> (
if #(l#1) < 3 then (
<< "number of points is " << #(l#1) << " but must be 3 or more; please try again" << endl) else
(
i:=0;
aE:=0;
bE:=0;
v:=l#1;
ev:={1,1};
scan(#v-2, i->(ev=join(ev,{0})));
ev={1,ev};  -- ev is the exceptional curve L-E_1-E_2
if ev != fundomB(l) then (
<< "Input class is not an exceptional class; ignore output" << endl) else (
d1:=l#0;
-- m1:=l#1#0;
m1:= (prmt(l#1))#0;
if d1-m1<m1+1 then (
aE=d1-m1;
bE=m1);
if d1-m1==m1+1 then (
aE=m1;
bE=d1-m1);
if d1-m1>m1+1 then (
<< "Using Grobner basis calculation to find splitting..." << endl << endl;
ii:=restrictionIdeal(l,0);
I:=(ii)_0;
J:=(ii)_1;
z:=res I;
aE=(degree (z_2)_0)_0-d1;
bE=(degree (z_2)_1)_0-d1;
));
{aE,bE}))

bundledecompnonverbose = (l) -> (
if #(l#1) < 3 then (
<< "number of points is " << #(l#1) << " but must be 3 or more; please try again" << endl) else
(
i:=0;
aE:=0;
bE:=0;
v:=l#1;
ev:={1,1};
scan(#v-2, i->(ev=join(ev,{0})));
ev={1,ev};  -- ev is the exceptional curve L-E_1-E_2
if ev != fundomB(l) then (
<< "Input class is not an exceptional class; ignore output" << endl) else (
d1:=l#0;
-- m1:=l#1#0;
m1:= (prmt(l#1))#0;
if d1-m1<m1+1 then (
aE=d1-m1;
bE=m1);
if d1-m1==m1+1 then (
aE=m1;
bE=d1-m1);
if d1-m1>m1+1 then (
ii:=restrictionIdeal(l,0);
I:=(ii)_0;
J:=(ii)_1;
z:=res I;
aE=(degree (z_2)_0)_0-d1;
bE=(degree (z_2)_1)_0-d1;
));
{aE,bE}))


-- lstprmt: arranges the elements of the list l1 in descending order,
-- and applies the same permutation to l2

lstprmt = (l1,l2) -> (  
tmpv1:=l1;
tmpv2:=l2;
v1:=l1;
v2:=l2;
i:=0;
j:=0;
k:=0;
scan(#l1, i->(scan(#l1, j->(
            if tmpv1#i < tmpv1#j then (
               if i < j then (
                  k=-1;
                  v1={};
                  v2={};
                  while(k<#l1-1) do (k=k+1;
                     if k==i then (v1=join(v1,{tmpv1#j});
                        v2=join(v2,{tmpv2#j})) else (
                        if k==j then (v1=join(v1,{tmpv1#i});
                           v2=join(v2,{tmpv2#i})) else (v1=join(v1,{tmpv1#k});
                              v2=join(v2,{tmpv2#k}))));
                  tmpv1=v1;
                  tmpv2=v2))))));
{v1,v2})


-- quad: performs a quadratic transform on a divisor class dE_0-(m_1E_1+...+m_nE_n).
-- Call it as quad({d,{m1,...,mn}}). The output is
-- {2d-m1-m2-m3,{d-m2-m3,d-m1-m3,d-m1-m2,m4,...,mn}}.

quad = (l) -> (
i:=0;
w:=l; 
if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
v:={w#0 - w#1#1 - w#1#2,w#0 - w#1#0 - w#1#2,w#0 - w#1#0 - w#1#1};
scan(#w#1, i->(if i>2 then v=join(v,{w#1#i})));
v={2*(w#0) - w#1#0 - w#1#1 - w#1#2,v};
v)


-- fundombothB: applies quad iteratively to l1 to reduce l1 to standard form
-- (this being l1 = {x,{y1,y2,...,yn}}, where x <= 1 and y1>=y2>= ...;
-- when l1 starts out being an exceptional class, it ends up as {1,{1,1,0,...,0}}).
-- domain of a certain group operation, and applies the same
-- group operation g to l2 (note: l1 and l2 should have the same number of entries). 
-- If l2 starts out in the fundamental domain,
-- and {l1',l2'}=fundombothB(l1,l2), then {l2,g^{-1}l}=fundombothB(l2',l).
-- This allows one to compute the action of g^{-1}.

fundombothB = (l1,l2) -> (
w1:=l1; 
w2:=l2;
v:={}; 
if #l1#1<3 then (
w1={l1#0,join(l1#1,{0,0,0})};
w2={l2#0,join(l2#1,{0,0,0})};
);
v=lstprmt(w1#1,w2#1);
w1={w1#0,v#0};
w2={w2#0,v#1};
while ((w1#0 < w1#1#0 + w1#1#1 + w1#1#2) and (w1#0 > 1)) do (
   w1=quad(w1);
   w2=quad(w2);
   v=lstprmt(w1#1,w2#1);
   w1={w1#0,v#0};
   w2={w2#0,v#1};
);
{w1,w2})

-- fundomB: Call it as fundomB({d,{m1,...,mn}}). The output is a new 
-- list {d',{m1',...,mn'}}; the class dE_0-(m_1E_1+...+m_nE_n) is
-- equivalent via Cremona transformations to d'E_0-(m_1'E_1+...+m_n'E_n),
-- where d' is the smallest positive integer possible.

fundomB = (l) -> (
(fundombothB(l,l))#0)

-- input: list {m1,...,mn} of integers and a 3 by n matrix M
-- output: list of integers in descending order, with matrix N whose columns 
-- are those of M permuted by the same permutation that arranged the integers

matprmt = (l,M) -> (
i:=0;
tmpprm:={};
scan(#l, i->(tmpprm=join(tmpprm,{i})));
myoutput:=lstprmt(l,tmpprm);
N:=M_(myoutput#1);
{myoutput#0,N})

-- input: l is an exceptional class, MM is a matrix of n=#l#1 points,
-- and flag is an integer; if flag=1, then extra steps are printed out.
-- Making MM an input allows use of a given set of points, if desired.
-- Output: the parametrization of the curve l

myloop = (l, MM, flag) -> (
mymaps:={};
tmpM:=id_(S^3);
npts:=#l#1;
M:=MM;
ev:={1,1};
i:=0;
scan(npts-2, i->(ev=join(ev,{0})));
ev={1,ev};  -- ev is the exceptional curve L-E_1-E_2
fe:={};
scan(npts, i->(fe=join(fe,{npts-i})));
fe={3*npts-2,fe}; -- element of fundamental domain
tfe:=(fundombothB(l,fe))_1; -- this element is a translate of previous fe;
                         -- translating fe back to fund dom takes ev to l
v:={};
f:={x,y,0*x};   -- this is the parametrization of ev, where ev is z=0 in P2
if flag==1 then << "test elt: " << tfe << " exc class: " << ev << endl;
if flag==1 then << "the points in P2: " << M << endl;
if flag==1 then << "Parametrization of exc class: " << toString f << endl;
line0:=det(submatrix(M,{1,2})|matrix{{a},{b},{c}});
line1:=det(submatrix(M,{0,2})|matrix{{a},{b},{c}});
line2:=det(submatrix(M,{0,1})|matrix{{a},{b},{c}});
g:=map(R,S,f);   -- this is the map S=k[P2]->R=k[P1] parametrizing ev
K:=kernel(g); -- compute equation of image in P2 of ev
I:=ideal(f#0,f#1,f#2);
J:=saturate I;
if flag==1 then << "Equation of image of exc class: " << toString K << endl;
while(tfe != fe) do (
   v=lstprmt(tfe#1,ev#1);
   M=(matprmt(tfe#1,M))#1;
   tfe={tfe#0,v#0};
   ev={ev#0,v#1};
   if (tfe#0 < tfe#1#0 + tfe#1#1 + tfe#1#2) then (
     line0=det(submatrix(M,{1,2})|matrix{{a},{b},{c}});
     -- << "line0: " << line0 << endl;
     line1=det(submatrix(M,{0,2})|matrix{{a},{b},{c}});
     -- << "line1: " << line1 << endl;
     line2=det(submatrix(M,{0,1})|matrix{{a},{b},{c}});
     -- << "line2: " << line2 << endl;
     -- << "param before sat: " << {g(line1*line2),g(line0*line2),g(line1*line0)} << endl;
     I={g(line1*line2),g(line0*line2),g(line1*line0)};
     J= oldergcd(oldergcd(I_0,I_1),I_2);  -- I typically has base points that need to be removed; J defines them
     -- << "J: " << J << endl;
     f={I_0//J,I_1//J,I_2//J}; -- this removes the base points
     -- << "param after sat: " << f << endl;
     g=map(R,S,f);  
     if flag==1 then K=kernel(g);  -- to see equations of the images of the curves, uncomment out: THIS IS SLOW
     tfe=quad(tfe);
     ev=quad(ev);
     -- << "After perm: " << M << endl;
     mymaps={}; -- for each point, define map given by evaluation at that point 
     scan(npts, i->(mymaps=join(mymaps,{map(S,S,{M_(0,i), M_(1,i), M_(2,i)})})));
     tmpM=id_(S^3);
     scan(npts, i->(if i>2 then tmpM=tmpM|matrix{{mymaps_i(line1*line2)},{mymaps_i(line0*line2)},{mymaps_i(line1*line0)}}));
     M=tmpM;
     -- uncomment the next line to see which points are on the given curve:
     -- scan(npts, i->(<< "pt " << i << " : " << mymaps_i(K_0) << endl ));  
     );
   if flag==1 then << "test elt: " << tfe << " exc class: " << ev << endl;
   if flag==1 then << "the points in P2: " << M << endl;
   if flag==1 then << "Parametrization of exc class: " << toString f << endl;
   if flag==1 then << "Equation of image of exc class: " << toString K << endl;  -- to see equations of the images of the curves, uncomment out
   );
f)

-- input: list of form l = {d,{m1,...,mn}} which reduces by 
-- quadratic transformations to {1,{1,1,0,...,0}}. Thus there is a
-- birational morphism f: P1 -> P2 parametrizing the curve {d,{m1,...,mn}}.
-- output: prints res of ideal on P1 of the linear forms in P2 pulled back via f
-- to the curve {d,{m1,...,mn}}. 
-- Note: flag is an integer; if flag=1, then extra steps are printed out.


restrictionIdeal = (l,flag) -> (
I:=ideal(0*x);
J:=ideal(0*x);
M:=matrix(table(3,#l#1-3,(i,j)->random(0,S)));
M=id_(S^3)|M;
f:=myloop(l,M,0);
I=ideal(f_0,f_1,f_2);
if flag==1 then << "Resolution of ideal of restrictions" << endl << betti res coker gens I << endl;
myM:=syz gens I;
J=ideal(myM_(0,0),myM_(1,0),myM_(2,0));
if flag==1 then << "Resolution of least degree relation ideal" << endl << betti res coker gens J << endl;
{I,J})


-- Run as GHI({6,4,5,6}) where the inputs are mults of points

GHI = (l)->(j:=0;
k:=0;
i:=0;
                   b:=ideal (matrix {{1}}**R);
                   scan(#l, i->(
                       f:=random(R^1,R^{-1});
                       g:=random(R^1,R^{-1});
                       I:=(ideal (f | g))^(l#i);
                       b=intersect(b,I)));
                   print betti res coker gens b;
i=0;
a:=hilbertFunction(i,b);
aa:=hilbertFunction(i+1,b);
while a < aa do (
j=hilbertFunction(i,R) - a;
if j>0 then (
if k==0 then << "hilb func starting in degree " << i << ": " << j << " ";
if k==1 then << j << " ";
k=1);
i=i+1;
a=aa;
aa=hilbertFunction(i+1,b);
);
if k==0 then << "hilb func starting in degree " << i << ": " << hilbertFunction(i,R) - a << " " << hilbertFunction(i+1,R) - aa << " " << hilbertFunction(i+2,R) - hilbertFunction(i+2,b) << " " << endl;
if k==1 then << hilbertFunction(i,R) - a << " " << hilbertFunction(i+1,R) - aa << " " << hilbertFunction(i+2,R) - hilbertFunction(i+2,b) << " " << endl;
b)



-- These routines have been debugged on MACAULAY 2, version 0.8.52
-- Brian Harbourne, October 12, 2000
-- (July 30, 2001: revised scripts unifbounds, ezunifHRalpha and ezunifHRalphaB)

-- findres: This computes the syzygy modules in any resolution
-- of the saturated homogeneous ideal defining any eight or fewer general 
-- fat points of P2. The hilbert function of the ideal is also found.
-- Call it as findres({m_1,...,m_n}) for n <= eight integers m_i.
-- Note that findres does not rely on Grobner bases, so it is fast by comparison.

findres = (l) -> (
if #l>8 then (
   << "This script works only for up to 8 points." << endl;
   << "Please try again with an input list of at most 8 integers." << endl)
else (
   i:=0;
   myflag2:=0;
   w2:={};
   dd1:=0;
   myker:=0;
   www:={};
   ww:=l;
-- the list l of multiplicities is, for simplicity, extended if need be
-- so that it has 8 elements.
   while(#ww < 8) do ww=join(ww,{0});
   n:=#ww;
   n=n-1;
   ww=zr(ww);  -- zero out negative elements of the list
   a1:=findalpha(ww);  -- find alpha, the least degree t such that I_t \ne 0
   d1:=a1-2;
   tau:=findtau(ww);
   v4:={}; -- list of number of syzygies in each degree t listed in v0
   v3:={}; -- list of dim of coker of \mu_t in each degree t listed in v0
   v2:={}; -- list of dim of ker of \mu_t in each degree t listed in v0
   v1:={}; -- list of Hilbert function values for each degree listed in v0
   v0:={}; -- list of degrees from alpha-2 to tau+2, where tau is the least 
         -- degree such that the fat points impose independent conditions
   while (d1 <= tau+2) do ( -- loop from alpha-2 to tau+2, computing v0, v1 and v2
-- append the current degree d1 to the list of degrees
      v0=join(v0,{d1});
-- append the value of the hilbert function in degree d1
-- to the list of values of the hilbert function
      v1=join(v1,{homcompdim(fundom({d1,ww}))});
-- now compute and append to the list v2 the dimension myker of the 
-- kernel of the map \mu_t : I_t\otimes k[P2]_1 \to I_{t+1} where t=d1 
-- and I is the ideal of the fat points subscheme
      if d1<a1 then (v2=join(v2,{0}));  -- d1<a1 means I_{d1}=0 so myker=0
      if d1>=a1 then (     -- for d1>=a1, compute myker 
         myflag2=0;  
         w2=ww;
         dd1=d1;  
         while(myflag2==0) do ( -- this loop implements the main theorem of [FHH]
                                -- which gives an algorithm for computing myker
            w2=zr(w2);  
            w2=prmt(w2);  
            if homcompdim({dd1,w2})==0 then myflag2=1 else (  
              if dd1*6 - (dot(w2,{3,2,2,2,2,2,2,2})) <= 2 then (  
                 dd1=dd1-6;  
                 w2 = {(w2#0)-3,(w2#1)-2,(w2#2)-2,(w2#3)-2,(w2#4)-2,
                 (w2#5)-2,(w2#6)-2,(w2#7)-2}) else (  
                 if dd1*5 - (dot(w2,{2,2,2,2,2,2,1,1})) <= 1 then (  
                    dd1=dd1-5;  
                    w2 = {(w2#0)-2,(w2#1)-2,(w2#2)-2,(w2#3)-2,(w2#4)-2,
                    (w2#5)-2,(w2#6)-1,(w2#7)-1}) else (  
                    if dd1*4 - (dot(w2,{2,2,2,1,1,1,1,1})) <= 1 then (  
                       dd1=dd1-4;  
                       w2 = {(w2#0)-2,(w2#1)-2,(w2#2)-2,(w2#3)-1,(w2#4)-1,
                       (w2#5)-1,(w2#6)-1,(w2#7)-1}) else (  
                       if dd1*3 - (dot(w2,{2,1,1,1,1,1,1,0})) <= 0 then (  
                          dd1=dd1-3;  
                          w2 = {(w2#0)-2,(w2#1)-1,(w2#2)-1,(w2#3)-1,(w2#4)-1,
                          (w2#5)-1,(w2#6)-1,(w2#7)}) else (  
                          if dd1*2 - (dot(w2,{1,1,1,1,1,0,0,0})) <= 0 then (  
                             dd1=dd1-2;  
                             w2 = {(w2#0)-1,(w2#1)-1,(w2#2)-1,(w2#3)-1,(w2#4)-1,
                             (w2#5),(w2#6),(w2#7)}) else (  
                             if dd1 - (dot(w2,{1,1,0,0,0,0,0,0}))< 0 then (  
                                dd1=dd1-1;  
                                w2 = {(w2#0)-1,(w2#1)-1,(w2#2),(w2#3),(w2#4),
                                (w2#5),(w2#6),(w2#7)}) else (  
                                myflag2=2))))))));
         if myflag2==1 then myker=0 else (  
            if dd1 - (dot(w2,{1,1,0,0,0,0,0,0})) == 0 then   
            myker=homcompdim({dd1-1,{(w2#0)-1,(w2#1),(w2#2),(w2#3),
               (w2#4),(w2#5),(w2#6),(w2#7)}})+  
            homcompdim({dd1-1,{(w2#0),(w2#1)-1,(w2#2),(w2#3),(w2#4),(w2#5),
               (w2#6),(w2#7)}}) else (  
               www={3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,
               3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,(w2#7)};  
               if {8*(w2#7)+3,www}=={dd1,w2} then myker=(w2#7)+1 else (
                  if homcompdim({dd1+1,w2})>3*(homcompdim({dd1,w2})) then 
                  myker=0 else myker=3*(homcompdim({dd1,w2}))-homcompdim({dd1+1,w2})))); 
         v2=join(v2,{myker}));  
      d1=d1+1);  
   scan(#v0, i->( -- this scan computes v3 from v2 and v1
         if i<2 then v3=join(v3,{0}) else (  
            if v2#(i-1)>-1 then v3=join(v3,{v1#i-3*(v1#(i-1))+v2#(i-1)}) else (  
               if v1#i-3*(v1#(i-1))>=0 then 
               v3=join(v3,{v1#i-3*(v1#(i-1))}) else v3=join(v3,{0})))));  
   scan(#v0, i->( -- this scan computes v4 from v3 and v1
         if i<3 then v4=join(v4,{0}) else 
         v4=join(v4,{v3#i-v1#i+3*(v1#(i-1))-3*(v1#(i-2))+v1#(i-3)})));    
   << "The output matrix has four columns. Column 1 indicates" << endl;
   << "each degree from alpha-2 (where alpha is the least" << endl;
   << "degree t such that I_t > 0 for the fat points ideal I)" << endl;
   << "to tau+2 (where tau is the least degree t such that the points" << endl;
   << "impose independent conditions in all degrees t or bigger)." << endl;
   << "Column 2 gives the value dim I_t of the Hilbert function" << endl;
   << "in each degree t listed in column 1. The resolution of I" << endl;
   << "is of the form 0 -> F_1 -> F_0 -> I -> 0, where F_1 and" << endl;
   << "F_0 are free S=k[P2] modules. Thus F_0=oplus_t S[-t]^{n_t}" << endl;
   << "and F_1=oplus_t S[-t]^{s_t} for integers s_t and n_t." << endl;
   << "Columns 3 and 4 give the values of n_t and s_t in each degree t" << endl;
   << "listed in column 1 (n_t and s_t are 0 in all other degrees)." << endl;
   transpose(matrix({v0,v1,v3,v4}))))  

-- This is the same as findres except suppresses printed output

findresReticent = (l) -> (
if #l>8 then (
   << "This script works only for up to 8 points." << endl;
   << "Please try again with an input list of at most 8 integers." << endl)
else (
   i:=0;
   myflag2:=0;
   w2:={};
   dd1:=0;
   myker:=0;
   www:={};
   ww:=l;
-- the list l of multiplicities is, for simplicity, extended if need be
-- so that it has 8 elements.
   while(#ww < 8) do ww=join(ww,{0});
   n:=#ww;
   n=n-1;
   ww=zr(ww);  -- zero out negative elements of the list
   a1:=findalpha(ww);  -- find alpha, the least degree t such that I_t \ne 0
   d1:=a1-2;
   tau:=findtau(ww);
   v4:={}; -- list of number of syzygies in each degree t listed in v0
   v3:={}; -- list of dim of coker of \mu_t in each degree t listed in v0
   v2:={}; -- list of dim of ker of \mu_t in each degree t listed in v0
   v1:={}; -- list of Hilbert function values for each degree listed in v0
   v0:={}; -- list of degrees from alpha-2 to tau+2, where tau is the least 
         -- degree such that the fat points impose independent conditions
   while (d1 <= tau+2) do ( -- loop from alpha-2 to tau+2, computing v0, v1 and v2
-- append the current degree d1 to the list of degrees
      v0=join(v0,{d1});
-- append the value of the hilbert function in degree d1
-- to the list of values of the hilbert function
      v1=join(v1,{homcompdim(fundom({d1,ww}))});
-- now compute and append to the list v2 the dimension myker of the 
-- kernel of the map \mu_t : I_t\otimes k[P2]_1 \to I_{t+1} where t=d1 
-- and I is the ideal of the fat points subscheme
      if d1<a1 then (v2=join(v2,{0}));  -- d1<a1 means I_{d1}=0 so myker=0
      if d1>=a1 then (     -- for d1>=a1, compute myker 
         myflag2=0;  
         w2=ww;
         dd1=d1;  
         while(myflag2==0) do ( -- this loop implements the main theorem of [FHH]
                                -- which gives an algorithm for computing myker
            w2=zr(w2);  
            w2=prmt(w2);  
            if homcompdim({dd1,w2})==0 then myflag2=1 else (  
              if dd1*6 - (dot(w2,{3,2,2,2,2,2,2,2})) <= 2 then (  
                 dd1=dd1-6;  
                 w2 = {(w2#0)-3,(w2#1)-2,(w2#2)-2,(w2#3)-2,(w2#4)-2,
                 (w2#5)-2,(w2#6)-2,(w2#7)-2}) else (  
                 if dd1*5 - (dot(w2,{2,2,2,2,2,2,1,1})) <= 1 then (  
                    dd1=dd1-5;  
                    w2 = {(w2#0)-2,(w2#1)-2,(w2#2)-2,(w2#3)-2,(w2#4)-2,
                    (w2#5)-2,(w2#6)-1,(w2#7)-1}) else (  
                    if dd1*4 - (dot(w2,{2,2,2,1,1,1,1,1})) <= 1 then (  
                       dd1=dd1-4;  
                       w2 = {(w2#0)-2,(w2#1)-2,(w2#2)-2,(w2#3)-1,(w2#4)-1,
                       (w2#5)-1,(w2#6)-1,(w2#7)-1}) else (  
                       if dd1*3 - (dot(w2,{2,1,1,1,1,1,1,0})) <= 0 then (  
                          dd1=dd1-3;  
                          w2 = {(w2#0)-2,(w2#1)-1,(w2#2)-1,(w2#3)-1,(w2#4)-1,
                          (w2#5)-1,(w2#6)-1,(w2#7)}) else (  
                          if dd1*2 - (dot(w2,{1,1,1,1,1,0,0,0})) <= 0 then (  
                             dd1=dd1-2;  
                             w2 = {(w2#0)-1,(w2#1)-1,(w2#2)-1,(w2#3)-1,(w2#4)-1,
                             (w2#5),(w2#6),(w2#7)}) else (  
                             if dd1 - (dot(w2,{1,1,0,0,0,0,0,0}))< 0 then (  
                                dd1=dd1-1;  
                                w2 = {(w2#0)-1,(w2#1)-1,(w2#2),(w2#3),(w2#4),
                                (w2#5),(w2#6),(w2#7)}) else (  
                                myflag2=2))))))));
         if myflag2==1 then myker=0 else (  
            if dd1 - (dot(w2,{1,1,0,0,0,0,0,0})) == 0 then   
            myker=homcompdim({dd1-1,{(w2#0)-1,(w2#1),(w2#2),(w2#3),
               (w2#4),(w2#5),(w2#6),(w2#7)}})+  
            homcompdim({dd1-1,{(w2#0),(w2#1)-1,(w2#2),(w2#3),(w2#4),(w2#5),
               (w2#6),(w2#7)}}) else (  
               www={3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,
               3*(w2#7)+1,3*(w2#7)+1,3*(w2#7)+1,(w2#7)};  
               if {8*(w2#7)+3,www}=={dd1,w2} then myker=(w2#7)+1 else (
                  if homcompdim({dd1+1,w2})>3*(homcompdim({dd1,w2})) then 
                  myker=0 else myker=3*(homcompdim({dd1,w2}))-homcompdim({dd1+1,w2})))); 
         v2=join(v2,{myker}));  
      d1=d1+1);  
   scan(#v0, i->( -- this scan computes v3 from v2 and v1
         if i<2 then v3=join(v3,{0}) else (  
            if v2#(i-1)>-1 then v3=join(v3,{v1#i-3*(v1#(i-1))+v2#(i-1)}) else (  
               if v1#i-3*(v1#(i-1))>=0 then 
               v3=join(v3,{v1#i-3*(v1#(i-1))}) else v3=join(v3,{0})))));  
   scan(#v0, i->( -- this scan computes v4 from v3 and v1
         if i<3 then v4=join(v4,{0}) else 
         v4=join(v4,{v3#i-v1#i+3*(v1#(i-1))-3*(v1#(i-2))+v1#(i-3)})));   
   transpose(matrix({v0,v1,v3,v4}))))  


-- findhilb: computes e(F_t(Z)) for alpha-1<= t <=tau+1, which gives 
-- a lower bound for the SHGH conjectural hilbert function
-- for a fat points subscheme involving general points of P2.
-- Call it as findhilb({m_1,...,m_n}) for integers m_i
-- specifying the multiplicities of the fat points in Z.
-- The conjecture is known to be correct for n<=9.

findhilb = (l) -> (
ww:=l;
if #l<3 then ww=join(l,{0,0,0});
n:=#ww;
n=n-1;
ww=zr(ww);
a1:=findalpha(ww);
tau:=findtau(ww);
d1:=a1-1;
<< "The output gives dim I_t, computed in degrees t from alpha(I)-1 to " << endl; 
<< "reg(I), where tau = reg(I)-1 is least degree such that" << endl;
<< "hilbert function of I equals hilbert polynomial of I." << endl;
if n>9 then (
   << "When more than 9 multiplicities are input," << endl;
   << "the output is a lower bound for dim I_t, which by the" << endl;
   << "SHGH conjecture should equal dim I_t." << endl);
<< endl;
<< " t      dim I_t" << "  (tau = " << tau << ")" << endl;
while (d1 <= tau+1) do (
   << " " << d1 << "     " << homcompdim(fundom({d1,ww})) << endl;
   d1=d1+1))

--input: l={n,m}, n = number of points, m = uniform multiplicity
--output: various bounds on alpha and tau

unifbounds = (l) -> (
n:=l#0;
m:=l#1;
ba:=bestrda(n);
bb:=bestrdb(n);
ea:=uniffindalpha(l);
t:=0;
<< "number of general points n of P2: " << n  << endl;
<< "multiplicity m of each point: " << m << endl;
<< endl;
if n<= 9 then (<< "Value of alpha: " << ea << endl)
else (
   << "Expected value of alpha (via SHGH conjecture): " << ea << endl;
   << " Note: The SHGH conjectural value of alpha is an upper bound." << endl);
<< "Lower Bounds on alpha:" << endl;
<< "  Roe's, via unloading: " << unifroealpha(l) << endl;
tmp:=unifezbhalpha(l);
<< "  Harbourne's, via Cor IV.i.2(a, b), using r="<<tmp#1<<" and d=";
<<tmp#2<<": "<< tmp#0 << endl;
t=ea;
while(t==unifbhalpha(l,ba#0,ba#1,t)) do t=t-1;
<< "  Harbourne's, via unloading, using r="<<ba#0<<" and d=";
<<ba#1<<": "<< t+1 << endl;
d:=0;
while(d*d < n) do d=d+1;
d=d-1;
r:=ceiling((n+d*d)/2);
<< "  Harbourne/Roe's first formula, using r="<<r<<" and d="<<
d<<": "<< ezunifHRalpha(l,r,d) << endl;
r=d*d;
while(r*r < n*d*d) do r=r+1;
r=r-1;
<< "  Harbourne/Roe's second formula, using r="<<r<<" and d="<<
d<<": "<<ezunifHRalphaB(l,r,d) << endl;
r=ba#0;
d=ba#1;
t=ea;
while(t==unifHRalpha(l,r,d,t)) do t=t-1;
t=t+1;
<<"  Harbourne/Roe's (via modified unloading), using r="<<r<<" and d="<<d;
<< ": "<<t<<endl;
r=bb#0;
d=bb#1;
t=ea;
while(t==unifHRalpha(l,r,d,t)) do t=t-1;
t=t+1;
<<"  Harbourne/Roe's (via modified unloading), using r="<<r<<" and d="<<d;
<< ": "<<t<<endl;
d=0;
while(d*d<=n) do d=d+1;
d=d-2;
r=d*(d+1);
t=ea;
while(t==unifHRalpha(l,r,d,t)) do t=t-1;
t=t+1;
<<"  Harbourne/Roe's (via modified unloading), using r="<<r<<" and d="<<d;
<< ": "<<t<<endl<<endl;
tt:=uniffindtau(l);
if n<= 9 then ( << "Value of tau: " << tt << endl)
else (
   << "Expected value of tau (via SHGH conjecture): " << tt << endl;
   << " Note: The SHGH conjectural value of tau is a lower bound." << endl);
<< "Upper Bounds on tau:" << endl;
<< "  Hirschowitz's: " << Hiuniftau(l) << endl;
<< "  Gimigliano's: " << Guniftau(l) << endl;
if n>4 then ( << "  Catalisano's: " << Cuniftau(l) << endl);
t=0;
while(t*(t+3)-n*m*(m+1) < 2*t*(m-1)-2) do t=t+1;
<< "  Ballico's: " << t << endl;
t=0;
while(9*(t+3)*(t+3) <= 10*n*(m+1)*(m+1)) do t=t+1;
<< "  Xu's: " << t << endl;
t=0;
tmp=0;
while(tmp*tmp < n) do tmp=tmp+1;
while(2*t < 2*m*tmp + tmp - 3) do t=t+1;
<< "  Harbourne/Holay/Fitchett's: " << t << endl;
<< "  Roe's, via unloading: " << unifroetau(l) << endl;
tmp=ezunifHRtau(l);
<< "  Harbourne/Roe's first formula, using r="<<tmp#1<<" and d="<<
tmp#2<<": "<< tmp#0 << endl;
<< "  Harbourne/Roe's second formula, using r="<<(tmp#2)*(tmp#2)<<" and d="<<
tmp#2<<": "<<ezunifHRtauB(l,(tmp#2)*(tmp#2),tmp#2) << endl;
r=ba#0;
d=ba#1;
t=unifHRtau(l,r,d,tt);
tmp=unifHRtau(l,bb#0,bb#1,tt);
if tmp<t then (
  t=tmp;
  r=bb#0;
  d=bb#1);
<<"  Harbourne/Roe's (via unloading), using r="<<r<<" and d="<<d<< ": "<<t<<endl;
if (ba#0)*(bb#0)<(ba#1)*(bb#1)*n then (
   t = -3 + ceiling((m+1)*(bb#1)*n/(bb#0));
   r=bb#0;
   d=bb#1) else (
   t = -3 + ceiling((m+1)*(ba#0)/(ba#1));
   r=ba#0;
   d=ba#1);
<< "  Via Ran's observation, and Harbourne's bound on alpha," << endl;
<< "      using r="<<r<<" and d="<<d<< ": "<<t<<endl)

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: various bounds on alpha and tau

bounds = (l) -> (
n:=#l;
ba:=bestrda(n);
bb:=bestrdb(n);
ea:=findalpha(zr(l));
t:=0;
tmp:=0;
<< "number of general points n of P2: " << n << endl;
<< "multiplicities of the points: " << l << endl;
<< endl;
if n<= 9 then ( << "Value of alpha: " << ea << endl)
else (
   << "Expected value of alpha (via SHGH conjecture): " << ea << endl;
   << " Note: The SHGH conjectural value of alpha is an upper bound." << endl);
<< "Lower bounds on alpha:" << endl;
<< "  Via Checking Psi: " << Psibound(l) << endl;
<< "  Roe's, via unloading: " << roealpha(l) << endl;
w:=ezbhalphaA(l);
<< "  Harbourne's, via Cor IV.i.2(a), using r="<<w#1<<", d=";
<<w#2<<": "<<w#0<<endl;
w=ezbhalphaB(l);
<< "  Harbourne's, via Cor IV.i.2(b), using r="<<w#1<<", d=";
<<w#2<<": "<<w#0<<endl;
w=ezbhalphaD(l);
<< "  Harbourne's, via Cor IV.i.2(d), using r="<<w#1<<", d="<<w#2;
<<", and j="<<w#3<<": "<<w#0<<endl;
r:=ba#0;
d:=ba#1;
t=ea;
while(t==HRalpha(l,r,d,t)) do t=t-1;
t=t+1;
tmp=ea;
while(tmp==HRalpha(l,bb#0,bb#1,tmp)) do tmp=tmp-1;
tmp=tmp+1;
if tmp>t then (
  t=tmp;
  r=bb#0;
  d=bb#1);
<<"  Harbourne/Roe's (via modified unloading), using r="<<r<<" and d="<<d;
<< ": "<<t<< endl << endl;
tt:=findtau(l);
if n<= 9 then ( << "Value of tau: " << tt << endl)
else (
   << "Expected value of tau (via SHGH conjecture): " << tt << endl;
   << " Note: The SHGH conjectural value of tau is a lower bound." << endl);
<< "Upper bounds on tau:" << endl;
<< "  Hirschowitz's: " << Hitau(l) << endl;
<< "  Gimigliano's: " << Gtau(l) << endl;
nn:=0;
i:=0;
w=prmt(zr(l));
scan(#l, i->(if w#i >0 then nn=nn+1));
if nn>4 then ( << "  Catalisano's: " << Ctau(l) << endl);
<< "  Roe's, via unloading: " << roetau(l) << endl;
r=ba#0;
d=ba#1;
t=HRtau(l,r,d,tt);
tmp=HRtau(l,bb#0,bb#1,tt);
if tmp<t then (
  t=tmp;
  r=bb#0;
  d=bb#1);
<< "  Harbourne/Roe's (via unloading), using r="<<r<<" and d="<<d<<": "<<t<<endl<<endl;
<<"One more bound on alpha; this one can be slow since it tries all r and d." << endl;
w=bestbhalpha(l,ea);
<< "  Harbourne's alpha lower bound (via unloading), using r=" << w#1 << " and d="; 
<< w#2 << ": " << w#0 << endl)

--decomp: prints a decomposition F=H+N for any divisor class F in Psi
--as described in the print statements of the script.
--input: decomp(l), where l={d,{m1,...,mn}}, signifying the divisor
--class F = dE_0-(m_1E_1+...+m_nE_n).

decomp = (l) -> (
<< "Let Psi be the subsemigroup of divisor classes generated by " << endl;
<< "exceptional classes and by -K. For any divisor class F, this " << endl;
<< "script determines if F is in Psi, and if so gives a decomposition" << endl;
<< "F=H+N, where N is a sum of exceptionals orthogonal to each other and to H" << endl;
<< "and H is in Psi but H.E >= 0 for all exceptionals E. The point of this" << endl;
<< "is that dim |F| = dim |H|, and conjecturally dim |H| = (H^2-H.K)/2." << endl << endl;
i:=0;
j:=0;
w:=l;
nef:={};
v:={};
ww:={};
ex:={};
mult:=0;
tmp:=fundom(l);
if tmp#0<tmp#1#0 or tmp#0<0 then (<< "Your class is not in Psi." << endl) else (
   if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
   d:=3*(#(w#1));  --define an element {d,v} in fundamental domain
   scan(#(w#1), i->(v=join(v,{(#(w#1)) - i})));
   ww=fundomboth(w,{d,v});
   if ww#0#1 == zr(ww#0#1) then (<< "N = 0" << endl ) else (
   << "N is a sum of the following fixed exceptional classes:" << endl;
   if (ww#0#0) - (ww#0#1#0) - (ww#0#1#1) < 0 then (
       scan(#(w#1), j->(if j<=1 then ex=join(ex,{1}) else ex=join(ex,{0})));
       ex={1,ex};
       ex=(fundomboth(ww#1,ex))#1;
       <<ex<<" is a fixed component of multiplicity ";
       <<(ww#0#1#0)+(ww#0#1#1)-(ww#0#0)<<endl;
ab=bundledecomp(ex);
<< "Its splitting constants: " << ab << endl;
       ex={(ww#0#0)-(ww#0#1#1),(ww#0#0)-(ww#0#1#0)};
       scan(#(w#1), j->(if j>1 then ex=join(ex,{ww#0#1#j})));
       ww={{2*(ww#0#0)-(ww#0#1#0)-(ww#0#1#1),ex},ww#1});
   scan(#(w#1), i->(
         if (ww#0#1)#i<0 then (
            ex={};
            mult=-(ww#0#1)#i;
            scan(#(w#1), j->(if j==i then ex=join(ex,{-1}) else ex=join(ex,{0})));
            ex={0,ex};
            ex=(fundomboth(ww#1,ex))#1;
            << ex << " is a fixed component of multiplicity " << mult << endl;
ab=bundledecomp(ex);
<< "Its splitting constants: " << ab << endl))));
nef=(fundomboth(ww#1,{ww#0#0,zr(ww#0#1)}))#1;
   << endl << "and H = " << nef << endl;
<< " which is Cremona equivalent to " << fundom(nef) << "." << endl))

decompnonverbose = (l) -> (
i:=0;
j:=0;
w:=l;
nef:={};
v:={};
ww:={};
ex:={};
mult:=0;
tmp:=fundom(l);
if tmp#0<tmp#1#0 or tmp#0<0 then (<< "Your class is not in Psi." << endl) else (
   if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
   d:=3*(#(w#1));  --define an element {d,v} in fundamental domain
   scan(#(w#1), i->(v=join(v,{(#(w#1)) - i})));
   ww=fundomboth(w,{d,v});
   if ww#0#1 == zr(ww#0#1) then ( ) else (
   if (ww#0#0) - (ww#0#1#0) - (ww#0#1#1) < 0 then (
       scan(#(w#1), j->(if j<=1 then ex=join(ex,{1}) else ex=join(ex,{0})));
       ex={1,ex};
       ex=(fundomboth(ww#1,ex))#1;
ab=bundledecompnonverbose(ex);
       ex={(ww#0#0)-(ww#0#1#1),(ww#0#0)-(ww#0#1#0)};
       scan(#(w#1), j->(if j>1 then ex=join(ex,{ww#0#1#j})));
       ww={{2*(ww#0#0)-(ww#0#1#0)-(ww#0#1#1),ex},ww#1});
   scan(#(w#1), i->(
         if (ww#0#1)#i<0 then (
            ex={};
            mult=-(ww#0#1)#i;
            scan(#(w#1), j->(if j==i then ex=join(ex,{-1}) else ex=join(ex,{0})));
            ex={0,ex};
            ex=(fundomboth(ww#1,ex))#1;
ab=bundledecompnonverbose(ex)))));
nef=(fundomboth(ww#1,{ww#0#0,zr(ww#0#1)}))#1))

-- This was obtained from the code for decomp.
-- It just returns the nef part.

nefpart = (l) -> (
i:=0;
j:=0;
w:=l;
nef:={};
v:={};
ww:={};
ex:={};
mult:=0;
tmp:=fundom(l);
if tmp#0<tmp#1#0 or tmp#0<0 then ( ) else (
   if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
   d:=3*(#(w#1));  --define an element {d,v} in fundamental domain
   scan(#(w#1), i->(v=join(v,{(#(w#1)) - i})));
   ww=fundomboth(w,{d,v});
   if ww#0#1 == zr(ww#0#1) then ( ) else (
   if (ww#0#0) - (ww#0#1#0) - (ww#0#1#1) < 0 then (
       scan(#(w#1), j->(if j<=1 then ex=join(ex,{1}) else ex=join(ex,{0})));
       ex={1,ex};
       ex=(fundomboth(ww#1,ex))#1;
       ex={(ww#0#0)-(ww#0#1#1),(ww#0#0)-(ww#0#1#0)};
       scan(#(w#1), j->(if j>1 then ex=join(ex,{ww#0#1#j})));
       ww={{2*(ww#0#0)-(ww#0#1#0)-(ww#0#1#1),ex},ww#1}));
nef=(fundomboth(ww#1,{ww#0#0,zr(ww#0#1)}))#1;
nef))

decompDeluxe = (l,ddd,nnn) -> (
ab:={};
i:=0;
j:=0;
gens:=0;
nef:={};
w:=l;
v:={};
ww:={};
ex:={};
mult:=0;
tmp:=fundom(l);
if tmp#0<tmp#1#0 or tmp#0<0 then (<< "Your class is not in Psi." << endl) else (
   if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
   d:=3*(#(w#1));  --define an element {d,v} in fundamental domain
   scan(#(w#1), i->(v=join(v,{(#(w#1)) - i})));
   ww=fundomboth(w,{d,v});
   if ww#0#1 == zr(ww#0#1) then ( ) else (
   if (ww#0#0) - (ww#0#1#0) - (ww#0#1#1) < 0 then (
       scan(#(w#1), j->(if j<=1 then ex=join(ex,{1}) else ex=join(ex,{0})));
       ex={1,ex};
       ex=(fundomboth(ww#1,ex))#1;
ab=bundledecomp(ex);
mult=(ww#0#1#0)+(ww#0#1#1)-(ww#0#0);
-- nef=(fundomboth(ww#1,{ww#0#0,zr(ww#0#1)}))#1;
-- << "1." << ab << endl;
-- << "1." << ex << endl;
-- << "1." << ddd << endl;
-- << "1." << {nef#0+1,nef#1} << endl;
-- << "1." << homcompdim({nef#0+1,nef#1}) << endl;
if mult > ex#0 then mult = ex#0;
-- << "1." << mult << endl;
if mult > ab#0 + 1 then gens=gens+floor((mult-ab#0)*(mult-ab#0-1)/2);
if mult > ab#1 + 1 then gens=gens+floor((mult-ab#1)*(mult-ab#1-1)/2);
-- gens=gens+ddd-homcompdim({nef#0+1,nef#1});
       ex={(ww#0#0)-(ww#0#1#1),(ww#0#0)-(ww#0#1#0)};
       scan(#(w#1), j->(if j>1 then ex=join(ex,{ww#0#1#j})));
       ww={{2*(ww#0#0)-(ww#0#1#0)-(ww#0#1#1),ex},ww#1});
   scan(#(w#1), i->(
         if (ww#0#1)#i<0 then (
            ex={};
            mult=-(ww#0#1)#i;
            scan(#(w#1), j->(if j==i then ex=join(ex,{-1}) else ex=join(ex,{0})));
            ex={0,ex};
            ex=(fundomboth(ww#1,ex))#1;
ab=bundledecomp(ex);
-- << "2." << l << endl;
-- << "2." << ab << endl;
-- << "2." << ex << endl;
-- << "2." << ddd << endl;
if mult > ex#0 then mult = ex#0;
-- << "2." << mult << endl;
if mult > ab#0 + 1 then gens=gens+floor((mult-ab#0)*(mult-ab#0-1)/2);
if mult > ab#1 + 1 then gens=gens+floor((mult-ab#1)*(mult-ab#1-1)/2)))));
-- << "2." << gens << endl;
-- << l << endl;
-- << ddd << endl;
-- << nnn << endl;
-- << gens << endl;
-- << "2." << {nef#0+1,nef#1} << endl;
-- << "2." << homcompdim({nef#0+1,nef#1}) << endl;
gens=gens+ddd-nnn;
-- << gens << endl;
gens))

decompDeluxenonverbose = (l,ddd,nnn) -> (
ab:={};
i:=0;
j:=0;
gens:=0;
nef:={};
w:=l;
v:={};
ww:={};
ex:={};
mult:=0;
tmp:=fundom(l);
if tmp#0<tmp#1#0 or tmp#0<0 then (<< "Your class is not in Psi." << endl) else (
   if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
   d:=3*(#(w#1));  --define an element {d,v} in fundamental domain
   scan(#(w#1), i->(v=join(v,{(#(w#1)) - i})));
   ww=fundomboth(w,{d,v});
   if ww#0#1 == zr(ww#0#1) then ( ) else (
   if (ww#0#0) - (ww#0#1#0) - (ww#0#1#1) < 0 then (
       scan(#(w#1), j->(if j<=1 then ex=join(ex,{1}) else ex=join(ex,{0})));
       ex={1,ex};
       ex=(fundomboth(ww#1,ex))#1;
ab=bundledecompnonverbose(ex);
mult=(ww#0#1#0)+(ww#0#1#1)-(ww#0#0);
-- nef=(fundomboth(ww#1,{ww#0#0,zr(ww#0#1)}))#1;
-- << "1." << ab << endl;
-- << "1." << ex << endl;
-- << "1." << ddd << endl;
-- << "1." << {nef#0+1,nef#1} << endl;
-- << "1." << homcompdim({nef#0+1,nef#1}) << endl;
if mult > ex#0 then mult = ex#0;
-- << "1." << mult << endl;
if mult > ab#0 + 1 then gens=gens+floor((mult-ab#0)*(mult-ab#0-1)/2);
if mult > ab#1 + 1 then gens=gens+floor((mult-ab#1)*(mult-ab#1-1)/2);
-- gens=gens+ddd-homcompdim({nef#0+1,nef#1});
       ex={(ww#0#0)-(ww#0#1#1),(ww#0#0)-(ww#0#1#0)};
       scan(#(w#1), j->(if j>1 then ex=join(ex,{ww#0#1#j})));
       ww={{2*(ww#0#0)-(ww#0#1#0)-(ww#0#1#1),ex},ww#1});
   scan(#(w#1), i->(
         if (ww#0#1)#i<0 then (
            ex={};
            mult=-(ww#0#1)#i;
            scan(#(w#1), j->(if j==i then ex=join(ex,{-1}) else ex=join(ex,{0})));
            ex={0,ex};
            ex=(fundomboth(ww#1,ex))#1;
ab= bundledecompnonverbose(ex);
-- << "2." << l << endl;
-- << "2." << ab << endl;
-- << "2." << ex << endl;
-- << "2." << ddd << endl;
if mult > ex#0 then mult = ex#0;
-- << "2." << mult << endl;
if mult > ab#0 + 1 then gens=gens+floor((mult-ab#0)*(mult-ab#0-1)/2);
if mult > ab#1 + 1 then gens=gens+floor((mult-ab#1)*(mult-ab#1-1)/2)))));
-- << "2." << gens << endl;
-- << l << endl;
-- << ddd << endl;
-- << nnn << endl;
-- << gens << endl;
-- << "2." << {nef#0+1,nef#1} << endl;
-- << "2." << homcompdim({nef#0+1,nef#1}) << endl;
gens=gens+ddd-nnn;
-- << gens << endl;
gens))


-- input: list l of multiplicities for fat points Z
-- output: least t such that F_t(Z) is in Psi

Psibound = (l) -> (
t:=findalpha(l);
tmp:=fundom({t,l});
while(tmp#0>=tmp#1#0 and tmp#0>=0) do (t=t-1;
   tmp=fundom({t,l}));
t+1) 

-- prmt: arranges the elements of the list l={m_1,...,m_n} in descending order
-- Call it as prmt(l) where l is a list of integers.

prmt = (l) -> (  
(prmtboth(l,l))#0)

-- prmtboth: arranges the elements of the list l1 in descending order,
-- and applies the same permutation to l2
-- Call it as prmt(l1,l2) where l1 and l2 are lists of integers.

prmtboth = (l1,l2) -> (  
tmpv1:=l1;
tmpv2:=l2;
v1:=l1;
v2:=l2;
i:=0;
j:=0;
k:=0;
scan(#l1, i->(scan(#l1, j->(
            if tmpv1#i < tmpv1#j then (
               if i < j then (
                  k=-1;
                  v1={};
                  v2={};
                  while(k<#l1-1) do (k=k+1;
                     if k==i then (v1=join(v1,{tmpv1#j});
                        v2=join(v2,{tmpv2#j})) else (
                        if k==j then (v1=join(v1,{tmpv1#i});
                           v2=join(v2,{tmpv2#i})) else (v1=join(v1,{tmpv1#k});
                              v2=join(v2,{tmpv2#k}))));
                  tmpv1=v1;
                  tmpv2=v2))))));
{v1,v2})

-- zr: replaces negative values in a list l by zeroes.
-- Call it as zr(l) where l is a list of integers.

zr = (l) -> (    
v:={};
i:=0;
scan(#l, i->(
      if l#i<0 then v=join(v,{0}) else v=join(v,{l#i})));
v)


-- fundomboth: applies fundom to l1 to reduce l1 to fundamental
-- domain of a certain group operation, and applies the same
-- group operation g to l2. If l2 starts out in the fundamental domain,
-- and {l1',l2'}=fundomboth(l1,l2), then {l2,g^{-1}l}=fundomboth(l2',l).
-- This allows one to compute the action of g^{-1}.

fundomboth = (l1,l2) -> (
w1:=l1; 
w2:=l2;
v:={}; 
if #l1#1<3 then w1={l1#0,join(l1#1,{0,0,0})};
if #l2#1<3 then w2={l2#0,join(l2#1,{0,0,0})};
v=prmtboth(w1#1,w2#1);
w1={w1#0,v#0};
w2={w2#0,v#1};
while ((w1#0 < w1#1#0 + w1#1#1 + w1#1#2) and (w1#0 >= 0)) do (
   w1=quad(w1);
   w2=quad(w2);
   v=prmtboth(w1#1,w2#1);
   w1={w1#0,v#0};
   w2={w2#0,v#1});
{w1,w2})

-- fundom: Call it as fundom({d,{m1,...,mn}}). The output is a new 
-- list {d',{m1',...,mn'}}; the class dE_0-(m_1E_1+...+m_nE_n) is
-- equivalent via Cremona transformations to d'E_0-(m_1'E_1+...+m_n'E_n),
-- where d' is the smallest positive integer possible.

fundom = (l) -> (
(fundomboth(l,l))#0)

-- homcompdim: computes e(F_t(Z)), the expected dimension of a component I_d 
-- of a fat points ideal I corresponding to a fat point subscheme Z of general
-- points taken with multiplicities m_1, ..., m_n. Call it as 
-- homcompdim({d,{m_1,...,m_n}}); the output is the SHGH conjectural 
-- dimension of I_d, which is the actual dimension if n < 10.

homcompdim = (l) -> (
h:=0;
i:=0;
w:=l; 
if #l#1<3 then w={l#0,join(l#1,{0,0,0})};
w=fundom(w);
d:=w#0;
w=fundom({d,zr(w#1)});
d=w#0;
v:=zr(w#1);
if d<0 then h=0 else (
   tmp:=0;
   scan(#v, i->(tmp = v#i*v#i+v#i+tmp));
   h=floor((d*d+3*d+2-tmp)/2);
   if h < 0 then h=0);
h)

-- findalpha: find alpha, the least degree t such that
-- I_t \ne 0, where I is the ideal corresponding to n general
-- points taken with multiplicities m_1, ...,m_n. Call it as 
-- findalpha({m_1,...,m_n}). The output is the SHGH conjectural value 
-- of alpha, which is the actual value if n < 10 and an upper bound otherwise.

findalpha = (l) -> (
i:=1;
w:=prmt(zr(l)); 
if #l<3 then w=join(l,{0,0,0});
d:=w#0;                -- alpha is at least the max mult
if (#w)<9 then (       -- if n<=8, to speed things, make an estimate
  while(i < (#w)) do (d=d+w#i;
     i=i+1);
  d=ceiling(d/3);
  if d < w#0 then d = w#0);
while (homcompdim({d,w}) < 1) do d=d+1;
d)

-- dot: computes a dot product of two lists l1 and l2 (of equal length)
-- of integers. Call it as dot(l1,l2).

dot = (l1,l2) -> (
i:=0;
dottot:=0;
scan(#l1, i->(dottot = dottot + (l1#i)*(l2#i)));
dottot)

-- input: l={m1,...,mn}, n >=1 (number of points), m1,... (the multiplicities)
-- output: the SHGH conjectured value of tau; this is the actual value if 
-- n < 10, and a lower bound otherwise.

findtau = (l) -> (
t:=findalpha(l);
if t > 0 then t = t-1;  -- tau is at least alpha - 1
n:=#l;
v:=l;
p:=dot(v,v);
K:={};
j:=0;
q:=0;
scan(#l, j->(q=q+l#j));
while(2*(homcompdim(join({t},{v}))) > t*t-p+3*t-q+2) do t=t+1;
t)

-- input: positive integer n
-- output: {r,d}, where r^2>=d^2n, n>=r, and nd/r is as big as possible

bestrda = (n) -> (
rootn:=0;
while(rootn*rootn<=n) do rootn=rootn+1;
rootn=rootn-1;
d:=1;
r:=0;
if rootn*rootn==n then r=rootn else r=rootn+1;
tmpr:=1;
tmpd:=1;
while(tmpd<=rootn) do (
   tmpr=tmpd*rootn;
   while(tmpr*tmpr<tmpd*tmpd*n) do tmpr=tmpr+1;
   if tmpr*d < tmpd*r then (
      r=tmpr;
      d=tmpd);
   tmpd=tmpd+1);
{r,d})

-- input: positive integer n
-- output: {r,d}, where r^2<=d^2n, n>=r, and r/d is as big as possible

bestrdb = (n) -> (
rootn:=0;
while(rootn*rootn<n) do rootn=rootn+1;
d:=1;
r:=0;
if rootn*rootn==n then r=rootn else r=rootn-1;
tmpr:=1;
tmpd:=1;
while(tmpd<=rootn) do (
   tmpr=tmpd*rootn;
   while(tmpr*tmpr>tmpd*tmpd*n) do tmpr=tmpr-1;
   if tmpr>n then tmpr=n;
   if tmpr*d > tmpd*r then (
      r=tmpr;
      d=tmpd);
   tmpd=tmpd+1);
{r,d})

-- input: l={m1,...,mn}, n= number of points, mi = multiplicity of ith point
-- output: Roe's algorithmic lower bound on alpha

roealpha = (l) -> (
i:=0;
v={};
w2:={};
w:=l;
i1:=2;
if #l<3 then w=join(w,{0,0,0});
while(i1<#w) do (
   v={1};
   scan(#w, i->(if i>0 then (if i<=i1 then v=join(v,{-1}) else v=join(v,{0}))));
   w=zr(prmt(w));
   while(dot(w,v)<0) do (
      w2={};
      scan(#w, i->(w2 = join(w2,{w#i+v#i})));
      w=zr(prmt(w2)));
   i1=i1+1);
w#0)

-- input: l={n,m}, n >=1 (number of points), m >=1 (uniform multiplicity)
-- output: Roe's algorithmic lower bound on alpha

unifroealpha = (l) -> (
i1:=0;
intchk:=0;
n:=l#0;
m:=l#1;
roebnd:=m;
q:=n-1; -- q keeps track during unloading of number of 
        -- points after the first with maximum multiplicity
if n>2 then (
   i1=2;
   while(i1<n) do (
      if q<i1 then intchk=i1*m-(i1-q) else intchk=i1*m;
      while(roebnd<intchk) do (
         roebnd=roebnd+1;
         if q<i1 then (q=n-i1+q-1;
            m=m-1) else (q=q-i1;
            if q==0 then (m=m-1;
               q=n-1));
         if q<i1 then intchk=i1*m-(i1-q) else intchk=i1*m);
      i1=i1+1)) else roebnd = m;
roebnd)

-- input: l={m1,...,mn}, n= number of points, mi = multiplicity of ith point
-- output: Roe's algorithmic upper bound on tau

roetau = (l) -> (
i:=0;
vv={};
vv1:={};
w2:={};
ww:=l;
i1:=1;
while(i1<#ww-1) do (
   vv={1};
   scan(#ww, i->(if i>0 then (if i-1<=i1 then vv=join(vv,{-1}) else vv=join(vv,{0}))));
   vv1={1};
   scan(#ww, i->(if i>0 then (if i<=i1 then vv1=join(vv1,{-1}) else vv1=join(vv1,{0}))));
   ww=zr(prmt(ww));
   while((dot(ww,vv)) < -1) do (
      w2={};
      scan(#ww, i->(w2 = join(w2,{ww#i+vv1#i})));
      ww=zr(prmt(w2)));
   i1=i1+1);
ww#0+ww#1-1)

-- input: l={n,m}, n >=1 (number of points), m >=1 (uniform multiplicity)
-- output: Roe's algorithmic upper bound on tau

unifroetau = (l) -> (
i:=1;
s:=0;
n:=l#0;
m1:=l#1;
m2:=0;
if n>1 then m2=m1;
n2:=n-1; -- n2 keeps track during unloading of number of 
         -- points with multiplicity m2
while(i < (n-1)) do (
   if (i+1) <= n2 then s=(i+1)*m2 else s=(i+1)*m2+n2-i-1;
   while((m1-s) < -1) do (
      m1=m1+1;
      if i < n2 then n2=n2-i else (
         n2=n-i+n2-1;
         m2=m2-1);
      if (i+1) <= n2 then s=(i+1)*m2 else s=(i+1)*m2+n2-i-1);
   i=i+1);
m1+m2-1)

-- input: l={n,m}, n >=1 (number of points), m >=1 (uniform multiplicity)
-- output: list {i,r,d}, with i being Harbourne's easy lower bound on 
-- alpha (via Cor IV.i.2 (a), (b)) computed using the best r and d

unifezbhalpha = (l) -> (
w:=bestrda(l#0);
i:=ceiling((l#1)*(l#0)*(w#1)/(w#0));  -- compute mnd/r rounded up
r:=w#0;
d:=w#1;
w=bestrdb(l#0);
j:=ceiling((l#1)*(w#0)/(w#1));  -- compute mr/d rounded up
if i<j then (i=j;
  r=w#0;
  d=w#1);
{i,r,d})

-- input: ({n,m},r,d), n >=1 (number of points), m >=1 (uniform multiplicity)
-- where r <= n, and 2r >= n+d^2
-- output: lower bound on alpha via formula of [HR] 

ezunifHRalpha = (l,r,d) -> (
n:=l#0;
m:=l#1;
t:=0;
q:=ceiling(n*m/r)-1;
while(((t+2)*(t+1)<=2*(m*n-r*q)) and t<d) do t=t+1;
t+q*d)

-- input: ({n,m},r,d), n >=1 (number of points), m >=1 (uniform multiplicity)
-- where r <= n, and d(d+1)/2 <= r <= nd^2
-- output: lower bound on alpha via formula of [HR2]. 

ezunifHRalphaB = (l,r,d) -> (
n:=l#0;
m:=l#1;
g:=(d-1)*(d-2)/2;
tmp:=floor((m*r+g-1)/d);
t:=ceiling(n*m/r)-1;
rr:=m*n-r*t;
s:=0;
while(((s+1)*(s+2) <= 2*rr) and s<d) do s=s+1;
s=s-1;
t=s+t*d;
if tmp<t then t=tmp;
t+1)

-- input: ({n,m},r,d), n >=1 (number of points), m >=1 (uniform multiplicity)
-- where r <= n, and r <= d^2
-- output: upper bound on tau via formula of [HR]. 

ezunifHRtauB = (l,r,d) -> (
n:=l#0;
m:=l#1;
g:=(d-1)*(d-2)/2;
q:=ceiling(n*m/r)-1;
rr:=m*n-r*q;
t:=q*d+ceiling((rr+g-1)/d);
tmp:=q*d+d-2;
if t<tmp then t=tmp;
t)

-- input: ({n,m},r,d,ea), n >=1 (number of points), m >=1 (uniform multiplicity)
-- ea is an estimate for alpha (for speed); it must be set equal to
-- a value no bigger than the eventual value of unifHRalpha (0, for example)
-- output: lower bound on alpha via modified unloading method of [HR], 
-- computed using given r and d

unifHRalpha = (l,r,d,ea) -> (
i:=ea-1;
n2:=0;  -- n2 keeps track of the number of points with maximum multiplicity
s:=0;
tmpi:=-1;
tmpm:=0;
g:=(d-1)*(d-2)/2;
while(tmpi<tmpm) do (
   i=i+1;
   tmpi=i;
   tmpm=l#1;
   n2=l#0;
   if r <= n2 then s=r*tmpm else s=r*tmpm-r+n2;
   while(((tmpi*d-s < g) and (tmpi >= d-2)) or (((tmpi+1)*(tmpi+2)<=2*s) and 
                                               (tmpi<d) and (tmpi>=0))) do (
         tmpi=tmpi-d;       
         if r<n2 then n2=n2-r else (          
             n2=l#0-r+n2;
             tmpm=tmpm-1;
             if tmpm <= 0 then (
                tmpm=0;
                n2=l#0));
         if r<=n2 then s=r*tmpm else s=r*tmpm-r+n2));
i)

-- input: ({m_1,...,m_n},r,d,ea), n >=1 (number of points), m_i >=1 (multiplicities)
-- ea is an estimate for alpha (for speed); it must be set equal to
-- a value no bigger than the eventual value of HRalpha (0, for example)
-- output: lower bound on alpha via modified unloading method of [HR], 
-- computed using given r and d

HRalpha = (l,r,d,ea) -> (
i:=ea-1;
j:=0;
n:=0;
g:=(d-1)*(d-2)/2;
v:=prmt(zr(l));
tmpv:=v;
ttmpv:={};
scan(#l, j->(if v#j>0 then n=n+1));
if n>0 then (
   ww:={};
   scan(#l, j->(if j<r then ww=join(ww,{1}) else ww=join(ww,{0})));
   tmpi:=-1;
   while(tmpi < tmpv#0) do (
      i=i+1;
      tmpi=i;
      tmpv=v;
      while(((tmpi*d-(dot(ww,tmpv))<g) and (tmpi >= d-2)) or 
            (((tmpi+1)*(tmpi+2)<=2*(dot(ww,tmpv))) and (tmpi<d) and (tmpi>=0))) do (
         tmpi=tmpi-d;
         ttmpv={};
         scan(#l, j->(ttmpv=join(ttmpv,{(tmpv#j)-(ww#j)})));
         tmpv=prmt(zr(ttmpv)))));
i)

-- input: ({n,m},r,d,ea), n >=1 (number of points), m >=1 (uniform multiplicity)
-- ea is an estimate for alpha (for speed); it must be set equal to
-- a value no bigger than the eventual value of unifbhalpha (0, for example)
-- output: Harbourne's algorithmic lower bound on alpha via unloading, 
-- using the given r and d. 

unifbhalpha = (l,r,d,ea) -> (
i:=ea-1;
n2:=0;  -- n2 keeps track of the number of points with maximum multiplicity
s:=0;
tmpi:=-1;
tmpm:=0;
while(tmpi<0) do (
   i=i+1;
   tmpi=i;
   tmpm=l#1;
   n2=l#0;
   if r <= n2 then s=r*tmpm else s=r*tmpm-r+n2;
   while((tmpi*d-s < 0) and (tmpi >= 0)) do (
      tmpi=tmpi-d;
      if r<n2 then n2=n2-r else (
         n2=l#0-r+n2;
         tmpm=tmpm-1;
         if tmpm <= 0 then (
            tmpm=0;
            n2=l#0));
      if r<=n2 then s=r*tmpm else s=r*tmpm-r+n2));
i)

-- input: ({n,m},r,d), n >=1 (number of points), m >=1 (uniform multiplicity)
-- where r <= n, and 2r >= n+d^2
-- output: lower bound on alpha via formula which agrees with that via unloading. 

ezunifBHalphaB = (l,r,d) -> (
n:=l#0;
m:=l#1;
q:=floor(n*m/r);
rr:=m*n-r*q;
t:=q-1+ceiling(rr/d);
tmp:=d*ceiling(n*m/r);
if tmp<t then t=tmp;
t+1)

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: list (aa,r,d), with aa being Harbourne's lower bound on 
-- alpha (via Cor IV.i.2(a)), computed using the best r and d

ezbhalphaA = (l) -> (
i:=0;
s:=0;
n:=0;
v:=prmt(zr(l));
scan(#l, i->(if v#i>0 then n=n+1));
w:=bestrda(n);
i=0;
while(i<n) do (
  s=s+v#i;
  i=i+1);
best:=ceiling(s*(w#1)/(w#0));
{best, w#0, w#1})

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: list (aa,r,d), with aa being Harbourne's lower bound on 
-- alpha (via Cor IV.i.2(b)), computed using the best r and d

ezbhalphaB = (l) -> (
i:=0;
s:=0;
n:=0;
v:=prmt(zr(l));
scan(#l, i->(if v#i>0 then n=n+1));
w:=bestrdb(n);
i=0;
while(i<n) do (
  s=s+v#i;
  i=i+1);
best:=ceiling(s*(w#0)/(n*(w#1)));
{best, w#0, w#1})

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: list (aa,r,d,j), with aa being Harbourne's lower bound on 
-- alpha (via Cor IV.i.2(d)), computed using the best r, d and j

ezbhalphaD = (l) -> (
w:={};
i:=0;
scan(#l, i->(if l#i>0 then w=join(w,{l#i})));
w=prmt(w);
n:=#w;
bnd:=0;
tmpbnd:=0;
r:=0;
d:=0;
j:=0;
tmpr:=0;
tmpd:=0;
tmpj:=0;
if n>0 then (
   while(tmpr<n) do (tmpr=tmpr+1;
      tmpd=0;
      while(tmpd*tmpd<tmpr) do (tmpd=tmpd+1;
         tmpj=0;
         while(tmpj<tmpd*tmpd) do (tmpj=tmpj+1;
             tmpbnd=lpa(w,tmpr,tmpd,tmpj);
             if tmpbnd>bnd then (
                 bnd=tmpbnd;
                 r=tmpr;
                 d=tmpd;
                 j=tmpj)))));
{bnd,r,d,j})

-- lpa computes the bound given in Cor IV.i.2(d); attempts
-- various solutions with the hope of approximating the optimal 
-- solution to the linear programming problem indicated by Thm IV.i.1.
-- lpa is called by ezbhalphaD

lpa = (l,r,d,j) -> (
i:=0;
n:=#l;
sum:=0;
sumb:=0;
bnd:=0;
if d*d >= r then  (
   scan(#l, i->(if i<r then sum=sum+l#i));
   bnd=ceiling(sum/d)) else (
   if j==0 then (
      scan(#l, i->(if i<d*d then sum=sum+l#i));
      bnd=ceiling(sum/d)) else (
      M:=floor((r-d*d)*(r-d*d+j)/j);
      scan(#l, i->(if i<d*d-j then sum=sum+l#i else (if i<M+r then sumb=sumb+l#i)));
      sumb=sumb*j/(r-d*d+j);
      if M<n-r then sumb=sumb+(l#(M+r))*(r-d*d-j*M/(r-d*d+j));
      bnd=ceiling((sumb+sum)/d)));
bnd)

-- input: ({m1,...,mn},r,d,ea), n >=1 (number of points), m1, ... >=1 
-- (the multiplicities), r and d positive integers, ea any value
-- not bigger than the eventual value of bhalpha; can be set to 0 
-- output: Harbourne's unloading lower bound on alpha, using given r and d

bhalpha = (l,r,d,ea) -> (
i:=ea-1;
j:=0;
n:=0;
v:=prmt(zr(l));
tmpv:=v;
ttmpv:={};
scan(#l, j->(if v#j>0 then n=n+1));
if n>0 then (
   ww:={};
   scan(#l, j->(if j<r then ww=join(ww,{1}) else ww=join(ww,{0})));
   tmpi:=-1;
   while(tmpi < tmpv#0) do (
      i=i+1;
      tmpi=i;
      tmpv=v;
      while((tmpi*d-(dot(ww,tmpv))<0) and (tmpi >= tmpv#0)) do (
         tmpi=tmpi-d;
         ttmpv={};
         scan(#l, j->(ttmpv=join(ttmpv,{(tmpv#j)-(ww#j)})));
         tmpv=prmt(zr(ttmpv)))));
i)

-- Find bhalpha using best possible r and d;
-- ea is an a priori estimate for alpha (for speed)
-- it must be set to a value >= than the actual value 
-- of alpha (e.g., ea=findalpha(l))

bestbhalpha = (l,ea) -> (
w:={};
i:=0;
scan(#l, i->(if l#i>0 then w=join(w,{l#i})));
w=prmt(w);
n:=#w;
bnd:=0;
tmpbnd:=0;
r:=0;
d:=0;
tmpr:=0;
tmpd:=0;
if n>0 then (
   while(tmpr<n) do (tmpr=tmpr+1;
      tmpd=0;
      while(tmpd*tmpd<tmpr) do (tmpd=tmpd+1;
             tmpbnd=ea;
             while(tmpbnd==bhalpha(w,tmpr,tmpd,tmpbnd)) do tmpbnd=tmpbnd-1;
             tmpbnd=tmpbnd+1;
             if tmpbnd>bnd then (
                 bnd=tmpbnd;
                 r=tmpr;
                 d=tmpd))));
{bnd,r,d})

-- input: ({m1,...,mn},r,d,et), n >=1 (number of points), m1, ... >=1 (the multiplicities),
-- r and d positive integers, et any lower bound for tau (used for speed; can be set to 0).
-- output: Harbourne/Roe's algorithmic upper bound on tau, 
-- with given r and d (assumes char = 0).

HRtau = (l,r,d,et) -> (
i:=et-1;
j:=0;
n:=0;
v:=prmt(zr(l));
tmpv:=v;
ttmpv:={};
scan(#l, j->(if v#j>0 then n=n+1));
if n>0 then (
   ww:={};
   g:=(d-1)*(d-2)/2; -- genus of plane curve of degree d
   scan(#l, j->(if j<r then ww=join(ww,{1}) else ww=join(ww,{0})));
   tmpi:=0;
   while(tmpv#0 > 0) do (
      i=i+1;
      tmpi=i;
      tmpv=v;
      while((tmpi*d-(dot(ww,tmpv))>=g-1) and (tmpi>=d-2) and (tmpv#0 >0)) do (
         tmpi=tmpi-d;
         ttmpv={};
         scan(#l, j->(ttmpv=join(ttmpv,{(tmpv#j)-(ww#j)})));
         tmpv=prmt(zr(ttmpv)))));
i)

-- input: ({n,m},r,d,et), n >=1 (number of points), m >=1 (the uniform multiplicity),
-- r and d positive integers, et any lower bound for tau (used for speed; can be set to 0)
-- output: Harbourne/Roe's algorithmic upper bound on tau, using given r and d
-- (assumes char = 0).

unifHRtau = (l,r,d,et) -> (
i:=et-1;
n2:=0;  -- n2 is the number of points with maximum multiplicity
s:=0;
tmpm:=1;
tmpi:=0;
g:=(d-1)*(d-2)/2; -- genus of plane curve of degree d
while(tmpm > 0) do (
   i=i+1;
   tmpi=i;
   tmpm=l#1;
   n2= l#0;
   if r<=n2 then s=r*tmpm else s=r*tmpm-r+n2;
   while((tmpi*d-s>=g-1) and (tmpi>=d-2) and (tmpm >0)) do (
      tmpi=tmpi-d;
      if r<n2 then n2=n2-r else (
          n2=l#0-r+n2;
          tmpm=tmpm-1;
          if tmpm<0 then (
              tmpm=0;
              n2= l#0));
          if r<=n2 then s=r*tmpm else s=r*tmpm-r+n2));
i)

-- input: l={n,m}, n >=1 (number of points), m >=1 (the uniform multiplicity)
-- output: list {a,r,d}, where a is Harbourne/Roe's formulaic upper 
-- bound on tau (char 0) computed using r and d

ezunifHRtau = (l) -> (
n:=l#0;
m:=l#1;
d:=0;
while (d*d <= n) do d=d+1;
d=d-1;
r:=d;
while (r*r < d*d*n) do r=r+1;
g:=(d-2)*(d-1)/2;
a:= ceiling((m*r+g-1)/d);
b:=-2+d*ceiling(m*n/r);
if a<b then a=b;
{a,r,d})

-- input: l={n,m}, n >=1 (number of points), m >=1 (the uniform multiplicity)
-- output: the SHGH conjectured value of alpha; this is the actual value if 
-- n < 10, and an upper bound otherwise.

uniffindalpha = (l) -> (
n:=l#0;
m:=l#1;
a:=-1;
if n==1 then a=m;
if n==2 then a=m;
if n==3 then a=ceiling(3*m/2);
if n==4 then a=2*m;
if n==5 then a=2*m;
if n==6 then a=ceiling(12*m/5);
if n==7 then a=ceiling(21*m/8);
if n==8 then a=ceiling(48*m/17);
if n==9 then a=3*m;
if n>9 then (
  while(a*a-n*m*m+3*a-n*m+2 <0) do a=a+m;
  a=a-m;
  while(a*a-n*m*m+3*a-n*m+2 <=0) do a=a+1);
a)

-- input: l={n,m}, n >=1 (number of points), m >=1 (the uniform multiplicity)
-- output: the SHGH conjectured value of tau; this is the actual value if 
-- n < 10, and a lower bound otherwise.

uniffindtau = (l) -> (
n:=l#0;
m:=l#1;
t:=-1;
if n==1 then t=m-1;
if n==2 then t=2*m-1;
if n==3 then t=2*m-1;
if n==4 then t=2*m;
if n==5 then t=ceiling((5*m-1)/2);
if n==6 then t=ceiling((5*m-1)/2);
if n==7 then t=ceiling((8*m-1)/3);
if n==8 then t=ceiling((17*m-1)/6);
if n==9 then t=3*m;
if n>9 then (
  while(t*t-n*m*m+3*t-n*m+2 <0) do t=t+m;
  t=t-m;
  while(t*t-n*m*m+3*t-n*m+2 <0) do t=t+1);
if t<0 then t=0;
t)

-- input: l={n,m}, n >=1 (number of points), m (the multiplicity of each point)
-- output: Hirschowitz's lower bound for tau

Hiuniftau = (l) -> (
n:=l#0;
m:=l#1;
t:=m;
s:=n*m*(m+1);
a:=ceiling((t+3)/2);
b:=ceiling((t+2)/2);
while(a*b*2 <= s) do (
  t=t+1;
  a=ceiling((t+3)/2);
  b=ceiling((t+2)/2));
t)

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: Hirschowitz's lower bound for tau

Hitau = (l) -> (
n:=#l;
i:=0;
w:=prmt(zr(l));
t:=w#0;
s:=0;
scan(#l, i->(s=s+(w#i)*((w#i)+1)));
a:=ceiling((t+3)/2);
b:=ceiling((t+2)/2);
while(a*b*2 <= s) do (
  t=t+1;
  a=ceiling((t+3)/2);
  b=ceiling((t+2)/2));
t)

-- input: l={n,m}, n >=1 (number of points), m (the multiplicity of each point)
-- output: Gimigliano's lower bound for tau

Guniftau = (l) -> (
n:=l#0;
m:=l#1;
t:=0;
while(t*(t+3)<2*n) do t=t+1;
m*t)

-- input: l={m1,...,mn}, n >=1 (number of points), m1, ... >=1 (the multiplicities)
-- output: Gimigliano's lower bound for tau

Gtau = (l) -> (
n:=0;
w:=prmt(zr(l));
scan(#l, i->(if w#i >0 then n=n+1));
t:=0;
s:=0;
i:=0;
while(t*(t+3)<2*n) do t=t+1;
scan(#l, i->(if i<t then s=s+w#i));
s)

-- input: l={n,m}, n >=5 (number of points), m>0 (the multiplicity of each point)
-- output: Catalisano's lower bound for tau

Cuniftau = (l) -> (
s:=l#0;
m:=l#1;
r:=0;
t:=0;
f:=0;
while(f*(f+1) <= 2*s) do f=f+1;
f=f-1;
while(2*r<2*s-f*(f+1)) do r=r+1;
d1:=0;
d:=f;
if r==0 then d1=f-1 else d1=f;
t=d1+(m-1)*d;
if 2*t+1 < 5*m then t=ceiling((5*m-1)/2);
if t<2*m-1 then t=2*m-1;
if r == f then (if s >= 9 then t=m*d1+1);
t)

-- input: l={m1,...,mn}, n >=5 (number of points), m1, ... >=1 (the multiplicities)
-- output: Catalisano's lower bound for tau

Ctau = (l) -> (
n:=0;
i:=0;
w:=prmt(zr(l));
scan(#l, i->(if w#i >0 then n=n+1));
vm:={};
vs:={};
i=0;
while(i < n-1) do (
  if w#i > w#(i+1) then (
  vs=join(vs,{i+1});
  vm=join(vm,{w#i}));
  i=i+1);
vs=join(vs,{n});
vm=join(vm,{w#(n-1)});
i=#vm - 1;
v:={vm#i};
while(i > 0) do (
  v=join({vm#(i-1) - vm#i},v);
  i=i-1);
vf:={};
vr:={};
scan(#vs, i->(
  f:=0;
  r:=0;
  while(f*(f+1) <= 2*(vs#i)) do f=f+1;
  f=f-1;
  while(2*r<2*(vs#i)-f*(f+1)) do r=r+1;
  vf=join(vf,{f});
  vr=join(vr,{r})));
t:=0;
if (vr#(#vr-1)) == 0 then t = - 1;
d1:=t+vf#(#vf-1);
scan(#vs, i->(
  t=t+(vf#i)*(v#i)));
if 2*t+1 < (w#0)+(w#1)+(w#2)+(w#3)+(w#4) then 
t=ceiling(((w#0)+(w#1)+(w#2)+(w#3)+(w#4)-1)/2);
if t<(w#0)+(w#1)-1 then t=(w#0)+(w#1)-1;
if (vr#0) == (vf#0) then (if s >= 9 then (if (w#0)==(w#(n-1)) 
    then (if (w#0)>1 then t=(w#0)*d1+1)));
if (vr#0) == 0 then (if s > 9 then (if (w#0)==(w#(n-2)) 
    then (if (w#(n-1)) == 1 then t=(w#0)*d1+1)));
t)

-- gcd not implemented for polynomials in Macaulay 2, version 0.9.95
-- so I've copied a gcd definition from Macaulay 2, version 0.9.2

-- oldergcd(RingElement,RingElement) := RingElement => (r,s) -> (
oldergcd = (r,s) -> (
     if r == 0 then s
     else if s == 0 then r
     else (
          z := syz( matrix{{r,s}}, SyzygyLimit => 1 );
          a := z_(0,0);
          if s%a != 0 then error "can't find gcd in this ring";
          t := s // a;
--          if isField coefficientRing ring S then (
               c := leadCoefficient t;
               t = t // c;
--               );
          t))