# This program is called Res6point.
# File1=testschemes is a list of test schemes Z  
# Change the line below to "verbose = 0" to
# disable printing out all of the negative curves.
# Each line of testschemes should be of the form:
# 10 1 2 2 3 3 4
# The 10 represents which config out of the 11 possible
# to use. The rest of the numbers are the 6 multiplicities.

# Run this as:
#  awk -f Res6point testschemes

BEGIN {verbose = 1
# Make an array of all 6 point (-1)-classes
E[1]="0  1  0  0  0  0  0"
E[2]="0  0  1  0  0  0  0"
E[3]="0  0  0  1  0  0  0"
E[4]="0  0  0  0  1  0  0"
E[5]="0  0  0  0  0  1  0"
E[6]="0  0  0  0  0  0  1"
E[7]="1 -1 -1  0  0  0  0"
E[8]="1 -1  0 -1  0  0  0"
E[9]="1 -1  0  0 -1  0  0"
E[10]="1 -1  0  0  0 -1  0"
E[11]="1 -1  0  0  0  0 -1"
E[12]="1  0 -1 -1  0  0  0"
E[13]="1  0 -1  0 -1  0  0"
E[14]="1  0 -1  0  0 -1  0"
E[15]="1  0 -1  0  0  0 -1"
E[16]="1  0  0 -1 -1  0  0"
E[17]="1  0  0 -1  0 -1  0"
E[18]="1  0  0 -1  0  0 -1"
E[19]="1  0  0  0 -1 -1  0"
E[20]="1  0  0  0 -1  0 -1"
E[21]="1  0  0  0  0 -1 -1"
E[22]="2  0 -1 -1 -1 -1 -1"
E[23]="2 -1  0 -1 -1 -1 -1"
E[24]="2 -1 -1  0 -1 -1 -1"
E[25]="2 -1 -1 -1  0 -1 -1"
E[26]="2 -1 -1 -1 -1  0 -1"
E[27]="2 -1 -1 -1 -1 -1  0"
En=27

C[1]="0  0  0  0  0  0  0"

C[2]="1 -1 -1 -1  0  0  0"

C[3]="1 -1 -1 -1 -1  0  0"

C[4]="1 -1 -1 -1 -1 -1  0"

C[5]="1 -1 -1 -1 -1 -1 -1"

C[6]="1 -1 -1 -1 -1  0  0"
C[7]="1 -1  0  0  0 -1 -1"

C[8]="1 -1 -1 -1  0  0  0"
C[9]="1  0  0  0 -1 -1 -1"

C[10]="1 -1 -1 -1  0  0  0"
C[11]="1 -1  0  0 -1 -1  0"

C[12]="1 -1 -1 -1  0  0  0"
C[13]="1 -1  0  0 -1 -1  0"
C[14]="1  0 -1  0 -1  0 -1"

C[15]="1 -1 -1 -1  0  0  0"
C[16]="1 -1  0  0 -1 -1  0"
C[17]="1  0 -1  0 -1  0 -1"
C[18]="1  0  0 -1  0 -1 -1"

C[19]="2 -1 -1 -1 -1 -1 -1"

Start[1]=1
Finish[1]=1
Start[2]=2
Finish[2]=2
Start[3]=3
Finish[3]=3
Start[4]=4
Finish[4]=4
Start[5]=5
Finish[5]=5
Start[6]=6
Finish[6]=7
Start[7]=8
Finish[7]=9
Start[8]=10
Finish[8]=11
Start[9]=12
Finish[9]=14
Start[10]=15
Finish[10]=18
Start[11]=19
Finish[11]=19
}

