/* 
   moment_two.c 8-5-1997 

   Only print out the first two moments.
------------------------------------------ 
  
   Execute the program as the followings.
   gcc -o 'execute file' moment_two.c -lm
   'execute file' m 'output file'
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct cnode{
  int j,    /* the values of b(j,k). */
      k; 
  double a; /* coefficients of the weight. */     
  struct cnode *N,*P;
}LIST;

typedef struct lnode{
  int **mat; /* DE matrix */
  int **se;  /* starting and ending index matrix */
  int row,  /* rows of the matrix. */     
      num,  /* 1's in side the matrix */ 
      bk;   /* the number of block inside the matrix. */
  LIST *C;
  struct lnode *N,*P;
}LINK;

typedef struct tnode{
  int *e;
  int row;  /* dimension of the matrix. */
  LIST *C;
  struct tnode *L,*R;
}TREE;

typedef struct{ /* comparing the blocks inside the matrix */
  int **t;      /* the submatirx of the block */
  int len;      /* length of the block */
}com;

typedef struct pnode{
  int s,t;  /* the values in R(s,t) */
  LIST *C;
  struct pnode *L,*R;
}PP;

FILE *fp;
LIST *W;
TREE *S;
PP *R;
int clump;
int **New_Array(int,int);
LINK *Lalloc(void);
TREE *Talloc(void);
LIST *Calloc(void);

main(argc, argv)
int argc;
char *argv[];
{
   double ca,tca;   
   int **nmat;
   int m,I,M,i,j,k,ii,jj,incre; 
   int colnum,cj,ck;
   LINK *Head;
   LIST *c;
   double acoeff(int,int);
   LIST *listfree(LIST *); 
   LINK *translate(LINK *,int **,double,int,int,int,int); 
   LINK *decom(LINK *);
   LIST *bcoef(LIST *);
   PP  *Pfree(PP *);
   TREE *Delete(TREE *);
   void print_tree(TREE *);
   void polyprint(PP *);

   

   if (argc != 3)
     {
       fprintf(stderr, "Usage: %s m outfilename\n", argv[0]);
       return 1;
     }


   if((fp=fopen(argv[2],"a+"))==NULL){     /* data file not found */
       fprintf(stderr, "Cannot open data file %s\n", argv[2]);
       exit(1);
   }

   m = atoi(argv[1]);

   clump=m-1;


   fprintf(fp,"m:=%d;\n",m);

   for(M=2;M>=1;M--){
      Head=NULL;
      S=NULL;
      R=NULL;
      W=NULL;
      if(M==1)
         fprintf(fp,"m1");
      else if(M==2)
         fprintf(fp,"m2");      
      fprintf(fp,":='\n");
      for(I=1;I<=M;I++){
         tca=acoeff(M,I);
         if(I==1){
             cj=m-1;
             ck=1;
             ca=tca;             
             colnum=clump;       
             nmat=New_Array(I,colnum);       
             for(i=0;i<colnum;i++)
                 nmat[0][i]=1;
             Head=translate(Head,nmat,ca,cj,ck,I,colnum);
         }else if(I==2){
             for(i=1;i<=clump;i++){         
                if(i==clump){
                    cj=2*m-3;
                    ck=2;
                }else if(i!=clump){
                    cj=m+i-1;
                    ck=1;
                }
                ca=tca;                
                colnum=clump+i;
                nmat=New_Array(I,colnum);
                for(ii=0;ii<I;ii++){
                    for(jj=0;jj<colnum;jj++){
                        nmat[ii][jj]=0;
                    }
                }
                for(ii=incre=0;ii<colnum && incre<clump;ii++,incre++){
                    nmat[0][ii]=1;
                    nmat[1][ii+i]=1;
                }
                Head=translate(Head,nmat,ca,cj,ck,I,colnum);
             }
         }
      }

      if(Head!=NULL)
         Head=decom(Head);      
      
      print_tree(S);

      S=Delete(S);
      
      W=bcoef(W);

      fprintf(fp,"(\n");
      for(c=W;c!=NULL;c=c->N){
         if(c->a>0)
            fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k);
         else if(c->a<0) 
            fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k);
      }
      fprintf(fp,")\n");
      
      W=listfree(W);
      polyprint(R);
     
      R=Pfree(R);     
      fprintf(fp,"';\n"); 
   }
   fclose(fp);
}



/* New_Array: make an allocation for a matrix. */
int **New_Array(int row, int col)
{ 
  int **A,i,j;

  A=(int **)malloc(row*sizeof(int*));
  if(A==NULL)
     return NULL;
  for(i=0;i<row ;i++){
      A[i]=(int *)malloc(col*sizeof(int));
      if(A[i]==NULL){
        for(j=0;j<i;j++)
            free(A[i]);
        free(A);
        return NULL;
      }
  }
  return A;
}

/* Lalloc: make a link node for decomposed matrices. */
LINK *Lalloc(void)
{
  return (LINK *)malloc(sizeof(LINK));
}

/* Talloc: make a tree node for decomposed matrices. */
TREE *Talloc(void)
{
  return (TREE *)malloc(sizeof(TREE));
}

/* Calloc: make a list node for coefficents j, k, and a. */
LIST *Calloc(void)
{
  return (LIST *)malloc(sizeof(LIST));
}

