/* 
    general_type.c 8-5-1997
   
    Modified to be good for both types of matrices.
   -------------------------------------------------------------
   This file prints out all the output execute the program as 
   the following.

   gcc -o 'execute file' general_type.c -lm
   'execute file' type 'input file' 'output file' 
   (type is the value of 1 or 2)
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.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 type;
   double aaa;
   FILE *fi;
   LINK *Head;
   LINK *read_matrix(FILE *,LINK *,int,double *);
   LINK *decom(LINK *);
   PP  *Pfree(PP *);
   TREE *Delete(TREE *);
   void print_tree(TREE *);
   void print_link(LINK *);
   void polyprint(PP *);

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

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

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

   type = atoi(argv[1]);

   Head=NULL;
   S=NULL;
   R=NULL;
   aaa=0.0;   
   Head=read_matrix(fi,Head,type,&aaa);
 /* print_link(Head); */


   if(Head!=NULL)
      Head=decom(Head);      
   else
      fprintf(fp,"Error in the input\n");
   
/*   fprintf(fp,"The simple link list\n");
   fprintf(fp,"-----------------\n");
   if(type==2)
      fprintf(fp,"%.0f\n",aaa); */

   print_tree(S);

   fprintf(fp,"checker=%.0f\n",checker); 
   
   S=Delete(S);
   
   fprintf(fp,"p:=proc()\n");

   if(type==2)
      fprintf(fp,"%.0f\n",aaa);
  
   polyprint(R);
   fprintf(fp,"end;\n");
   
   R=Pfree(R);     

   fclose(fi);
   fclose(fp);
}


/* read_matrix: read data into matrix and record mrow, mcol */
LINK *read_matrix(FILE *ff,LINK *p,int tt,double *aaa)
{
    int i,j,k,err;
    int row,col;
    double a;
    int **d;
    LINK *new;
    LINK *newnode(LINK*,int **,double,int,int);
    LINK *createlink(LINK *,LINK *);
    LINK *translate(LINK *,int **,double,int,int);

    fscanf(ff,"%d",&k);
    fprintf(fp,"# of matrices=%d\n",k);
    while(k>0){
        err=fscanf(ff,"%d %d %lf",&row,&col,&a);
        if(err==EOF){
           fprintf(fp,"End of file, check the value of k\n");
           exit(0);
	}
        d=New_Array(row,col);
        if(d==NULL){
           fprintf(fp,"No room for d(read_matrix)\n");
           exit(0);
	}
	fprintf(fp,"row=%d,col=%d,wt=%.0f\n",row,col,a);
        for(i=0;i<row;i++){
           for(j=0;j<col;j++){
               fscanf(ff,"%d ",&d[i][j]);
               fprintf(fp,"%d ",d[i][j]);
	   }
           fputc('\n',fp);
	}
        if(tt==1){
           new=NULL;
           new=newnode(new,d,a,row,col);
           p=createlink(p,new);
        }else if(tt==2){
           *aaa += a; /* add a to aaa */
           p=translate(p,d,a,row,col); 
        }
        k--;
   }
 
   return p;
}

LINK *newnode(LINK *new,int **d,double a,int row,int col)
{
    int i,j;
    int **block(LINK *);

    new=Lalloc();
    if(new==NULL){
       fprintf(fp,"No room for new(new_node)\n");
       exit(0);
    }
    new->row=row;
    new->bk=0;
    new->num=0;
    new->N=new->P=NULL;
    new->mat=New_Array(new->row,3);
    if(new->mat==NULL){
       fprintf(fp,"No room for new->mat(new_node)\n");
       exit(0);
    }
 /* transfer the matrix into DE form */
 /*
      b: b[i][0] the number of zeros before the leading 1 
         b[i][1] the number of 1's
         b[i][2] = b[i][0]+b[i][1]-b[i-1][0]b[i-1][1]
         and b[0][2]=0
         the above b[i][2] is only good for inside the block..........
 */      

    for(i=0;i<row;i++){
       new->mat[i][0]=0;
       for(j=0;j<col;j++){
          if(d[i][j]==0)
             new->mat[i][0]++;
          else
             break;
       }
       new->mat[i][1]=0;
       for(;j<col;j++){
          if(d[i][j]==1)
             new->mat[i][1]++;
          else
             break;
       }
       new->num += new->mat[i][1];
       if(i==0)
         new->mat[i][2]=0;
       else{
         if(new->mat[i][0]==new->mat[i-1][0]+new->mat[i-1][1]) /* nonoverlapped */
            new->mat[i][2]=0;
         else
            new->mat[i][2]=new->mat[i][0]+new->mat[i][1]-new->mat[i-1][0]-new->mat[i-1][1];
       }
    }
    new->se=block(new); 
    new->a=a;

    return new;
}

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


LINK *translate(LINK *rt,int **d,double ca,int row,int col)
{
   int asize,i,j,k,ii,jj;
   int **a,*delcl,*indices,**temp;
   int star,location,jump,count,ncol;
   LINK *NEW;
   double Comb(double,double);
   LINK *newnode(LINK *,int **,double,int,int);
   LINK *createlink(LINK *,LINK *);
   void permute(int,int,int *,int,int **,int *);
  
  
/*
   fprintf(fp,"ca=%.0f\n",ca);
   for(ii=0;ii<row;ii++){
        for(jj=0;jj<col;jj++){
           fprintf(fp,"%d ",d[ii][jj]);
        }
        fputc('\n',fp);
   }
*/


   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 */
           NEW=NULL;
           NEW=newnode(NEW,temp,ca,k,ncol);  
           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 *);
   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);

/*
          if(p->a>0)
              fprintf(fp,"+%.0f*",p->a);    
          else if(p->a<0)
              fprintf(fp,"%.0f*",p->a);    
	  if(p->a!=0){
             for(i=0;i<p->row;i++)
                fprintf(fp,"(%d)",p->e[i]);
             fputc('\n',fp); 
	  }
*/	  
          checker += p->a;
          polynomials(p); 
       print_tree(p->L);
    }
}


void print_link(LINK *p)
{   
    int i,j;

    if(p!=NULL)
       fprintf(fp,"The Input list\n");
    else
       fprintf(fp,"No Input list\n");
    fprintf(fp,"----------------------\n");
    while(p!= NULL){
          fprintf(fp,"%.0f\n",p->a);              
          fprintf(fp,"row=%d,bk=%d\n",p->row,p->bk);
          for(i=0;i<p->row;i++){
             for(j=0;j<3; j++){
                fprintf(fp,"%d ",p->mat[i][j]);
             }
             fputc('\n',fp); 
          }
          fprintf(fp,"----------------------\n");
       p=p->N;
    }
}

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