{config=$1
testscheme=$2" "$3" "$4" "$5" "$6" "$7
print ".........................................."
print "Configuration: "config
print ".........................................."
print " "

# Find the minus 1 curves
k=0
for(i=1;i<=En;i++) {
  j=Start[config]
  test=0
  while(test==0 && j<=Finish[config]) {
    if(dot(E[i],C[j]) < 0)
      test=1
    j=j+1
    }
  if(test==0) {
    k=k+1
    M1[k]=i
    M1n=k
    }
  }

if(verbose==1)
{
print "The exceptional curves are:"
for(i=1;i<=M1n;i++)
 print "    "E[M1[i]]
print "The other negative curves are:"
for(i=Start[config];i<=Finish[config];i++)
 print "    "C[i]
}

print "Scheme: Z = "testscheme
Ftmp=0" "testscheme
flg=1
flag=0
while(flg==1) {
flg=0
for(i=1;i<=M1n;i++)
  if(dot("1 0 0 0 0 0 0",E[M1[i]])==0 && dot(Ftmp,E[M1[i]]) > 0) {
    Ftmp=subtract(Ftmp,mult(E[M1[i]],-1))
    flg=1
    flag=1}
for(i=Start[config];i<=Finish[config];i++)
  if(dot("1 0 0 0 0 0 0",C[i])==0 && dot(Ftmp,C[i]) > 0) {
    Ftmp=subtract(Ftmp,mult(C[i],-1))
    flg=1
    flag=1}
}
if(flag==1) {
split(Ftmp, D, " ")
testscheme=D[2]" "D[3]" "D[4]" "D[5]" "D[6]" "D[7]
print "Your Z does not satisfy the proximity inequalities."
print "Subtracting off components of blow ups to get a Z"
print "which does, turns your Z into: Z = "testscheme
}

h0=0
deg=-1
fflag=1
while(h0==0) {
deg=deg+1
testclass=deg" "mult(testscheme,-1)
split(h0andh1(testclass), hh, " ")
h0=hh[1]
h1=hh[2]
fflag=fixedflag(testclass)}
alpha=deg
print "alpha: "alpha

while(fflag==1) {
deg=deg+1
testclass=deg" "mult(testscheme,-1)
fflag=fixedflag(testclass)}
print "FCfree in degree: "deg
delta=deg+1

h1=1
deg=alpha-2
while(h1>0) {
deg=deg+1
testclass=deg" "mult(testscheme,-1)
split(h0andh1(testclass), hh, " ")
h0=hh[1]
h1=hh[2]}
tau=deg
print "tau: "tau


for(j=alpha;j<=delta+2;j++) {
testclass=j" "mult(testscheme,-1)
split(h0andh1(testclass), hh, " ")
hilb[j-alpha+1]=hh[1]
}


for(j=alpha;j<=delta+1;j++) {
testclass=j" "mult(testscheme,-1)
myclass=subtract(getFreepart(testclass),"-1 0 0 0 0 0 0")
split(h0andh1(myclass), hh, " ")
hilbfree[j-alpha+2]=hh[1]
}

gens[1]=hilb[1]

for(j=alpha+1;j<=delta+2;j++) {
gens[j-alpha+1]=hilb[j-alpha+1]-hilbfree[j-alpha+1]
if(hilbfree[j-alpha+1]>3*hilb[j-1-alpha+1])
  gens[j-alpha+1]=gens[j-alpha+1]+hilbfree[j-alpha+1]-3*hilb[j-1-alpha+1]
}

syz[1]=gens[1]-hilb[1]
syz[2]=gens[2]-hilb[2]+3*hilb[1]
syz[3]=gens[3]-hilb[3]+3*hilb[2]-3*hilb[1]

for(j=alpha+3;j<=delta+2;j++) {
syz[j-alpha+1]=gens[j-alpha+1]-hilb[j-alpha+1]+3*hilb[j-alpha]-3*hilb[j-alpha-1]+hilb[j-alpha-2]
}

print "deg    hilb_I  gens    syz    hilb_Z"
for(j=alpha;j<=delta+2;j++) 
print j"       "hilb[j-alpha+1]"       "gens[j-alpha+1]"      "syz[j-alpha+1]"       "(j+1)*(j+2)/2-hilb[j-alpha+1]

# end of testschemes loop
print ".........................................."
print " "

}


function fixedflag(F,    flg,i)    # return 1 if fixed part is nonempty, 0 otherwise
{flg=0
for(i=1;i<=M1n;i++)
  if(dot(F,E[M1[i]]) < 0) {
    flg=1}
for(i=Start[config];i<=Finish[config];i++)
  if(dot(F,C[i]) < 0) {
    flg=1}
return flg
}

function getFreepart(F, Ftmp,n,D,flg,i)
{Ftmp=F
n=split(F, D, " ")
flg=0
while(D[1]>=0 && flg==0) {
flg=1
for(i=1;i<=M1n;i++)
  if(dot(Ftmp,E[M1[i]]) < 0) {
    Ftmp=subtract(Ftmp,E[M1[i]])
    flg=0}
for(i=Start[config];i<=Finish[config];i++)
  if(dot(Ftmp,C[i]) < 0) {
    Ftmp=subtract(Ftmp,C[i])
    flg=0}
split(Ftmp, D, " ")
}
return Ftmp
}


# F must be a class for at least 3 blow ups
function h0andh1(F,     n, D, Ftmp, i, vdim, dim, flg)
{Ftmp=F
n=split(F, D, " ")
vdim=(D[1]+1)*(D[1]+2)
for(i=2;i<=n;i++)
  vdim=vdim-(-D[i]+1)*(-D[i])
i=int(vdim/2)-1
while(i*2<vdim)
  i=i+1
vdim=i # this just divided by 2

# Remove fixed components
flg=0
split(Ftmp, D, " ")
while(D[1]>=0 && flg==0) {
flg=1
for(i=1;i<=M1n;i++)
  if(dot(Ftmp,E[M1[i]]) < 0) {
    Ftmp=subtract(Ftmp,E[M1[i]])
    flg=0}
for(i=Start[config];i<=Finish[config];i++)
  if(dot(Ftmp,C[i]) < 0) {
    Ftmp=subtract(Ftmp,C[i])
    flg=0}
split(Ftmp, D, " ")
}

# The divisor has now been reduced 
if(D[1]<0)
  dim=0
else {
   dim=(D[1]+1)*(D[1]+2)
   for(i=2;i<=n;i++) {
     dim=dim-(-D[i]+1)*(-D[i])}
i=int(dim/2)-1
while(i*2<dim)
  i=i+1
dim=i # this just divided by 2
# h1=(dim-vdim)
# h0=dim
}
return dim" "(dim-vdim)
}


function dot(F,G,    n, x, i, m, y, s)
{$0=F
n=NF
for(i=1;i<=n;i++)
  x[i]=$i
$0=G
m=NF
for(i=1;i<=n;i++)
  y[i]=$i
s=x[1]*y[1]
for(i=2;i<=n;i++)
  s=s-x[i]*y[i]
return s
}

function subtract(F,G,  n,i,x,y,z)
{$0=F
n=NF
for(i=1;i<=n;i++)
  x[i]=$i
$0=G
for(i=1;i<=n;i++)
  y[i]=$i
z=x[1]-y[1]
for(i=2;i<=n;i++)
  z=z" "x[i]-y[i]
return z
}


function mult(F,c,  n,i,x,y,z)
{$0=F
n=NF
for(i=1;i<=n;i++)
  x[i]=c*$i
z=x[1]
for(i=2;i<=n;i++)
  z=z" "x[i]
return z
}