/* acoeff: the recursion of a(k,i)=i*a(k-1,i)+i*a(k-1,i-1). */
double acoeff(int a,int b)
{
    if(a==1 && b==1)
       return 1.0;
    else if (b > a)
        return 0.0;
    else if (a<=0 || b<=0)
        return 0.0;
    else        
        return b*acoeff(a-1,b)+b*acoeff(a-1,b-1);
} 


LINK *translate(LINK *rt,int **d,double ca,int cj,int ck,int row,int col)
{
  int asize,i,j,k,ii,jj;
  int **a,**b,*delcl,*indices,**temp;
  int star,location,jump,count,ncol;
  LINK *NEW;
  LIST *cnew;
  double Comb(double,double);
  int **block(LINK *);
  LINK *createlink(LINK *,LINK *);
  LIST *combine(LIST *,LIST *);
  void permute(int,int,int *,int,int **,int *);
  
  
  cnew=Calloc();
  if(cnew==NULL){     
     fprintf(fp,"no room for temp(cplist)\n");
     exit(0);
  }
  cnew->j=cj;
  cnew->k=ck;
  cnew->a=ca;
  cnew->N=cnew->P=NULL;
  if(W==NULL)
     W=cnew;
  else
     W=combine(W,cnew);

  delcl=(int *)malloc(col*sizeof(int));
  if(delcl==NULL){
     fprintf(fp,"No room for delcl(translate)\n");
     exit(0);
  }
  indices=(int *)malloc(row*sizeof(int));
  if(indices==NULL){
     fprintf(fp,"No room for indices(translate)\n");
     exit(0);
  }
  for(k=1;k<=row;k++){
    ca *= -1.0;
    location=0;
    star=k;
    asize=Comb(row,star); 
 /* This is achievable job since the maximum row is 4.*/
    a=New_Array(asize,row);
    if(a==NULL){
      fprintf(fp,"No room for a(translate)\n");
      exit(0);
    }
    for(i=0;i<row;i++)
        indices[i]=0;
    permute(row,star,indices,0,a,&location);
    temp=New_Array(k,col);
    if(temp==NULL){
      fprintf(fp,"No room for temp(translate)\n");
      exit(0);
    }
    for(ii=0;ii<asize;ii++){
       for(jj=i=0;jj<row &&i<k;jj++){
           if(a[ii][jj]==1){
              for(j=0;j<col;j++){
                  temp[i][j]=d[jj][j];
	      }
              i++;
           }
       }
       /* delete the reduntant columns in temp matrix */
       ncol=0;         /* initial the value of ncol. */
       for(j=0;j<col;j++)
           delcl[j]=0;
       for(i=0;i<col;i++){
           for(j=jump=0;j<k;j++){
               jump += temp[j][i];
           }
           if(jump==0)
              delcl[i]=1;
       }   
       for(i=count=0;i<col;i++){
           if(delcl[i]==1)
              continue;
           for(j=0;j<k;j++)
               temp[j][count]=temp[j][i];
           count++;
       }
       ncol=count;  /* the dimension of new matrix will be k*ncol */
       /* transfer the matrix into DE form */
       b=New_Array(k,3);
       if(b==NULL){
          fprintf(fp,"No room for b(translate)\n");
          exit(0);
       }
       /*
          b: b[i][0] the number of zeros before the leading 1 
             b[i][1] the number of 1's=clump
             b[i][2] = b[i][0]+b[i][1]-b[i-1][0]b[i-1][1]
                     = b[i][0]-b[i-1][0]
         and b[0][2]=0
             the above b[i][2] is only good for inside the block..........
       */
       for(i=0;i<k;i++){                  
	   b[i][0]=0;        /* initial the values of b[i][0]=0 */
	   b[i][1]=clump;    /* since it calculates moments, b[i][1]=clump */
	   for(j=0;j<ncol;j++){   /* count the zeros before leading 1 */
	       if(temp[i][j]==0)
		  b[i][0]++;
	       else
		  break;
	   }
	   if(i==0)
	      b[i][2]=0;
	   else{
              if(b[i][0]==b[i-1][0]+b[i-1][1])  /* if nooverlap b[i][2]=0 */
                 b[i][2]=0;
              else
	         b[i][2]=b[i][0]-b[i-1][0];
	   }
       }
       /* create the new term into link list */
       NEW=Lalloc();
       if(NEW==NULL){
          fprintf(fp,"No room for NEW(translate)\n");
          exit(0);
       }
       NEW->mat=b; /* cp the DE matrix to node */
       NEW->row=k;
       NEW->num=0; /* count 1's inside the matrix */
       for(i=0;i<NEW->row;i++)
           NEW->num += NEW->mat[i][1];
       NEW->bk=0;   /* initialize the block be 0 */
       NEW->N=NEW->P=NULL;
       NEW->se=block(NEW);  /* examine the blocks inside the matrix */
       NEW->C=Calloc();
       if(NEW->C==NULL){
          fprintf(fp,"No room for NEW->C(translate)\n");
          exit(0);
       }
       NEW->C->j=cj;
       NEW->C->k=ck;
       NEW->C->a=ca;
       NEW->C->N=NEW->C->P=NULL;
       rt=createlink(rt,NEW); 
    }
    for(i=0;i<k;i++)
       free(temp[i]);
    free(temp);
    for(i=0;i<asize;i++)
        free(a[i]);
    free(a);
  }
  free(indices);
  free(delcl);
  return rt;
}

