/****************************************************
This program asks for 7 nonnegative integer inputs, m_1, ..., m_7.
Given 7 general points p_1, ..., p_7 of P2, we define
the homogeneous ideal I (over an algebraically closed field k)
in the homogeneous coordinate ring R=k[P2] generated by all forms
vanishing to order at least m_i at p_i. This program then computes: 
various numerical characters for I; the hilbert function of I;
and the number v_d of generators in each degree d in a 
minimal homogeneous set of generators for I, and the number 
of generators in each degree for the syzygy module. (Note that 
the minimal free resolution of I over R has the form
0 -> F_1 -> F_0 -> I -> 0, where F_0 = \sum_d R[-d]^{v_d}, so
given F_0 and the hilbert function of I, one can determine F_1.)

The numerical characters printed out are: alpha (the least 
degree in which I is nonzero); beta (the least degree greater 
than or equal to alpha in which the forms in I of that degree 
do not have a common divisor of positive degree); and tau (the 
least degree such that the hilbert function of I first and forever 
after equals the hilbert polynomial of I). In addition, the 
Campanella (see [J. Alg, 1986]) upper and lower bounds for the 
minimal number of generators of I is also printed out, 
together with the actual total number of generators.
Then the values of the hilbert function (i.e., the dimension of 
the homogeneous component I_d of I of degree d) for degrees d 
from alpha to tau+1 are printed out, followed for each degree by 
an expression of the form lb<= g1 + g2 (maxrkValue)<= ub,
where lb (ub) is Campanella's lower (upper) bound for the number 
of generators in that degree, g1 is the number of 
generators in that degree coming from the loss of fixed components 
as the degree increases, and g1 + g2 is the actual
number of generators in that degree.

For degree d with alpha < d <= tau +1, fcgens in degree d turns 
out to be (h^0(X, F_d)-h^0(X,e_0+(F_{d-1})_f)), where X is 
the plane blown up at the n general points, e_0 is the divisor 
class of a line on X, e_i for i>0 are the divisor classes of the 
exceptional loci of the blown up points, F_d=de_0-m_1e_1-\cdots-m_1e_n) 
is the divisor class corresponding to I_d, and (F_d)_f is the class of the
mobile (or free) part of |F_d|. Also, neffgens in degree d is the 
dimension of the cokernel of 
H^0(X,(F_{d-1})_f)\otimes H^0(X, e_0)\to H^0(X, e_0+(F_{d-1})_f),
while maxrkValue is h^0(X,(F_{d-1})_f) minus the dimension of  
H^0(X, e_0)\to H^0(X, e_0+(F_{d-1})_f); thus the maximum of
maxrkValue and 0 is the dimension the cokernel would have 
assuming this map had maximal rank. See my papers (B. Harbourne) 
on resolutions of fat point subschemes for details.

Finally, the number of generators in each degree for 
F_0 and F_1 are also printed out.

March 1, 1998 B. Harbourne   bharbourne@unl.edu
************************************************/


#include "stdio.h"
#include "time.h"
/* If you increase MAXMULT you may need to increase the MAXINDEX values. */
#define MAXMULT 300 
#define MAXINDEX1  320 

/* PICNMBR is 1 plus the number of points; it is the Picard number
of the surface obtained by blowing up the points. The number of points should 
be 7. If you change PICNMBR, you must remove the "7 point specific code".
Search for the string "7 point specific code" to find where the code
which must be removed starts and stops. With that removed, the program
should still work for values 1<=PICNMBR<=6. For 7<=PICNMBR with that code removed,
the program will compute numbers of generators and syzygies ASSUMING max rank holds
for H^0(X, F)\otimes H^0(X, L)\to H^0(X, F+L), where L is the class of a line
and F is effective and numerically effective on the blow up X of P2 at the points.
This assumption is however false in general. Moreover, for 10<=PICNMBR,
the program's values for the hilbert function are based on the conjecture 
that H^1(X, F)=0 for a neff and effective F when X is a general blow up
of P2, and that the only integral curves on X of negative self-intersection
are exceptional curves. */

