/****************************************************
This program asks for 8 nonnegative integers, m_1, ..., m_8
from standard input.
Given 8 general points p_1, ..., p_8 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
syzygies in each degree. The output is to standard output,
in html.

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 bound for the
minimal number of generators of I is also printed out.


Sep 30, 1998 revised Aug 8, 1999 B. Harbourne   bharbourne1@unl.edu
************************************************/

#include "time.h"
#include "stdio.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 8. If you change PICNMBR, you must remove the "8 point specific code".
Search for the string "8 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 8<=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 9

void nuelimfcs( int class2[]);
void numypermute( int class2[],  int class3[]);
void nuonestep( int class2[],  int class3[]);
void nugotofd( int class2[],  int class3[]);
long int nufindExpectedH( int class2[]);
long int nucombin( int d);
long int nufindH0( int class2[]);
long int nufindH1( int class2[]);
long int nulstar( int class2[]);
long int nuqstar( int class2[]);
long int numxrkcok( int class2[]);
void numakemonotone( int class2[]);
int nuesix( int class2[]);
int nuefive( int class2[]);
int nuefour( int class2[]);
int nuethree( int class2[]);
int nuetwo( int class2[]);
int nueone( int class2[]);
void nusix( int class2[]);
void nufive( int class2[]);
void nufour( int class2[]);
void nuthree( int class2[]);
void nutwo( int class2[]);
void nuone( int class2[]);
void nufreepart( int class2[]);
int nuonesteptoample( int class2[]);
int numinusdotK( int class2[]) ; 
int nuabs( int x) ; 
long int nucok( int class2[]);

main()
{  int mult[PICNMBR], stopit, class4[PICNMBR], nisha, kira;
   int alpha, beta, tau;
   long int h[MAXINDEX1], conds;
   long int sumL, newnu[MAXINDEX1], mxrkvalue;
   int q, i, j, k, ii;
   long int testflag, numsyzgen[MAXINDEX1+1], f1, f2;	


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

q=scanf( "%d %d %d %d %d %d %d %d", &mult[0],&mult[1],&mult[2],&mult[3],&mult[4],&mult[5],&mult[6],&mult[7]);   


/* remove negative multiplicities */
for(i=0; i<PICNMBR; i+=1)   
	if(mult[i]<0)
		mult[i]=0;

alpha=0;
class4[0]=0;
for(i=0; i<PICNMBR; i+=1)   
	class4[i+1]=-mult[i];
	
/* Find alpha: least deg in which ideal is nonzero */
while(nufindH0(class4)==0) 
	class4[0]++;

alpha=class4[0];

class4[0]=0;

/* Find tau: least deg >= alpha in which h1 is zero */
while(nufindH1(class4)>0)   
	class4[0]++;
tau=class4[0];

/* Find beta: least degree with no fixed components in that degree */
class4[0]=alpha;
i=alpha;
nufreepart(class4);
j=class4[0];
while(j<i++)   
	{class4[0]=i;
	 for(k=0; k<PICNMBR; k+=1)   
	 	class4[k+1]=-mult[k];
	 nufreepart(class4);
	 j=class4[0];}
beta=j;


printf("Content-type: text/html%c%c",10,10);
printf("<html>\n<head>\n<title>Generators for Ideals of 8 Fat General Points on P2</title>\n</head>");
printf("<body>");
printf("<pre>");

/**********************START 8 point specific code******************/
printf("\nGiven 8 general points of the projective plane\n");
printf("taken with multiplicities %d, %d, %d, %d, %d, %d, %d, %d\n", mult[0],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+%dp_8)\n", mult[0],mult[1],mult[2],mult[3],mult[4],mult[5],mult[6],mult[7]);
/**********************END 8 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 bounds on the number of \n");
printf("homogeneous generators for each degree in which F_0 has generators.\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("the actual number of generators in that degree.)\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
{

for(i=alpha; i<tau+2; i+=1)
{     class4[0]=i-1;
		for(k=0; k<PICNMBR; k+=1)   
			class4[k+1]=-mult[k];
      newnu[i-alpha]=nucok(class4);
}

for(i=0; i<PICNMBR; i+=1)   
	 class4[i+1]=-mult[i];

for(i=alpha; i<tau+2; i+=1)
{     class4[0]=i;
      h[i-alpha+1]=nufindH0(class4);
}


/**********************END 8 point specific code******************/

sumL=0L;
for(i=alpha; i<tau+2; i+=1)
      sumL=sumL+newnu[i-alpha];

printf("<br>\n");
printf("The actual total number of generators is: %ld\n", sumL);
printf("<br>\n");

printf("[Note:The max rank value below gives the number of generators, under the assumption\n");
printf("that the rank of the multiplication map is maximal. (A negative value is\n");
printf("just minus the dimension of the kernel, so under the maximal rank assumption\n");
printf("this means that the map is onto and there are no generators; on the other hand\n");
printf("a nonnegative value means the map is injective, with the nonnegative value\n");
printf("then just given the number of generators under the maximal rank assumption.)]\n");
printf("<br>\n");
printf("<br>\n");

for(i=alpha; i<tau+2; i+=1)
{for(j=0; j<PICNMBR; j+=1)   
	class4[j+1]=-mult[j];
	class4[0]=i-1;
class4[1]=class4[1]-1;  
nisha=nufindH1(class4);
for(j=0; j<PICNMBR; j+=1)   
	class4[j+1]=-mult[j];
	class4[0]=i-2;
class4[1]=class4[1]+1;
kira=nufindH1(class4);

      printf("dim I_%d = %ld   number of gens =%ld  (max rank value would be: %ld, q* = %ld, l* = %ld )\n", i, h[i-alpha+1], newnu[i-alpha], h[i-alpha+1]-3*h[i-alpha], nisha, kira);}

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

class4[0]=0;
for(i=0; i<PICNMBR; i+=1)   
	class4[i+1]=-mult[i];
conds=1L+nufindH1(class4);

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

      if(j<tau+2-alpha)
         {for(k=0;k<=j;k+=1)
             f1=f1+newnu[k]*nucombin(j-k+2);
          for(k=0;k<=j;k+=1)
             f2=f2+numsyzgen[k]*nucombin(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+newnu[k]*nucombin(j-k+2);
          for(k=0;k<=j;k+=1)
             f2=f2+numsyzgen[k]*nucombin(j-k+2);
          numsyzgen[j]=f1-f2- (nucombin(j+alpha+2)-conds);
          if(numsyzgen[j] != 0)
            testflag=testflag+numsyzgen[j];
          }
     }

 /*  Note that 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 */
}  /* end of main */