/* permutation of n elements with m of them 1's and the other 0's */
/* This is equivalent to the m combination of n indices, if we label 
   each of the n elements with a different index.  */
void permute(int n,int m,int *indIndices,int whichIndex,int **a,int *location)
{
   int i;
 
   if (whichIndex >= m)
     {
       for (i = 0; i < n; i++)
          {
            int j;
            int value = 0;  /* first assumes the value as 1 */

            /* if i is not in the indIndices, then it's 0 */
            for (j = 0; !value && j < whichIndex; j++)
               {
                  value = i == indIndices[j];
               } 
            a[*location][i]=value;
            if(i==n-1)
               (*location)++;
          }
       return;
     }
   for (i = whichIndex == 0 ? 0 : indIndices[whichIndex - 1] + 1;
        i <= n - (m - whichIndex);
        i++)
      {
        indIndices[whichIndex] = i;
        permute(n, m, indIndices, whichIndex + 1,a,location);
      }
   return;
}

double Fact(double n)
{
   if(n==0.0)
     return 1.0;
   else
     return n*Fact(n-1.0);
}

double Comb(double a, double b)
{
   int i,ans=1;
   double Fact(double);

   if(a < b)
     return 0.0;
   else{
     for(i=b;i>=1.0;i--,a--)
       ans=ans*a;    
     ans=ans/Fact(b); 
     return ans;
   }
}

/* createlink: add the node to the linklist */
LINK *createlink(LINK *rt,LINK *NEW)
{
   int i,j,match;
   LINK *cur,*prev;
   LINK *addchecking(LINK *,LINK *,LINK *,LINK *);
   LIST *combine(LIST *,LIST *);
   TREE *addtree(TREE *,LINK *);


/* check if the new matrix is simple matrix or not */
   if(NEW->bk==0){
      /*simple matrix call the shorthand fun. */
       S=addtree(S,NEW); 
       return rt;                
   }else{
       if(rt==NULL){
          return NEW;
       }else{
          for(prev=NULL,cur=rt;cur!=NULL;prev=cur,cur=cur->N){
            /* larger size of matrix first be decomposed */
              if(NEW->row<cur->row)
                 continue;
              else if(NEW->row>cur->row)
                 return addchecking(rt,prev,cur,NEW);
              else if(NEW->row==cur->row){
            /* more 1's inside the matrix first be decomposed */
                 if(NEW->num<cur->num)
                    continue;
                 else if(NEW->num>cur->num)
                    return addchecking(rt,prev,cur,NEW);
                 else if(NEW->num==cur->num){
            /* larger size of block is decomposed first */
                    if(NEW->se[NEW->bk-1][2]<cur->se[cur->bk-1][2])
                       continue;
                    else if(NEW->se[NEW->bk-1][2]>cur->se[cur->bk-1][2])
                       return addchecking(rt,prev,cur,NEW);
                    else if(NEW->se[NEW->bk-1][2]==cur->se[cur->bk-1][2]){
                       /*check if the identical term */
                       for(i=0,match=1;i<NEW->row && match==1;i++){
                           for(j=0;j<3 && match==1;j++){
                               if(NEW->mat[i][j]!=cur->mat[i][j])
                                  match=0;
                           }
                       }
                       if(match==0) /* Not identical, move to next node to check. */
                          continue;
                       else{       /* Identical, add coeffs to LIST */
                          cur->C=combine(cur->C,NEW->C);
                          for(i=0;i<NEW->bk;i++)   /* free the node */
                              free(NEW->se[i]);
                          free(NEW->se);
                          for(i=0;i<NEW->row;i++)
                              free(NEW->mat[i]);
                          free(NEW->mat);
                          free(NEW);
                          return rt;
	               }
                    }  
                 }
              }
          } 
          /* end of for-loop */
          /*connect the new node to the end of list */
          prev->N=NEW;
          NEW->P=prev;
          return rt;
       } /* end of else-when rt!=NULL */
   }
}

/* addchecking: add the node into the appropriate place */
LINK *addchecking(LINK *rt,LINK *prev,LINK *cur,LINK *NEW)
{
   if(prev==NULL){
      NEW->N=cur;
      cur->P=NEW;
      return NEW;
   }else{
      NEW->N=cur;
      NEW->P=prev;
      prev->N=NEW;
      cur->P=NEW;
      return rt;
   }
}

/* combine: combine the two list */
LIST *combine(LIST *old,LIST *new)
{
    int jp;
    LIST *p,*c,*temp;
    LIST *listcheck(LIST *,LIST *,LIST *,LIST *);

   /* old has at least one node in the LIST */
    while(new!=NULL){
        temp=new;
        new=new->N;
        if(new!=NULL)
           new->P=NULL;
        temp->N=NULL;
        for(p=NULL,c=old,jp=0;c!=NULL;p=c,c=c->N){
            if(temp->k<c->k)
               continue;
            else if(temp->k>c->k){
               old=listcheck(old,p,c,temp);
               jp=1;
               break; 
            }else if(temp->k==c->k){
               if(temp->j<c->j)
                  continue;
               else if(temp->j>c->j){
                  old=listcheck(old,p,c,temp);
                  jp=1;
                  break;
               }else if(temp->j==c->j){
                  c->a += temp->a;  
                  jp=1;
                  free(temp);
                  break;
               }

            }
        } 
        if(jp==0){
          p->N=temp;
          temp->P=p;
	}
    }
    return old; 
}