#define PICNMBR 8 

void elimfcs( int class2[]);
void mypermute( int class2[],  int class3[]);
void onestep( int class2[],  int class3[]);
void gotofd( int class2[],  int class3[]);
int FindAlpha( int mult[]);
long int findExpectedH( int class2[]);
long int combin( int d);

main()
{  int mult[PICNMBR], class2[PICNMBR], class3[PICNMBR], Line[PICNMBR];
   int alpha, Galpha[PICNMBR], newL[PICNMBR];
   long int h[MAXINDEX1], conds, lambda, t, h2[MAXINDEX1];
   long int dh[MAXINDEX1], ddh[MAXINDEX1], dddh[MAXINDEX1];
   long int eps[MAXINDEX1], ub[MAXINDEX1], lb[MAXINDEX1];
   long int neffnu[MAXINDEX1], fcnu[MAXINDEX1], sumL, nu[MAXINDEX1];
   int w, x, y, z, beta, tau, F[PICNMBR], a, sum;
   int F0[PICNMBR], F1[PICNMBR], F2[PICNMBR], F3[PICNMBR];
   int F4[PICNMBR], F5[PICNMBR], F6[PICNMBR], F7[PICNMBR]; 
   int q, i, flag, j, k, stopit;
   long int testflag, numsyzgen[MAXINDEX1+1], f1, f2;

printf("This program will ask you for 7 multiplicities.\n");
printf("It then computes a resolution of the corresponding\n");
printf("ideal of fat points. \n");
printf("**** ENTER 1 TO CONTINUE, 0 TO EXIT ********\n");
stopit=0;

if((q=scanf( " %d", &i))==1 && i==1)
{    
while(stopit==0)
{    
/* initialization of class of a line*/
for(i=1; i<PICNMBR; i+=1) 
    Line[i]=0;
Line[0]=1;

/* get the multiplicities */
for(i=0; i<PICNMBR; i+=1)   
mult[i]=0;

printf("Enter the multiplicities of your %d points; i.e.,\n", PICNMBR-1);
printf("enter %d integers, separated by spaces.\n", PICNMBR-1);
printf("(Your inputs will be constrained to be between 0 and %d.)\n", MAXMULT);
for(i=1; i<PICNMBR; i+=1)   
    q=scanf( " %d", &mult[i]); 
for(i=0; i<PICNMBR; i+=1)   
   if(mult[i]>MAXMULT)
        mult[i]=MAXMULT;


/* remove negative multiplicities */  
elimfcs(mult);

/* compute expected number of conditions imposed by those multiplicities */
conds=-findExpectedH(mult)+1L;


/* Find alpha: least deg in which ideal is nonzero */
alpha= FindAlpha(mult);


/* Find translate Galpha of Falpha in fund dom
     and translate newL of Line under same transf  */
     /* initialize class2 and class3 */
for(i=1; i<PICNMBR; i+=1)   
      class2[i]=mult[i];
class2[0]=alpha;
for(i=0; i<PICNMBR; i+=1)   
      class3[i]=Line[i];
gotofd(class2, class3);
for(i=0; i<PICNMBR; i+=1)   
      {Galpha[i]=class2[i];
       newL[i]=class3[i];}


/* Find beta: least degree with no fixed components in that degree */
z=0;
w=0;
for (i=1; i<PICNMBR; i+=1)
  {
   x=Galpha[i];
   y=newL[i];
   if(x<0)
      {
       x=-x;
       w=x/y;
       if(x>w*y)
          w=w+1;
      }
   if(w>z)
      z=w;
  }
beta=alpha+z;

/* Find tau: the least degree in which the conditions are indep */
if(alpha==beta)
{
     tau=beta;
     if(combin(alpha+1) == conds)
     tau=alpha - 1;
}
else
{    z=beta - alpha - 1;
     flag=0;
     tau=beta - 1;
     for(i=1; i<PICNMBR; i+=1)
         {
          x=Galpha[i];
          y=newL[i];
          if (x<0)
          {   x=-x;
              if(x*1L > (z*y)*1L+1L)
                  flag=1;
          }
      }
    if(flag == 1)
        tau=beta;
}


printf("\n");
/**********************START 7 point specific code******************/
printf("\nGiven 7 general points of the projective plane\n");
printf("taken with multiplicities %d, %d, %d, %d, %d, %d, %d,\n", mult[1],mult[2],mult[3],mult[4],mult[5],mult[6],mult[7]);
printf("a resolution of the ideal I(%dp_1+%dp_2+%dp_3+%dp_4+%dp_5+%dp_6+%dp_7)\n", mult[1],mult[2],mult[3],mult[4],mult[5],mult[6],mult[7]);
/**********************END 7 point specific code******************/
printf("over the homogeneous coordinate ring R=k[x,y,z] of P2 is given.\n\n");

printf("The resolution takes the form 0 -> F_1 -> F_0 -> I -> 0;\n");
printf("among the results returned are the number of homogeneous generators\n");
printf("for each degree in which F_0 has generators; likewise for F_1.\n");
printf("Also returned are:\n");
printf("     alpha: the least degree in which I has a nonzero element.\n");
printf("     beta : the least degree d for which the elements of the \n");
printf("            homogeneous component I_d has no nontrivial common divisor.\n");
printf("     tau  : the least degree d for which the fat points impose\n");
printf("            independent conditions in degree d.\n");
printf("In addition, the value dim I_d of the Hilbert function of I in\n");
printf("each degree d from alpha to tau+1 is given, together with:\n");
printf("Campanella's lower bounds lb; the number g1 of generators forced \n");
printf("by the occurrence of nontrivial common divisors in I_d; the number\n");
printf("g2 of generators remaining in degree d; then maxrkValue (which\n");
printf("gives the value g2 would have if a certain maximal rank property held);\n");
printf("and finally Campanella's upper bounds ub.\n");

printf("\n");
printf("Ideal characters:    alpha = %d, beta = %d, tau = %d.\n", alpha, beta, tau);
printf("Campanella/Dubreil upper bound for the number of gens = alpha+beta-tau = %d.\n", alpha+beta-tau);

/* MAXINDEX1 need not be bigger than tau-alpha+1, but to be safe keep it even smaller */
if(tau-alpha>MAXINDEX1)
printf("Multiplicities too large; array bounds exceeded. EXIT EXIT EXIT\n");
else
{

/* Find the hilbert function */
h[0]=0;
for(i=1; i<PICNMBR; i+=1)   
      class2[i]=Galpha[i];
elimfcs(class2);
h[1]=findExpectedH(class2);
for(i=alpha+1; i<tau+2; i=i+1)   
{    x=newL[0];
     a=Galpha[0];
     F[0]=a + (i-alpha)*x;
     for(j=1; j<PICNMBR; j+=1)   
         F[j]=(i-alpha)*newL[j]+Galpha[j];
     elimfcs(F);
     h[i-alpha+1]=findExpectedH(F);
}

/* Here find fcnu; note h2=h^0(Freepart(Galpha+i*newL) + newL)  */
h2[1]=0;
for(i=alpha; i<tau+2; i=i+1)   
{    x=newL[0];
     a=Galpha[0];
     F[0]=a + (i-alpha+1)*x;
     for(j=1; j<PICNMBR; j+=1)   
         F[j]=(i-alpha)*newL[j]+Galpha[j];
     elimfcs(F);
     for(j=1; j<PICNMBR; j+=1)   
         F[j]=F[j]+newL[j];
     h2[i-alpha+2]=findExpectedH(F);
}
fcnu[1]=h[1];
for(i=alpha+1; i<tau+2; i+=1)
    fcnu[i-alpha + 1]=h[i-alpha + 1]-h2[i-alpha + 1];

/* Find neffnu's assuming max rank*/
neffnu[1]=0L;
for(i=alpha; i<tau+1; i+=1)
{   lambda= h2[i - alpha + 2]-3*h[i - alpha + 1];
    if(lambda < 0L)
	    lambda=0L;
    neffnu[i - alpha + 2]= lambda;
}

/**********************START 7 point specific code******************/

/* Here check for the exceptional cases in which
t>lambda but max rank still holds. This can only
happen in degree alpha+1, and only then when
tau>=alpha. This is because there are no gens
in degrees beyond tau+1, and because the exceptional
cases only occur for divisor classes F with
F-e_0 not effective.*/

flag=0;
if(tau>=alpha)
{   x=newL[0];
    for(i=0; i<PICNMBR; i+=1)
    F[i]=Galpha[i];
    elimfcs(F);
    for(j=0; j<PICNMBR; j+=1)
    class2[j]=newL[j];
    gotofd(class2, F);
	/* Order the multiplicities of F */
    for(k=1; k<PICNMBR-1; k+=1)
    for(j=k+1; j<PICNMBR; j+=1)
	if(F[j]>F[k])
	   {y=F[j];
        F[j]=F[k];
        F[k]=y;}        
    F0[0]=0;
    F0[1]=0;
    F0[2]=0;
    F0[3]=0;
    F0[4]=0;
    F0[5]=0;
    F0[6]=0;
    F0[7]=0;
    F1[0]=4;
    F1[1]=2;
    F1[2]=2;
    F1[3]=2;
    F1[4]=1;
    F1[5]=1;
    F1[6]=1;
    F1[7]=1;
    F2[0]=7;
    F2[1]=3;
    F2[2]=3;
    F2[3]=3;
    F2[4]=3;
    F2[5]=2;
    F2[6]=2;
    F2[7]=2;
    F3[0]=10;
    F3[1]=4;
    F3[2]=4;
    F3[3]=4;
    F3[4]=4;
    F3[5]=4;
    F3[6]=3;
    F3[7]=3;
    F4[0]=13;
    F4[1]=5;
    F4[2]=5;
    F4[3]=5;
    F4[4]=5;
    F4[5]=5;
    F4[6]=5;
    F4[7]=4;
    F5[0]=16;
    F5[1]=6;
    F5[2]=6;
    F5[3]=6;
    F5[4]=6;
    F5[5]=6;
    F5[6]=6;
    F5[7]=6;
    F6[0]=5;
    F6[1]=2;
    F6[2]=2;
    F6[3]=2;
    F6[4]=2;
    F6[5]=2;
    F6[6]=2;
    F6[7]=1;
    F7[0]=8;
    F7[1]=3;
    F7[2]=3;
    F7[3]=3;
    F7[4]=3;
    F7[5]=3;
    F7[6]=3;
    F7[7]=3;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F0[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F1[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F2[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F3[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F4[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F5[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F6[j])
		      flag=flag+1;
    if(flag<8)
       flag=0;
    if(flag<8)
	  for(j=0; j<8; j+=1)
	    if(F[j] == F7[j])
		      flag=flag+1;
}


/* Adjust neffnu's to account for occurrence of nonmax ranks */
for(i=alpha; i<tau+1; i+=1)
{   x=newL[0];
    a=Galpha[0];
    F[0]= a + (i-alpha)*x;
    for(j=1; j<8; j+=1)
       F[j]=(i-alpha)*newL[j]+Galpha[j];
    elimfcs(F);
    for(j=0; j<8; j+=1)
       class2[j]=newL[j];
    gotofd(class2, F);
    sum=3*F[0];
    for(j=1; j<8; j+=1)
       sum=sum - F[j];
    t=0L;
    for(j=1; j<8; j+=1)
        if(sum == F[j]) 
		  t=t + 1L;
   
    if(i==alpha)
      if(flag<8)
        if(t>neffnu[i - alpha + 2])
           neffnu[i - alpha + 2]= t;
    if(i>alpha)
        if(t>neffnu[i - alpha + 2])
           neffnu[i - alpha + 2]= t;
}

/**********************END 7 point specific code******************/

sumL=0L;
for(i=alpha; i<tau+2; i+=1)
{     nu[i - alpha+1]=fcnu[i - alpha + 1]+neffnu[i - alpha+1];
      sumL=sumL+nu[i - alpha+1];
}	 


/* Here we determine the Campanella Bounds */
/* Note that the indexing for most of the arrays has i=1
correspond to degree alpha, but for numsyzgen[],
dh[], ddh[], dddh[], ubl[] and l[], which I borrowed
from a program I wrote earlier, the indexing has i=0
correspond to alpha. */

    dh[0]=h[1];
    ddh[0]=dh[0];
    dddh[0]=ddh[0];
    for(j=1; j<= tau+1-alpha; j+=1)
      {dh[j]=h[j+1]-h[j];
       ddh[j]=dh[j]-dh[j-1];
       dddh[j]=ddh[j]-ddh[j-1];  
       }  


    for(j=0; j<= beta-1-alpha; j+=1)
       eps[j]=1L;

    for(j=beta-alpha; j<= tau+1-alpha; j+=1)
       eps[j]=2L;

         ub[0]=ddh[0];
         lb[0]=dddh[0];

    for(j=1; j<= tau+1-alpha; j+=1)
       { ub[j]=ddh[j]- eps[j-1];
         lb[j]=dddh[j];
        if(lb[j]<0) 
           lb[j]=0;
        if(ub[j]<0) 
           ub[j]=0;
        }
	
	if(lb[beta-alpha]<1L)
	     lb[beta-alpha]=1L;  

t=0L;
for(j=0; j<= tau+1-alpha; j+=1)
t=t+lb[j];

printf("Campanella lower bound for the number of gens = %ld.\n", t);
printf("Actual total number of gens =  %ld.\n", sumL);
printf("\n");
printf("dim I_d:       lb <= g1 + g2 (maxrkValue) <= ub\n");
for(i=alpha; i<tau+2; i+=1)
    printf("dim I_%d = %ld: %ld <= %ld + %ld (maxrkValue= %ld) <= %ld\n", i, h[i-alpha+1], lb[i-alpha], fcnu[i-alpha+1], neffnu[i-alpha+1], h2[i - alpha + 1]-3*h[i - alpha], ub[i-alpha]);
printf("\n");
for(i=alpha; i<tau+2; i+=1)
      printf("Number of gens of F_0 in degree  %d is %ld\n", i, nu[i-alpha+1]);

j=0;
testflag=0L;
for(k=alpha;k<alpha+MAXINDEX1+1;k+=1)
numsyzgen[k-alpha]=0L;

while(testflag<sumL-1L)
     {f1=0L;
      f2=0L;
      j+=1;

      if(j<tau+2-alpha)
         {for(k=0;k<=j;k+=1)
             f1=f1+nu[k+1]*combin(j-k+2);
          for(k=0;k<=j;k+=1)
             f2=f2+numsyzgen[k]*combin(j-k+2);
          numsyzgen[j]=f1-f2-h[j+1];
          if(numsyzgen[j] != 0)
            testflag=testflag+numsyzgen[j];
          }
       
      if(j>tau+1-alpha)
         {for(k=0;k<=tau+1-alpha;k+=1)
             f1=f1+nu[k+1]*combin(j-k+2);
          for(k=0;k<=j;k+=1)
             f2=f2+numsyzgen[k]*combin(j-k+2);
          numsyzgen[j]=f1-f2- (combin(j+alpha+2)-conds);
          if(numsyzgen[j] != 0)
            testflag=testflag+numsyzgen[j];
          }
     }
 /*    omega=j;  */
              
     printf("\n");

     for(k=1;k<=j;k+=1)
        if(numsyzgen[k]>0)
            printf("Number of gens of F_1 in degree %d is %ld\n", k+alpha, numsyzgen[k]);

} /* end if checking for exceeding of array bounds */

printf("**** ENTER 1 TO RUN AGAIN, 0 TO EXIT ********\n");


if(!((q=scanf( " %d", &i))==1 && i==1))
stopit=1;

}/* end of program while loop  */

} /* end of if checking user's intention to run program */
else
printf("User has terminated program run.\n");
}  /* end of main */


void elimfcs( int class2[])
{ int i;
for(i=1;i<PICNMBR; i+=1)
if(class2[i]<0)
class2[i]=0;
}

void mypermute( int class2[],  int class3[])
{  int max, max3, test, test3, i;

   max=class2[1];
   max3=class3[1];
   for(i=2; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test>max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[1]=max;
      class3[1]=max3;
    }

   max=class2[2];
   max3=class3[2];
   for(i=3; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test>max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[2]=max;
      class3[2]=max3;
    }
  
   max=class2[3];
   max3=class3[3];
   for(i=4; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test>max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[3]=max;
      class3[3]=max3;
    }
  
}


void onestep( int class2[],  int class3[]) 
{   int a, b, c, d; 
    a=class2[0];
    b=class2[1];
    c=class2[2];
    d=class2[3];
    if(a-b-c-d<0)
    {   class2[0]=2*a-b-c-d;
        class2[1]=a-c-d;
        class2[2]=a-b-d;
        class2[3]=a-b-c;

        a=class3[0];
        b=class3[1];
        c=class3[2];
        d=class3[3];
        class3[0]=2*a-b-c-d;
        class3[1]=a-c-d;
        class3[2]=a-b-d;
        class3[3]=a-b-c;
    }
}

void gotofd( int class2[],  int class3[])
{  int effflag;
   int beforedegree, afterdegree;
   effflag=0;
   while(effflag==0)
    {
       beforedegree=class2[0];
       mypermute(class2,class3);
       onestep(class2,class3);
       afterdegree=class2[0];
       if(afterdegree==beforedegree)
           effflag=1;
    }
}


int FindAlpha( int mult[])
{   int i, effflag;
    int alpha, beforedegree, afterdegree, temp[PICNMBR], dummy[PICNMBR], class2[PICNMBR];

  for(i=1; i<PICNMBR; i+=1)
      class2[i]= mult[i];
  mypermute(class2, dummy);
  alpha=class2[1];
  class2[0]=alpha;

  for(i=0; i<PICNMBR; i+=1)
     temp[i]=class2[i];

  effflag=0;
  while (effflag <1)
  {  while (effflag==0)
     {
        beforedegree=class2[0];
        mypermute(class2, dummy);
        elimfcs(class2);
        onestep(class2, dummy);
        afterdegree=class2[0];
        if (afterdegree<0)
            effflag=-1;
        if (afterdegree==beforedegree)
            effflag=1;
    }
    if (effflag<1)
    { effflag=0;
      alpha=alpha + 1;
      for(i=0; i<PICNMBR; i+=1)
         class2[i]=temp[i];
      class2[0]=alpha;
      for(i=0; i<PICNMBR; i+=1)
         temp[i]=class2[i];
    }
  }
  return alpha;
}

long int findExpectedH( int class2[])
{   int i;
    long int sum, x;
    x=class2[0];
    sum=x*x+3*x;
    for(i=1; i<PICNMBR; i+=1)
        {x=class2[i];
        sum=sum-x*x-x;}
    sum=sum/2L + 1L;
    return sum;
}


long int combin(int d)
{ long int x, y;
    y=d;
    if(d<2)
       x=0L;
    else
       x=y*(y-1L)/2L;
    return x;
}