long int nucok( int class2[])
{  int class3[PICNMBR];
   int i, testflag, j, s, f;
   long int cok, ker;
	
for(i=0; i<PICNMBR; i+=1)   
	class3[i]=class2[i];
cok=numxrkcok(class2);
ker=0L;
if(cok<0L)
	{ker=-cok;
	cok=0L;}
	
testflag=0;

if(nufindH0(class2)<=1)
	testflag=1;
while(testflag==0)
{nufreepart(class2);
if(nuonesteptoample(class2)==1)
	{if(nueone(class2)==0)
		testflag=2;
	else
		{if(class2[8]==0)
			testflag=1;
		else
			{if(numinusdotK(class2)==1)
				testflag=1;
			else
				{if(nuesix(class2)==1)
					nusix(class2);
				else
					{if(nuefive(class2)==1)
						nufive(class2);
					else
						testflag=3;
					}
				}
			}
		}
	}
if(nufindH0(class2)<=1)
		testflag=1;
}

numakemonotone(class2);
if(testflag==1)
	if(numxrkcok(class2)<0)
		cok=cok-(numxrkcok(class2)+ker);

if(testflag==2)
	cok=cok+(nuqstar(class2)+nulstar(class2)-numxrkcok(class2))-ker;

i=class2[0];
j=-1;
if(testflag==3)
	{while(i>=0)
		{i=i-8;
		j=j+1;}
	if((nuabs(class2[0]-j*8-3)+nuabs(class2[1]+j*3+1)+nuabs(class2[2]+j*3+1)+nuabs(class2[3]+j*3+1)+nuabs(class2[4]+j*3+1)+nuabs(class2[5]+j*3+1)+nuabs(class2[6]+j*3+1)+nuabs(class2[7]+j*3+1)+nuabs(class2[8]+j))==0)
		cok=cok+j+1-ker;
	else
	{
	if(numxrkcok(class2)<0)
		cok=cok-(numxrkcok(class2)+ker);
	}
	}
return cok;
}

void nusix( int class2[])
{   class2[0]=class2[0]-6;
	class2[1]=class2[1]+3;
	class2[2]=class2[2]+2;
	class2[3]=class2[3]+2;
	class2[4]=class2[4]+2;
	class2[5]=class2[5]+2;
	class2[6]=class2[6]+2;
	class2[7]=class2[7]+2;
	class2[8]=class2[8]+2;
}