/* listcheck: add the node into the appropriate place */
LIST *listcheck(LIST *rt,LIST *prev,LIST *cur,LIST *NEW)
{
   if(prev==NULL){
      NEW->N=cur;
      cur->P=NEW;
      return NEW;
   }else{
      NEW->N=cur;
      NEW->P=prev;
      prev->N=NEW;
      cur->P=NEW;
      return rt;
   }
}

/* block: check the blocks inside the matrix and return the number of block and matrix */
int **block(LINK *b)
{

   int i,j,k,l,kk,ll,end;
   int *bk=&b->bk;  /* initialize to be 0 */
   int **temp;
   com *c;
   int Num_Cmp1( const void * , const void * );
   int Num_Cmp2( const void * , const void * );

   if(b->row==1)
     return NULL;
   /* record the matrix where is not overlapped with other rows */ 
   temp=New_Array(b->row,3);      
   if(temp==NULL){
      fprintf(fp,"No room for temp(block)\n");
      exit(0);
   }
   /*
      Good for b->row > 2 
      When b->mat[i][2]=0 and b->mat[i+1][2]!=0 it starts the block 
           b->mat[i][2]!=0 and b->mat[i+1][2]=0 it ends the block 
      b->se[i][0] records the starting row
      b->se[i][1] records the ending row
      b->se[i][2] records the number of rows inside the block
   */ 
   for(k=0,i=end=1;i<b->row;i++){
       if(b->mat[i][2]==0 && b->mat[i-1][2]==0){  
	  for(ll=0;ll<3;ll++)   /* cp the nonoverlapped row */
	    temp[k][ll]=b->mat[i-1][ll];
	  k++;                                    
       }else if(b->mat[i][2]!=0 && b->mat[i-1][2]!=0){  /* still inside the block */
	  continue;
       }else if(b->mat[i][2]==0 && b->mat[i-1][2]!=0){ /* it ends the block */
	  b->se[*bk][1]=i-1;             /* the ending row */  
	  b->se[*bk][2]=i-b->se[*bk][0]; /* the number of rows inside this block */
	  (*bk)++;                       /* increment to next block */
	  end=1;                         /* index for ending the block */
       }else if(b->mat[i][2]!=0 && b->mat[i-1][2]==0){ /* it starts the block */
	  if(*bk==0){
	     b->se=New_Array((int)b->row/2,3); /* The maximum blocks will be row/2 */
	     if(b->se==NULL){
		fprintf(fp,"No room for b->se(block)\n");
		exit(0);
	     }
	  }
	  b->se[*bk][0]=i-1; /* the staring row of the first block */
	  end=0;             /* index for starting the block */
       }
       if(i==b->row-1 && b->mat[i][2]==0){  /* examine the last row */ 
	 for(ll=0;ll<3;ll++)
	    temp[k][ll]=b->mat[i][ll];
	 k++;
       }
   }
   if(end==0){    /* completes the record of the block */
      b->se[*bk][1]=b->row-1;
      b->se[*bk][2]=b->row-b->se[*bk][0];
      (*bk)++;
   }
  /* if it is simple, sort(increasing) the number of 1's among rows */ 
   if(*bk==0 && k==b->row){  
      qsort(b->mat,b->row,sizeof(int *),Num_Cmp1);
      for(i=0;i<b->row;i++){
	  if(i==0)
	     b->mat[i][0]=0;
	  else
	     b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1];
      }
      b->se=NULL;
   }else{       /* There is at least one block inside the matrix */
      for(i=*bk;i<(int)b->row/2;i++) /* Free the un-used part */
	 free(b->se[i]);
      /* if k!=0 then put those unoverlapped rows to the front and sort them
         Also compare the size of the blocks and sort them at the back */
      if(k!=0){      
         if(b->se[0][0]!=k){
	    for(i=0,kk=k;i<*bk;i++){
	        l=kk;
	        for(j=b->se[i][0];j<=b->se[i][1];j++){
	            for(ll=0;ll<3;ll++)
		        temp[kk][ll]=b->mat[j][ll];
	            kk++;
	        }
	        b->se[i][0]=l;
	        b->se[i][1]=kk-1;
	    }
	    for(i=0;i<b->row;i++){
	        for(ll=0;ll<3;ll++)
	            b->mat[i][ll]=temp[i][ll];
	    }
            qsort(b->mat,b->se[0][0],sizeof(int *),Num_Cmp1);
	    for(i=0;i<b->se[0][0];i++){
	        if(i==0)
	           b->mat[i][0]=0;
	        else
	           b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1];
	    }
	 }
      }    
      if(*bk>=2){     
         c=(com *)malloc(*bk*sizeof(com));
         if(c==NULL){
            fprintf(fp,"No c space(block)\n");
            exit(0);
         }
         for(i=0;i<*bk;i++){
             c[i].t=New_Array(b->se[i][2],3);
             if(c[i].t==NULL){
                fprintf(fp,"No c.t space(block)\n");
	        exit(0);
	     }
	     for(j=0,kk=b->se[i][0];kk<=b->se[i][1];j++,kk++){
	         for(ll=0;ll<3;ll++)
		     c[i].t[j][ll]=b->mat[kk][ll];
	     }
	     c[i].len=b->se[i][2];
	 }       
         qsort(c, *bk, sizeof(com), Num_Cmp2);
	 for(i=0,j=b->se[0][0];i<*bk;i++){
	     for(kk=0;kk<c[i].len;j++,kk++){
		 for(ll=0;ll<3;ll++)
		     b->mat[j][ll]=c[i].t[kk][ll];
	     }
             /* rearrange the values inside b->se */
	     b->se[i][2]=c[i].len;       
             if(i!=0)
	        b->se[i][0]=b->se[i-1][1]+1;
             else
                b->se[0][0]=k;
             b->se[i][1]=b->se[i][0]+b->se[i][2]-1; 
             for(kk=0;kk<c[i].len;kk++)
                 free(c[i].t[kk]);
	     free(c[i].t); 
          
	 }
         free(c);  
      } 
      /* rearrange the values in b->mat[i][0] */   

       for(i=0,j=b->se[0][0];i<*bk;i++){	   
          for(kk=0;kk<b->se[i][2];kk++,j++){
	      if(kk==0){
                 if(i==0 && k==0)
                    b->mat[0][0]=0;
                 else
                    b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1];                
	      }else
	          b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1]+b->mat[j][2]-b->mat[j][1];
	  }
      }
   }
   for(i=0;i<b->row;i++)
       free(temp[i]); 
   free(temp);
   return b->se;
}

