/**********************************************************************************
This program asks for two positive integer inputs, n and m=mult, corresponding to
the homogeneous ideal I defining n general points of projective 2 space P2
(over an algebraically closed field k). This program then computes various
numerical characters for the m-th symbolic power I^{(m)} of I, the hilbert
function of I^{(m)}, and the number v_d of generators in each degree d in a 
minimal homogeneous set of generators for I^{(m)}. Note that the minimal free 
resolution of I^{(m)} over the homogeneous coordinate ring R=k[P2] of P2 has the form
0 -> F_1 -> F_0 -> I^{(m)} -> 0, where F_0 = \sum_d R[-d]^{v_d}. From this and the 
hilbert function of I^{(m)}, this program computes a minimal set of generators for F_1.

The numerical characters printed out are: alpha, the least degree in which I^{(m)} is
nonzero; beta, the least degree greater than or equal to alpha in which the forms 
in I^{(m)} of that degree do not have a common divisor of positive degree; delta,
which is harder to explain (I^{(m)} gives rise to a sheaf of ideals S on P2, and
each d-th twist S(d) is f_*(F_d) for a certain invertible sheaf F_d on the surface X
obtained by blowing up the n general points defined by I, where f: X -> P2 is the blow up;
then delta is the least degree d such that F_d is numerically effective); and tau, the 
degree such that the hilbert function of I^{(m)} first and forever after equals the hilbert
polynomial of I^{(m)}. The values of the hilbert function (i.e., the dimension of I^{(m)}_d) 
for degrees d from alpha to tau+1 are also printed out, along with Campanella's lower and 
upper bounds for the minimal number of generators of I^{(m)} in each degree. The actual 
number of generators in each degree are printed between Campanella's bounds as a sum; 
the first summand is the number of generators in that degree coming from the loss of 
fixed components as the degree increases, and the second summand, sig[d], gives the number that 
don't. [The number v_d of generators in each degree alpha < d <= tau +1 turns out to
be v_d=(h^0(X, F_d)-h^0(X,e_0+(F_{d-1})_f))+sig[d], 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(e_1+\cdots+e_n) 
is the divisor class corresponding to (I^{(m)})_d, and (F_d)_f is the class of the
mobile (or free) part of |F_d|. See my papers (B. Harbourne) on resolutions of fat point 
subschemes for details.] Also, the values of dim((R/I^{(m)})_d) are printed out 
for alpha <= d <= tau+1. Finally, the number of generators in each degree 
for F_0 and F_1 are also printed out.

A few last comments about the program: 100-element arrays are used to hold the values
of the hilbert function, etc., in degrees alpha <= d <= tau + 1, so if
tau+2-alpha is greater than 100 the program terminates after computing alpha, beta,
tau and delta. Also, the actual state of affairs when the number of points is greater than 9 
is only conjectural (as of now, 7-31-95); the results this program gives for 10 or more points 
are thus conjectural. For 9 or fewer points the results are justified in "Generators for 
Symbolic Powers of Ideals Defining General Points of P2", B. Harbourne, preprint, 1995.
***************************************************************************************/


#include "stdio.h"
#define USERLIMIT 100 /* increase this value to handle bigger problems */ 