void nufive( int class2[])
{   class2[0]=class2[0]-5;
	class2[1]=class2[1]+2;
	class2[2]=class2[2]+2;
	class2[3]=class2[3]+2;
	class2[4]=class2[4]+2;
	class2[5]=class2[5]+2;
	class2[6]=class2[6]+2;
	class2[7]=class2[7]+1;
	class2[8]=class2[8]+1;
}

void nufour( int class2[])
{   class2[0]=class2[0]-4;
	class2[1]=class2[1]+2;
	class2[2]=class2[2]+2;
	class2[3]=class2[3]+2;
	class2[4]=class2[4]+1;
	class2[5]=class2[5]+1;
	class2[6]=class2[6]+1;
	class2[7]=class2[7]+1;
	class2[8]=class2[8]+1;
}

void nuthree( int class2[])
{   class2[0]=class2[0]-3;
	class2[1]=class2[1]+2;
	class2[2]=class2[2]+1;
	class2[3]=class2[3]+1;
	class2[4]=class2[4]+1;
	class2[5]=class2[5]+1;
	class2[6]=class2[6]+1;
	class2[7]=class2[7]+1;
}

void nutwo( int class2[])
{   class2[0]=class2[0]-2;
	class2[1]=class2[1]+1;
	class2[2]=class2[2]+1;
	class2[3]=class2[3]+1;
	class2[4]=class2[4]+1;
	class2[5]=class2[5]+1;
}

void nuone( int class2[])
{   class2[0]=class2[0]-1;
	class2[1]=class2[1]+1;
	class2[2]=class2[2]+1;
}


int nuesix( int class2[])
{  return class2[0]*6+class2[1]*3+class2[2]*2+class2[3]*2+class2[4]*2+class2[5]*2+class2[6]*2+class2[7]*2+class2[8]*2;
}

int nuefive( int class2[])
{  return class2[0]*5+class2[1]*2+class2[2]*2+class2[3]*2+class2[4]*2+class2[5]*2+class2[6]*2+class2[7]+class2[8];
}

int nuefour( int class2[])
{  return class2[0]*4+class2[1]*2+class2[2]*2+class2[3]*2+class2[4]+class2[5]+class2[6]+class2[7]+class2[8];
}

int nuethree( int class2[])
{  return class2[0]*3+class2[1]*2+class2[2]+class2[3]+class2[4]+class2[5]+class2[6]+class2[7];
}

int nuetwo( int class2[])
{  return class2[0]*2+class2[1]+class2[2]+class2[3]+class2[4]+class2[5];
}

int nueone( int class2[])
{  return class2[0]+class2[1]+class2[2];
}

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

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

void numypermute(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;
    }
  
   max=class2[4];
   max3=class3[4];
   for(i=5; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test<max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[4]=max;
      class3[4]=max3;
    }
  
   max=class2[5];
   max3=class3[5];
   for(i=6; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test<max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[5]=max;
      class3[5]=max3;
    }
  
   max=class2[6];
   max3=class3[6];
   for(i=7; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test<max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[6]=max;
      class3[6]=max3;
    }
  
   max=class2[7];
   max3=class3[7];
   for(i=8; i<PICNMBR; i+=1)
   {  test=class2[i];
      test3=class3[i];
      if(test<max)
      {  class2[i]=max;
         max=test;
         class3[i]=max3;
         max3=test3;
       }
      class2[7]=max;
      class3[7]=max3;
    }
  
}