int Num_Cmp1( const void * S1, const void * S2)
{
   int ** T1 = (int **) S1,
       ** T2 = (int **) S2;

   if(*(*T1+1) < *(*T2+1))
      return -1;
   else if (*(*T1+1) > *(*T2+1))
      return 1;
   else
      return 0;
}

int Num_Cmp2( const void * S1, const void * S2)
{
   com * T1 = (com *) S1,
       * T2 = (com *) S2;

   if(T1->len < T2->len)
      return -1;
   else if (T1->len > T2->len)
      return 1;
   else
      return 0;
}



TREE *addtree(TREE *p,LINK *I)
{
   int i,match;
   LIST *combine(LIST *,LIST *);

   if(p==NULL){
       p=Talloc();
       if(p==NULL){
          fprintf(fp,"No room for p(addtree)\n");
          exit(0);
       }
       p->e=(int *)malloc(I->row*sizeof(int));
       if(p->e==NULL){
          fprintf(fp,"No room for p->e(addtree)\n");
          exit(0);
       }
       for(i=0;i<I->row;i++)
          p->e[i]=I->mat[i][1];
       p->row=I->row;
       p->C=I->C; 
       p->L=p->R=NULL;
       for(i=0;i<I->row;i++)
           free(I->mat[i]);
       free(I->mat);
       /* I->se is NULL, no need to free it */
       free(I);
   }else if(I->row>p->row){  /* if row is greater, goes to Left Tree */
       p->L=addtree(p->L,I);
   }else if(I->row<p->row){  /* if row is lesser, goes to Right Tree */
       p->R=addtree(p->R,I);
   }else if(I->row==p->row){
       for(i=match=0;i<I->row;i++){
           if(I->mat[i][1]==p->e[i]){
              continue;
           }else if(I->mat[i][1]>p->e[i]){
              match=1;
              break;
	   }else if(I->mat[i][1]<p->e[i]){
              match=2;
              break;
	   }
       }
       if(match==1){       /* Not identical, insert to Left Tree */
           p->L=addtree(p->L,I);
       }else if(match==2){
           p->R=addtree(p->R,I); 
       }else if(match==0){
            /* combine two link list of coefficients */
            p->C=combine(p->C,I->C);
           for(i=0;i<I->row;i++)
               free(I->mat[i]);
           free(I->mat);
           free(I);
       }
   }
   return p;
}

/* decom: decompose the matrices inside the link list */
LINK *decom(LINK *p)
{
    int i,j,col,len,tg;
    int *loc;
    LIST *cr;
    int *coef(LINK *,int *,int *,int *);
    LINK *organ(LINK *,int,int,double);
    LINK *Lfree(LINK *);

    while(p!=NULL){   
      /* We only look at the last big block and
             put corresponding columns to coeffs inside loc */
          col=p->mat[p->row-1][0]+p->mat[p->row-1][1];
          loc=(int *)malloc(col*sizeof(int));
          if(loc==NULL){
 	     fprintf(fp,"No room for loc(decom)\n");
	     exit(0);
          }
          /* initialize the values of array loc */
          for(i=0;i<col;i++) 
	      loc[i]=0;
          tg=-1;  /* initial target=-1 return 0 means tg is 0 vector
                                              1 means tg is (1,0...0) */
          len=0; /* the number of nonzero coeffs in this block */
         
          loc=coef(p,loc,&tg,&len);
         
          /*
             replace the target vector into matrix and produce the new node 
             the corresponding coeffs are 1 and -1 continuedly 
             i.e. odd's have coef=1 even's have -1.
          */
          for(i=0;i<len;i++){
              p=organ(p,loc[i],tg,fmod(1.0+i,2.0)); 
          }
          free(loc);
          if(p->N!=NULL){
             p=p->N;
             p->P=Lfree(p->P);
          }else{
             p=Lfree(p);  
          }
   }
   return p;
}

/* Lfree: To free the node inside the link list */
LINK *Lfree(LINK *p)
{
    int i;
    LIST *cr;

    for(i=0;i<p->bk;i++)
        free(p->se[i]);
    free(p->se);
    for(i=0;i<p->row;i++)
        free(p->mat[i]);
    free(p->mat);
    cr=p->C;
    while(cr!=NULL){
       if(cr->N!=NULL){
          cr=cr->N;
          free(cr->P);
          cr->P=NULL;
       }else{
          free(cr);
          cr=NULL;
       }
    }
    free(p);
    return NULL;
}