main()
{  long int totlb, totub, totactual, fcnu, testflag, f1, f2, k, omega;
   long int n, mult, e0, e1, alpha, tau, delta, beta, temp, d, F, c0, c1, j;
   long int h[USERLIMIT], hfreeplus[USERLIMIT];
   long int cok[USERLIMIT], lb[USERLIMIT];
   long int ub[USERLIMIT], dh[USERLIMIT], ddh[USERLIMIT];
   long int dddh[USERLIMIT], eps[USERLIMIT];
   long int numgen[USERLIMIT], numsyzgen[USERLIMIT], sig[USERLIMIT];
   int q, maxmult;

/* Later, e0 and e1 will be used to hold the degree and multiplicity
of a uniform abnormal curve (apres Nagata) for n general points of P2;
here they are just initialized to 0, for the default case of an n with no
abnormal curve.  */


    e0=0;
    e1=0;

printf("\nGiven n general points of the projective plane,\n");
printf("and a multiplicity m, we will compute a resolution\n");
printf("over the homogeneous coordinate ring R=k[x,y,z] of P2\n");
printf("for the m-th symbolic power I^(m) of the homogeneous\n");
printf("ideal I defining the n points. Note that the result is still\n");
printf("conjectural for n>9.\n\n");

printf("The resolution takes the form 0 -> F_1 -> F_0 -> I^(m) -> 0;\n");
printf("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("Other information about I^(m) is also printed out.\n\n");

  printf( "\nEnter the number n of points: " ); 
   n=0;
scanf( "%ld", &n ); 
   while(getc(stdin)!='\n');
   if(n>9999)
      n=0;
   while(n<1) 
      {n=0;
       printf("Please enter a number between 1 and 9999 ... :\n");
       scanf( "%ld", &n );
       while(getc(stdin)!='\n');
       if(n>9999)
          n=0;
   }

  
    if (n>9) 
         printf("\nWARNING: Results for 10 or more generic points are conjectural!!!!!\n\n");

maxmult=0;
while(n*(maxmult+1)*(maxmult+1)<9999*9999)
maxmult++;

printf( "Enter the uniform multiplicity m at the points: " );
   mult=0;
scanf( "%ld", &mult ); 
   if(mult> maxmult)
      mult=0;
   while(getc(stdin)!='\n');
   while(mult<1) 
      {mult=0;
       printf("Please enter a number between 1 and %d ... :\n", maxmult   );
       scanf( "%ld", &mult );
       while(getc(stdin)!='\n');
       if(mult>maxmult)
          mult=0;
   }
               
/* Compute alpha, beta, tau and delta, depending on the number of points, n.
   Note that the computation depends on certain theorems not evident in the 
   program; it is not done from the definitions.	
*/

    if ( n==1)
         {alpha=mult;
          tau=mult-1;
          delta=mult;
          beta=delta;
         }

    if ( n==2)
         {alpha=mult;
          tau=2*mult-1;
          delta=2*mult;
          beta=delta;
          e0=1;
          e1=1;
         }

    if ( n==3)
         {alpha=(3*mult+1)/2;
          tau=2*mult-1;
          delta=2*mult;
          beta=delta;
          e0=3;
          e1=2;
         }

    if ( n==4)
         {alpha=2*mult;
          tau=2*mult;
          delta=2*mult;
          beta=delta;
         }


    if ( n==5)
         {alpha=2*mult;
          tau=5*mult/2;
          delta=(5*mult+1)/2;
          beta=delta;
          e0=2;
          e1=1;
         }

    if ( n==6)
         {alpha=(12*mult+4)/5;
          tau=5*mult/2;
          delta=(5*mult+1)/2;
          beta=delta;
          e0=12;
          e1=5;
         }

    if ( n==7)
         {alpha=(21*mult+7)/8;
          tau=(8*mult+1)/3;
          delta=(8*mult+2)/3;
          beta=delta;
          e0=21;
          e1=8;
         }

    if ( n==8)
         {alpha=(48*mult+16)/17;
          tau=(17*mult+4)/6;
          delta=(17*mult+5)/6;
          beta=delta;
          e0=48;
          e1=17;
         }

    if ( n==9)
         {alpha=3*mult;
          tau=3*mult;
          delta=3*mult;
          beta=delta+1;
         }

    if ( n>9)
         {tau=1;
          temp=6-n*(mult*mult+mult);
          while(temp<0)
             {temp=temp+2*tau+4;
              tau += 1;
             }
           alpha=tau;
           if(temp==0)
             alpha=tau+1;
           beta=alpha;
           if(temp==2)
             beta=alpha+1;
           delta=1;
           temp=1-n*mult*mult;
           while(temp<0)
              {temp += 2*delta+1;
              delta += 1;
              }
          }

	printf("delta= %ld  beta= %ld  tau= %ld  alpha= %ld\n", delta, beta, tau, alpha);      
    printf("\n");	
          
     
  if(tau-alpha+2>USERLIMIT)
    {printf("First generator in degree %ld, last in degree at least %ld.",alpha,beta);
   printf("\nUser Limits exceeded; try smaller example.\n\n");}
  else
   {sig[0]=0;
    for (d=0; d<=tau+1-alpha; d+=1) 
      {F=d+alpha;
       c0=F;
       c1=mult;
       if(n<10)
           while(c0*e0<n*c1*e1) {
            c0=c0-e0;
            c1=c1-e1;}
       else
            if(c0*c0+3*c0-n*c1*c1-n*c1==0)
               {c0=0;
                c1=0;
                }

 /* For each degree d>=alpha, the fixed components must be subtracted off */

          printf( "Free part in degree %ld has degree %ld and mult %ld\n",F,c0,c1); 

          sig[d+1]=0;

/* The number of generators not coming from the fixed components or computable from
maximal rank considerations are dealt with here, based on results from  "Generators for 
Symbolic Powers of Ideals Defining General Points of P2", B. Harbourne, preprint, 1995.*/


       if((n==7) && (c0*3==c1*8) && (c1>8))
          sig[d+1]=7;
       if((n==8) && (c0*6==c1*17) && (c1>53))
          sig[d+1]=48;
       if((n==8) && ((c0-3)*6==(c1-1)*17) && (c1>36))
          sig[d+1]=16; 

 /* h[d]=h^0(X,F_d) is the array holding the values of the hilbert function;
hfreeplus[d] is the array holding the values of h^0(X, (F_d)_m+e_0);
and cok[d]=dim((R/I^{(m)})_d)  */

       h[d]=(c0*c0+3*c0-n*c1*c1-n*c1)/2+1;
       hfreeplus[d]=h[d]+c0+2;
       cok[d]=(F*F+3*F)/2+1-h[d];             
      }

    printf("\n");	
/* To compute Campanella bounds on the number of generators in ech degree,
difference operators must be applied to the hilbert function, h.
Thus dh is the arry of first differences, ddh the second differences etc. */

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

/* In Campanella's result, numbers of generators are compared with a complete
intersection; eps[d] is an array accounting for the number of generators
in a complete intersection: 1 less for alpha <= d <beta, and 2 less
for beta <= d.  Then ub[d] and lb[d] are arrays holding Campanella's bounds. */

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

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

         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]<1)
     lb[beta-alpha]=1;  

/* totlb contains the sum of Campanella's lower bounds,
totub contains the sum of Campanella's upper bounds,
and totactual is the total actual number of generators. 
For each degree, fcnu is the number coming from loss of fixed components,
and sig[d] the rest. */

totlb=0;
totub=0;

totactual=0;

for(j=0;j<=tau+1-alpha;j+=1)
{totlb+=lb[j];
totub+=ub[j];
if(j==0)
        fcnu=h[j];
    else
       {fcnu=h[j]-hfreeplus[j-1];
        if((3*(h[j-1])+fcnu<h[j]) &&  sig[j]==0)
           sig[j]=h[j]-3*(h[j-1])-fcnu;
        }
        numgen[j]=fcnu+sig[j];
        totactual+=numgen[j];

        printf("deg=%ld h^0=%ld nu: %ld <= %ld + %ld <= %ld Hilb=%ld\n", j+alpha, h[j], lb[j], fcnu, sig[j], ub[j], cok[j]);
     }

/* Here the number of generators in each degree for F_1 is computed. The method is
to compare the difference in dimensions of (F_0)_d with (I^{(m)})_d
and with that part of (F_1)_d coming from already determined generators (if any)
of degree less than d. The routine stops when the number of generators 
accumulated for F_1 is 1 less than the number for F_0. */   

     j=-1;
     testflag=0;
     numsyzgen[0]=0;
     while(testflag<totactual-1)
     {f1=0;
      f2=0;
      j+=1;

      if(j<tau+2-alpha)
         {for(k=0;k<=j;k+=1)
             f1=f1+numgen[k]*( ((j-k)*(j-k)+3*(j-k))/2 +1);
          for(k=0;k<=j-1;k+=1)
             f2=f2+numsyzgen[k]*( ((j-k)*(j-k)+3*(j-k))/2 +1);
          numsyzgen[j]=f1-f2-h[j];
          if(numsyzgen[j] != 0)
            testflag=testflag+numsyzgen[j];
       }
       

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

     for(k=0;k<=tau+1-alpha;k+=1)
        if(numgen[k]>0)
            printf("Number of gens of F_0 in degree %ld is %ld\n", k+alpha, numgen[k]);

     printf("\n");

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