void nuonestep(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 nugotofd(int class2[], int class3[])
{  int effflag;
   int beforedegree, afterdegree;
   effflag=0;
   while(effflag==0)
    {
       beforedegree=class2[0];
       numypermute(class2,class3);
       nuonestep(class2,class3);
       afterdegree=class2[0];
       if(afterdegree==beforedegree)
           effflag=1;
    }
}


long int nufindExpectedH( 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 nufindH0( int class2[])
{   int i, class3[PICNMBR], class4[PICNMBR];
    long int myh0;
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=0;
	for(i=0; i<PICNMBR; i+=1)
          class4[i]=class2[i];
	nugotofd(class4, class3);
	nuelimfcs(class4);
	if(class4[0]<0)
	      myh0=0;
    else
	      myh0=nufindExpectedH(class4);
   if(myh0<0)
          myh0=0;
    return myh0;
}
*/


long int nufindH0( int class2[])
{   int i, j, class3[PICNMBR], class4[PICNMBR];
    long int myh0;
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=0;
	for(i=0; i<PICNMBR; i+=1)
          class4[i]=class2[i];
	nugotofd(class4, class3);
	nuelimfcs(class4);
	if(class4[0]+class4[1]<0)
	      myh0=0;
    else
         {if(class4[0]+class4[1]+class4[2]<0)
              {j=class4[0]+class4[1]+class4[2];
              class4[0]=class4[0]+j;
              class4[1]=class4[1]-j;
              class4[2]=class4[2]-j;
               }
	      myh0=nufindExpectedH(class4);
         }
    return myh0;
}
/*
long int nufindH0( int class2[])
{   int i, class3[PICNMBR], class4[PICNMBR];
    long int myh0;
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=0;
	for(i=0; i<PICNMBR; i+=1)
          class4[i]=class2[i];
	nugotofd(class4, class3);
	nuelimfcs(class4);
	if(class4[0]-class4[1]<0)
	      myh0=0;
    else
         {if(class4[0]+class4[1]+class4[2]<0)
              {
              class4[0]=class4[0]+(class4[0]+class4[1]+class4[2]);
              class4[1]=class4[1]+(class4[0]+class4[1]+class4[2]);
              class4[2]=class4[2]+(class4[0]+class4[1]+class4[2]);
               }
	      myh0=nufindExpectedH(class4);
         }
    return myh0;
}
  
*/


long int nufindH1( int class2[])
{   int i, class3[PICNMBR], class4[PICNMBR], class5[PICNMBR];
    long int myh1;
	for(i=0; i<PICNMBR; i+=1)
          {class3[i]=class2[i];
          class4[i]=class2[i];
          class5[i]=1-class2[i];}
    class5[0]=-4+class5[0];
	myh1=nufindH0(class3)-nufindExpectedH(class4)+nufindH0(class5);
    return myh1;
}

long int nuqstar( int class2[])
{   int i, class3[PICNMBR], class4[PICNMBR];
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=class2[i];
	for(i=0; i<PICNMBR; i+=1)
          class4[i]=0;
    numypermute(class3,class4);
    class3[1]=class3[1]-1;
    return nufindH1(class3);
}

long int nulstar( int class2[])
{   int i, class3[PICNMBR], class4[PICNMBR];
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=class2[i];
	for(i=0; i<PICNMBR; i+=1)
          class4[i]=0;
    numypermute(class3,class4);
    class3[1]=class3[1]+1;
    class3[0]=class3[0]-1;
    return nufindH1(class3);
}

long int numxrkcok( int class2[])
{   int i, class3[PICNMBR];
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=class2[i];
    class3[0]=class3[0]+1;
    return nufindH0(class3)-3*nufindH0(class2);
}

void numakemonotone( int class2[])
{   int i, class3[PICNMBR];
	for(i=0; i<PICNMBR; i+=1)
          class3[i]=0;
 	numypermute( class2, class3);
 }


void nufreepart( int class2[]) 
{   int testflag;
testflag=0;
while(testflag==0)	
{numakemonotone(class2);
if(nuesix(class2)<0)
	nusix(class2);
 else
   {if(nuefive(class2)<0)
		nufive(class2);
	else
		{if(nuefour(class2)<0)
			nufour(class2);
		else
			{if(nuethree(class2)<0)
				nuthree(class2);
			else
				{if(nuetwo(class2)<0)
					nutwo(class2);
				else
					{if(nueone(class2)<0)
						nuone(class2);
					else
						{if(class2[8]>0)
							class2[8]=0;
						else
						    testflag=1;}}}}}}}
						
 }

int nuonesteptoample( int class2[])   /* assumes class2 is nef, nontrivial  */
{int x;
x=0;
numakemonotone(class2);
if(nuesix(class2)==0)
	nusix(class2);
 else
   {if(nuefive(class2)==0)
		nufive(class2);
	else
		{if(nuefour(class2)==0)
			nufour(class2);
		else
			{if(nuethree(class2)==0)
				nuthree(class2);
			else
				{if(nuetwo(class2)==0)
					nutwo(class2);
				else
					x=1;
}}}}
return x;}

int numinusdotK( int class2[]) 
{int x, i;
x=0;
for(i=0; i<PICNMBR; i+=1)   
	x=x+class2[i];
x=x+2*class2[0];
return x;}


int nuabs( int x)
{int y;
y=x;
if(x<0)
	y=-x;
return y;}