/* coef: find the location of coefficients and return the taget */
int *coef(LINK *p,int *c,int *tar,int *l)
{
    int mr,k,sum,i,j;
    int *fill(LINK *,int,int *,int *);

    mr=k=p->row-1; /* start to mark the set at the last row of the block */
    for(sum=0;k>=p->se[p->bk-1][0];k--){
       sum=sum+p->mat[k][2];
       if(sum>=p->mat[mr][1]){
	  c=fill(p,mr,c,l);
          /* if equal continue to mark set 
                not equal means stop and tg=0 and return */
	  if(sum==p->mat[mr][1]){  
	     mr=k-1;   /* the next marking row and initialize sum=0 */
	     sum=0;
	  }else if (sum>p->mat[mr][1]){  /* break here, no need to continue mark set */
	     *tar=0;   /* target is  0 vector */
	     return c;
	  }
       }
       /* if m* is the first row of the block */
       if(mr==p->se[p->bk-1][0]){  
	  *tar=1;     /* target is (1,0,...0) */
	  return c;
       }
    }
    /* if m* is not the first row of the block */
    c=fill(p,mr,c,l);
    *tar=0;       /* target is  0 vector */
    return c;
}

/* fill: find the ending and starting columns of current marking row */
int *fill(LINK *p,int mr,int *c,int *l)
{
    if(mr==p->row-1){
       c[(*l)++]=p->mat[mr][0]+p->mat[mr][1]-1; /* the last ending column */
    }
    c[(*l)++]=p->mat[mr][0]; /* the starting column of current one  */
    c[(*l)++]=p->mat[mr][0]-1;   /* the ending column of next one */
    return c;
}

/* organ: organize the new produced matrix
       (also delet the redundant rows) and to the link list */
LINK *organ(LINK *p,int c,int tg,double r)
{
     int i,j,k,kk,l,jp;
     int **d,*delrw;
     LINK *NEW;
     LIST *cur;
     LIST *cplist(LIST *,LIST *,double);

    /* d is the replaced matrix accroding to tg */
     d=New_Array(p->row,3);
     if(d==NULL){
	fprintf(fp,"No room for d(organ)\n");
	exit(0);
     }
     for(i=0;i<p->row;i++)
	 for(j=0;j<3;j++)
	     d[i][j]=p->mat[i][j];
    /* depending on tg to change the values inside matrix 
       c is the current replaced column 
       Spiritly we move the new replaced column to the first column
    */
     kk=p->se[p->bk-1][0];
     for(i=kk;i<p->row;i++){
         /* when the replaced column is not inside the leading 0's */  
	 if(c>d[i][0]-1){   
            /* when the replaced column is inside the lasting 0's */  
 	    if(c>d[i][0]+d[i][1]-1){
               if(tg==0)
                  d[i][0] += 1; 
	       if(tg==1){
                 /* add 1 to the first row of the block and 0 in the rest */ 
		  if(i==kk)  
		     d[i][1] += 1;
		  else if(i!=kk)
		     d[i][0] += 1;
	       }
	    }else if(c<=d[i][0]+d[i][1]-1){  
            /* when the replaced column is inside the 1's 
               when tg=1 the first row not change,
               others, the number of 1's reduce 1 
                                     0's add 1             */  
	       if(tg==0 || (tg==1 && i!=kk)){
		  d[i][0] += 1;
		  d[i][1] -= 1;
	       }
	    }
	 }

         if(i==kk||(i!=kk && d[i][0]==d[i-1][0]+d[i-1][1]))
	    d[i][2]=0;
	 else if(i!=kk && d[i][0]<d[i-1][0]+d[i-1][1])
	    d[i][2]=d[i][0]+d[i][1]-d[i-1][0]-d[i-1][1];

         /* when the replaced column is inside the leading 0's, no need to change
            so we jump out of the loop */    
	 if(c<=d[i][0]-1) 
	    break;
     }

    /* Delete the redundant rows */
     delrw=(int *)malloc(p->se[p->bk-1][2]*sizeof(int));
     if(delrw==NULL){
	fprintf(fp,"No room for delrw(organ)\n");
	exit(0);
     }
     for(i=0;i<p->se[p->bk-1][2];i++)
	delrw[i]=0;
     for(i=1,k=kk+1;k<p->row;k++,i++){
	if(delrw[i]==1)
	   continue;
	if(d[k-1][0]==d[k][0]){
           if(d[k-1][1]<=d[k][1])
	      delrw[i]=1;
	}else if(d[k-1][0]<d[k][0]){
           if(d[k-1][0]+d[k-1][1]==d[k][0]+d[k][1])
              delrw[i-1]=1;
           else if(d[k-1][0]+d[k-1][1]<d[k][0]+d[k][1])
              continue;
	}
     }

     /* reduced the redudant rows and l be the exact rows left 
        Here needs to be very careful about d[k][2].          
        if k==kk (the first row to be deleted) d[k+1][2]=0;
        if k==p->row-1 (the last row to be deleted) no change needed;
        For other k's, if d[k+1][0]>=d[k-1][0]+d[k-1][1] (nonoverlapped)
                       then d[k+1][2]=0;
                       else d[k+1][2] += d[k][2]
     */
     for(i=0,l=k=kk;k<p->row;i++,k++){
	if(delrw[i]==1){            
           if(k==kk)
              d[k+1][2]=0;
           else if(k<p->row-1){         
              if(d[k+1][0] >= d[k-1][0]+d[k-1][1])
                 d[k+1][2]=0;
              else 
                  d[k+1][2] += d[k][2];
           } 
	   continue;
	}
        for(j=0;j<3;j++)
	     d[l][j]=d[k][j];
        l++;
     }
     /* rearrange the values in d[k][0] */
     if(kk!=0)
       d[kk][0]=d[kk-1][0]+d[kk-1][1];
     else
       d[kk][0]=0;
     for(i=kk+1 ;i<l;i++){
         if(d[i][2]!=0)
	    d[i][0]=d[i-1][0]+d[i-1][1]+d[i][2]-d[i][1];
         else
            d[i][0]=d[i-1][0]+d[i-1][1];
     }
     /* Now the dimension of inserted matrix is l*3 */       

     /* Insert the new into node and add to the link*/
     NEW=Lalloc();
     if(NEW==NULL){
        fprintf(fp,"No room for NEW(organ)\n");
        exit(0);
     }
     NEW->mat=New_Array(l,3);
     if(NEW->mat==NULL){
        fprintf(fp,"No room for NEW->mat(organ)\n");
        exit(0);
     }

     for(i=0;i<l;i++){
         for(j=0;j<3;j++){
             NEW->mat[i][j]=d[i][j];             
         }
     }

     NEW->row=l;
     NEW->num=0;
     for(i=0;i<NEW->row;i++)
        NEW->num += NEW->mat[i][1];
     /* Here is wasting time to do it, Might need to fix it later 
        Still hard to say which is better */
     NEW->bk=0;   /* initialize the block be 0 */
     NEW->N=NEW->P=NULL;
     NEW->se=block(NEW);  /* examine the blocks inside the matrix */
     /* copy down the list of p->C into NEW->C */ 
     NEW->C=NULL;
    
     NEW->C=cplist(NEW->C,p->C,r);

     p=createlink(p,NEW);  

     return p;
     
}

