/* 
   scan.c 8-05-1997 
   
   gcc -o 'execute file' scan.c -lm
   'execute file' p k 'output file' 
   (p is the number of 1's in the each row
    k in the number of columns in a matrix)
*/

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

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. */
  double a; /* coefficients of the weight. */      
  struct lnode *N,*P;
}LINK;

typedef struct tnode{
  int *e;
  int row;  /* dimension of the matrix. */
  double a;
  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) */
  double a;
  struct pnode *L,*R;
}PP;

FILE *fp;
TREE *S;
PP *R;
double checker=0.0;

int **New_Array(int,int);
LINK *Lalloc(void);
TREE *Talloc(void);

main(argc, argv)
int argc;
char *argv[];
{
   int **nmat;
   int p,k; 
   LINK *Head;
   LINK *scanmatrix(LINK*, int, int);
   LINK *decom(LINK *);
   PP  *Pfree(PP *);
   TREE *Delete(TREE *);
   void print_tree(TREE *);
   void polyprint(PP *);

   time_t start_time,end_time,total_time;

   
   time(&start_time);

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


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

   p = atoi(argv[1]);

   k = atoi(argv[2]);

   Head=NULL;
   S=NULL;
   R=NULL;
   Head=scanmatrix(Head,p,k);


   if(Head!=NULL)
      Head=decom(Head);      
   else
      fprintf(fp,"Error in the scanmatrix function\n");

    print_tree(S); 
    
   fprintf(fp,"checker=%.0f\n",checker); 
   
          
   S=Delete(S);
   fprintf(fp,"p:=proc()\n");
   polyprint(R);
   fprintf(fp,"end;\n");
   R=Pfree(R); 
   
   time(&end_time);
   total_time=end_time - start_time;
 
   fprintf(fp,"Total run time = %d \n",total_time);

   fclose(fp);
}

/* scanmatrix: create a scan matrix to decompose */ 
LINK *scanmatrix(LINK *H,int p,int k)
{
    int i,j;   
    int **block(LINK *);

    fprintf(fp,"p=(m-1)=%d,k=(n-1)=%d\n",p,k);
    H=Lalloc();
    if(H==NULL){
	fprintf(fp,"No room for H(scanmatrix)\n");
	exit(0);
    }
    H->row=k-p+1;
    H->bk=0;
    H->num=H->row*p;
    H->N=H->P=NULL;
    H->mat=New_Array(H->row,3);
    if(H->mat==NULL){
       fprintf(fp,"No room for H->mat(scanmatrix)\n");
       exit(0);
    }
    for(i=0;i<H->row;i++){
       if(i)
	  H->mat[i][2]=1;
       else
	  H->mat[i][2]=0;
       H->mat[i][0]=i;
       H->mat[i][1]=p;
    }
  
    H->se=block(H); 
    H->a=1.0;
    return H;
}


/* 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));
}

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 *);
   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 */
                          cur->a += NEW->a;
                          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;
   }
}


/* 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;   


   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->a=I->a; 
       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  coefficients */
           p->a += I->a;
           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;    
    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;

    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);    
    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;

    /* 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 */ 

     if(r==1.0)             
        NEW->a=p->a;
     else if(r==0.0)
        NEW->a=p->a*(-1.0);
     p=createlink(p,NEW);  

     return p;
     
}



/*
 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,newans;
   int *b;   
   double Fact(double);
   PP *addpoly(PP *,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]);
       newans= ans*p->a;       
       R=addpoly(R,newans,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,double a,int s,int t)
{

   PP *palloc(void);  

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


void polyprint(PP *p)
{
    

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

/*
 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)
{

    free(p->e); 
    free( p);
    return NULL;
}

void print_tree(TREE *p)
{
    int i;
    void polynomials(TREE *);

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

       checker += p->a;
       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){
          free(p);          
       }   
   }
   return p;
}