/* cplist: copy the list of old to new */
LIST *cplist(LIST *new,LIST *old,double sign)
{
    LIST *c,*temp,*tail;

    for(c=old;c!=NULL;c=c->N){
        temp=Calloc();
        if(temp==NULL){     
           fprintf(fp,"no room for temp(cplist)\n");
           exit(0);
        }
        temp->j=c->j;
        temp->k=c->k;
    /* the odd term with coef 1 and the even term with -1 */
        if(sign==1.0)             
           temp->a=c->a;
        else if(sign==0.0)
           temp->a=c->a*(-1.0);
        temp->N=temp->P=NULL;
      
        if(new==NULL){ /* the starting node in the list */
           new=temp;
        }else{         /* add the node to the last of the list */
           tail=new;
           while(tail->N!=NULL)
             tail=tail->N;            
           temp->P=tail;
           tail->N=temp;
	}
      
    }
    return new;
}
   
/*
 polynomials: transfer the matrix notation into modified polynomials
p(interaction)=\sum_{0<=k<=l}{k \choose k1,k2,...,km}R(a,k)
and R(a,k)={n \choose k}(1-sum ai)^(n-k)\prod ai^ki      
*/
void polynomials(TREE *p)
{
   int i,j,k,ktotal;
   double ans;
   int *b;
   LIST *c,*new;
   double Fact(double);
   LIST *cplist(LIST *,LIST *,double);
   PP *addpoly(PP *,LIST *,double,int,int);

   /* the array to store the values of k's inside the multinomials. */
   b=(int *)malloc(p->row*sizeof(int)); 
   for(i=0;i<p->row;i++)
       b[i]=0;
  /* from 0,...,0 to the values-1's inside p->e[i]
     start from the first one(add 1) up to its max (value-1) and set to 0
     then the second one (add 1) the move the first one(add 1) to its max
     and then set to 0; then move up 1 in the second one till it reached 
     its max and set to 0 on the second one and move to the third and 
     repeat the previous work again then move to the fourth and so on. */
   i=0;
   while(1){
       for(k=ktotal=0;k<p->row;k++)
           ktotal += b[k];
       /* calculate {ktotal \choose b[0] ... b[p->row-1]} */
       for(ans=Fact(ktotal),k=0;k<p->row;k++)
           ans /= Fact(b[k]);
       new=NULL;
       new=cplist(new,p->C,1.0);
       for(c=new;c!=NULL;c=c->N){
           c->a *= ans;
       }
       
       R=addpoly(R,new,ans,ktotal,p->row); 

       if(b[i]+1<p->e[i]) /* the first one increase step by step */
          b[i] += 1;
       else{
          b[i]=0;
          for(j=i+1;j<p->row;j++){ 
	      b[j] += 1;
	      if(b[j]<p->e[j]) /* not reach its max */
	         break;
	      else   /* reach the max set back to 0 and move to next one */
	         b[j]=0;
          }
          if(j==p->row)  /* increment reaches the last component */
	     break;
          i=0;   /* back to the first one component */
       }
   }
}

    
PP *palloc(void)
{
  return (PP *)malloc(sizeof(PP));
}

PP *addpoly(PP *o,LIST *d,double kcom,int s,int t)
{


   LIST *combine(LIST *,LIST *);
   PP *palloc(void);  
   LIST *c;

   if(o==NULL){
      o=palloc();
      if(o==NULL){
         fprintf(fp,"no room for o(addpoly)\n");
         exit(0);
      }
      o->C=NULL;
      o->C=d;
      o->s=s;
      o->t=t;
      o->L=o->R=NULL;
   }else if(t>o->t){
      o->L=addpoly(o->L,d,kcom,s,t);
   }else if(t<o->t){
      o->R=addpoly(o->R,d,kcom,s,t);
   }else if(t==o->t){
      if(s>o->s){
	 o->L=addpoly(o->L,d,kcom,s,t);
      }else if(s<o->s){
	 o->R=addpoly(o->R,d,kcom,s,t);
      }else if(s==o->s){
	 o->C=combine(o->C,d);
      }
   }
   return o;
}


void polyprint(PP *p)
{

    LIST *c;
    LIST *bcoef(LIST *);

   if(p!=NULL){
     polyprint(p->R);

     /* Now b(j,k)=b(j+1,k)+b(j+1,k-1) */

     p->C=bcoef(p->C); 
     fprintf(fp,"+(\n");
     for(c=p->C;c!=NULL;c=c->N){
         if(c->a>0)
            fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k);
         else if(c->a<0) 
            fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k);
     }
     fprintf(fp,")*");
     fprintf(fp,"R(%d,%d)\n",p->s,p->t);
     polyprint(p->L);
   }
}

/*
   bcoef: use  b(j,k)=b(j+1,k)+b(j+1,k-1) 
   find the largest value of j amongst k and when k=0 jump out 
   the loop. (the largest j is always in the first node among the
   same k. Traverse the list and find the last node and execute
   the fact and delete the node from the list and add two more
   new nodes into the list. Continue this process until it reaches
   the first node and move to next k.
*/
LIST *bcoef(LIST *p)
{
    int jj,kk,num;
    LIST *tail,*head,*end,*temp1,*temp2,*temp3;
    LIST *combine(LIST *,LIST *);
   
    for(head=p,tail=p->N; tail!=NULL;){
        num=0;
        jj=head->j;
        kk=head->k;
        if(kk==0)
           break;
        while(tail!=NULL){
           if(tail->k==kk){
              num++;
              if(tail->N==NULL)
                 end=tail;
              tail=tail->N;
	   }else{
              end=tail->P;
              break;
	   }
	}
        if(num==0){
           head=end->N;
           if(head!=NULL)
              tail=head->N;
           else
              tail=NULL;
	}else{
            tail=end;
            while(head->j!=tail->j){
               temp1=Calloc();
               if(temp1==NULL){
                  fprintf(fp,"No room for temp1(polyprint)\n");
                  exit(0);
	       }
               temp1->j=tail->j+1;
               temp1->k=tail->k;
               temp1->a=tail->a;
               temp1->P=NULL;
               temp2=Calloc();
               if(temp2==NULL){
                  fprintf(fp,"No room for temp2(polyprint)\n");
                  exit(0);
	       }
               temp2->j=tail->j+1;
               temp2->k=tail->k-1;
               temp2->a=tail->a;
               temp2->P=temp1;
               temp1->N=temp2;
               temp2->N=NULL;
               combine(p,temp1);
               temp3=tail->P;
               temp3->N=tail->N;
               tail->N->P=temp3;
               free(tail);
               tail=temp3;
            }
            if(head->N!=NULL){
               head=head->N;
               tail=head->N;
	    }else{
               tail=NULL;
	    }
	}
    }
    return p; 
}


/*
 Delete the item in the given tree
 Delete the item when there is no child.
*/
TREE  *Delete(TREE *p)
{
    int i;    
    TREE *temp;
    TREE *Tfree(TREE *);

    if (p!= NULL){
       if(p->R==NULL && p->L==NULL){
          p=Tfree(p);
          return p;
       }else if(p->R==NULL){
          temp=p;
          p=p->L;
          p=Delete(p);
          temp=Tfree(temp);
          return temp;
       }else if(p->L==NULL){
          temp=p;
          p=p->R;
          p=Delete(p);
          temp=Tfree(temp);
          return temp;
       }else{
          temp=p;
          p->L=Delete(p->L);
          p->R=Delete(p->R);
          temp=Tfree(temp);
          return temp;
       }   
    } else {
       return p;
    }
    
}


TREE *Tfree(TREE *p)
{
    
    LIST *listfree(LIST *);

    
    free(p->e); 
    p->C=listfree(p->C);
    free( p);
    return NULL;
}

LIST *listfree(LIST *p)
{
    LIST *c;
    
    c=p;
    while(c!=NULL){
        if(c->N!=NULL){
           c=c->N;
           free(c->P);
           c->P=NULL;
        }else{
           free(c);
           c=NULL;
	}
    }
    return NULL;
}


void print_tree(TREE *p)
{
    int i;
    LIST *c;
    void polynomials(TREE *);
    
    if(p!= NULL){
       print_tree(p->R);
          polynomials(p);          
       print_tree(p->L);
    }
}


PP  *Pfree(PP *p)
{

    if(p!= NULL){
       p->R=Pfree(p->R);
       p->L=Pfree(p->L); 
       if(p->R==NULL && p->L==NULL){
          p->C=listfree(p->C);
          free(p);          
       }   
   }
   return p;
}