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

static void ReadInfo(int *pNumRow, int *pNumCol, int *pM, int *pN);
static void ReadTable (int NumRow, int NumCol, int **table, int *rowSum, 
		       int *colSum, int **s0, int *s0RowSum);
static int Reciprocal(int **a, int NumCol);
static void ReorderRowCol (int **table, int **s0, int *rowSum, int *colSum, 
			   int NumRow, int NumCol, 
			   int *initialPermute, int *initialColPermute, 
			   int *s0WhichCol, int *s0WhichRow, 
			   int *partialColSum);
static void  InitialSetting (double *pSampleMean, double *pSampleVar, 
			     double *pPValue);
static void GenerateATable (int **t, int **s0, int NumRow, int NumCol, 
			    int *rowSum, int *colSum, int *partialColSum, 
			    int *s0RowSum, int *s0WhichCol, int *s0WhichRow, 
			    int *initialPermute, int *initialColPermute, 
			    double *pCurrentSample);
static void UpdatePValue (double currentSample, int **t, 
			  int NumCol, double reciprocal0, 
			  double *pSampleMean, 
		          double *pSampleVar, double *pPValue);
static void PrintOutput (FILE *outFile, int N, double sampleMean, 
			 double sampleVar, double pValue);
static double  NChooseM(int n, int m);
static double Mean(double x[], int n);
static int RandomSample2(double weight[], int n);
static double Drafting(int w[], int start, int end, int n, int a[],
		       int s[], int NumRow);
static int ChooseOneInteger(double w[], int start, int end, int NumRow);
static double R(int k, double w[], int start, int end, int NumRow);
static double T(int i, double w[], int start, int end);
static double RandomReal(double low, double high);
static int RandomInteger(int low, int high);

main()
{      
        int NumRow, NumCol, M, N, numEstimate, numSample, i;
	int *rowSum, *colSum, *initialPermute, *initialColPermute, 
	  *s0WhichCol, *s0WhichRow, *partialColSum, *s0RowSum;	
        FILE *outFile;
	double sampleMean, sampleVar, currentSample, pValue, reciprocal0, 
	  cvsquare;       
	int **table, **t, **s0;

	ReadInfo(&NumRow, &NumCol, &M, &N);   

	//allocate space for pointers
	table=malloc(NumRow*sizeof(int *));
	t=malloc(NumRow*sizeof(int *));
	s0=malloc(NumRow*sizeof(int *));
	rowSum=malloc(NumRow*sizeof(int));
	colSum=malloc(NumCol*sizeof(int));
	s0RowSum=malloc(NumRow*sizeof(int));
	initialPermute=malloc(NumRow*sizeof(int));
	initialColPermute=malloc(NumCol*sizeof(int));
	s0WhichCol=malloc(NumRow*sizeof(int));
	s0WhichRow=malloc(NumCol*sizeof(int));
	partialColSum=malloc(NumCol*sizeof(int));
        for (i=0; i<NumRow; i++) {  
	  table[i]=malloc(NumCol*sizeof(int));
	  t[i]=malloc(NumCol*sizeof(int));
	  s0[i]=malloc(NumCol*sizeof(int));
        }	
       //end of allocation

	srand((int) time(NULL));

	ReadTable(NumRow, NumCol, table, rowSum, colSum, s0, s0RowSum);

	reciprocal0=Reciprocal(table, NumCol);

	ReorderRowCol(table, s0, rowSum, colSum, NumRow, NumCol, 
		      initialPermute, initialColPermute, s0WhichCol, 
		      s0WhichRow, partialColSum);

	outFile=fopen("result.txt", "w");
	if (outFile != NULL) {

	  for (numEstimate=0; numEstimate<M; numEstimate++) {
	    InitialSetting(&sampleMean, &sampleVar, &pValue);

	    for (numSample=0; numSample<N; numSample++) {

	      GenerateATable(t, s0, NumRow, NumCol, rowSum, colSum, 
			     partialColSum, s0RowSum, s0WhichCol, s0WhichRow, 
			     initialPermute, initialColPermute, 
			     &currentSample);

	      UpdatePValue(currentSample, t, NumCol, reciprocal0, 
			   &sampleMean, &sampleVar, &pValue);
	    }

	    PrintOutput(outFile, N, sampleMean, sampleVar, pValue);
	  }
	}
	else { printf("File creating failed!\n");}
	fclose(outFile);
}	           

static void ReadInfo (int *pNumRow, int *pNumCol, int *pM, int *pN)
{
        int NumRow, NumCol, M, N;

	printf("\nPlease enter the number of rows in the table?\n");
	scanf("%d", &NumRow);
	printf("\nPlease enter the number of columns in the table?\n");
	scanf("%d", &NumCol);
	printf("\nHow many estimates do you need?\n");
	scanf("%d", &M);
	printf("\nHow many samples should each estimate be based on?\n");
	scanf("%d", &N);	

	*pNumRow=NumRow;
	*pNumCol=NumCol;
	*pM=M;
	*pN=N;
}

static void ReadTable (int NumRow, int NumCol, int **table, int *rowSum, 
		       int *colSum, int **s0, int *s0RowSum)
{
  int i, j;
  FILE *inFile;
  char nameInFile[256];

	printf("\nPlease enter the data file name:\n");
	scanf("%s", nameInFile);
	while((inFile=fopen(nameInFile, "r"))==NULL)
	  {
	    printf("\nCan't read file. Please enter the data file name:\n");
	    scanf("%s", nameInFile);
	  }

        for (i=0; i<NumRow; i++) {  
          for (j=0; j<NumCol; j++) {
	    fscanf(inFile, "%d", &table[i][j]);	    
          }	  
        }	

	for (i=0; i<NumRow; i++) { 
	  rowSum[i]=0;
	  for (j=0; j<NumCol; j++) {
	    rowSum[i]=rowSum[i]+table[i][j];
	  }
	}
	for (j=0; j<NumCol; j++) { 
	  colSum[j]=0;
	  for (i=0; i<NumRow; i++) {
	    colSum[j]=colSum[j]+table[i][j];
	  }
	}

	for (i=0; i<NumRow; i++) { 
	  s0RowSum[i]=1;
	  for (j=0; j<NumCol; j++) {
	    if (i==j) {
	      s0[i][j]=1;
	    }
	    else s0[i][j]=0;
	  }
	}
}

static int Reciprocal(int **a, int NumCol)
{

	int i,j;
	int sum;

	sum=0;

	for (j=1; j<NumCol; j++) {
	  for (i=0; i<j; i++) {
	    sum = sum+a[i][j]*a[j][i];
	  }
	}

	return (sum);
}

static void ReorderRowCol (int **table, int **s0, int *rowSum, int *colSum,
			   int NumRow, int NumCol, 
			   int *initialPermute, int *initialColPermute, 
			   int *s0WhichCol, int *s0WhichRow, 
			   int *partialColSum)
{
  int i, k, m, n, temp, kk, temp1;

	for (i=0; i< NumRow; i++) {
	  initialPermute[i]=i;
	}      
	for (i=0; i< NumCol; i++) {
	  initialColPermute[i]=i;
	}

	//change order of columns
	for (m=0; m<NumCol; m++) {
	  k=m;
	  for (n=m+1; n<NumCol; n++) {
	    if (colSum[n]>colSum[k])
	      k=n;
	  }
	  if (k != m) {
	    temp=colSum[m];
	    colSum[m]=colSum[k];
	    colSum[k]=temp;

	    temp=initialColPermute[m];
	    initialColPermute[m]=initialColPermute[k];	
	    initialColPermute[k]=temp;

	    for (kk=0; kk<NumRow; kk++) {
	      temp1=table[kk][m];
	      table[kk][m]=table[kk][k];
	      table[kk][k]=temp1;

	      temp1=s0[kk][m];
	      s0[kk][m]=s0[kk][k];
	      s0[kk][k]=temp1;
	    }
	  }
	}

	for (m=0; m<NumRow; m++) {
	  s0WhichCol[m]=NumCol;
	  for (n=0; n<NumCol; n++) {
	    if (s0[m][n]==1) {
	      s0WhichCol[m]=n;
	      break;
	    }
	  }
	}

	//change order of rows
	for (m=0; m<NumRow; m++) {
	  k=m;
	  for (n=m+1; n<NumRow; n++) {
	    if ((rowSum[n]>rowSum[k]) || (rowSum[n]==rowSum[k] && s0WhichCol[n]<s0WhichCol[k]))
	      k=n;
	  }
	  if (k != m) {
	    temp=rowSum[m];
	    rowSum[m]=rowSum[k];
	    rowSum[k]=temp;

	    temp=initialPermute[m];
	    initialPermute[m]=initialPermute[k];	
	    initialPermute[k]=temp;

	    temp=s0WhichCol[m];
	    s0WhichCol[m]=s0WhichCol[k];	
	    s0WhichCol[k]=temp;	  

	    for (kk=0; kk<NumCol; kk++) {
	      temp1=table[m][kk];
	      table[m][kk]=table[k][kk];
	      table[k][kk]=temp1;

	      temp1=s0[m][kk];
	      s0[m][kk]=s0[k][kk];
	      s0[k][kk]=temp1;
	    }
	  }
	}	

	for (n=0; n<NumCol; n++) {
	  s0WhichRow[n]=NumRow;
	  for (m=0; m<NumRow; m++) {
	    if (s0[m][n]==1) {
	      s0WhichRow[n]=m;
	      break;
	    }
	  }
	}

	partialColSum[NumCol-1]=colSum[NumCol-1];
	for (n=NumCol-2; n>=0; n--) {
	  partialColSum[n] = partialColSum[n+1]+colSum[n];
	}
}

static void  InitialSetting (double *pSampleMean, double *pSampleVar, 
			     double *pPValue)
{    
	    *pSampleMean=0;
	    *pSampleVar=0;
	    *pPValue=0;
}

static void GenerateATable (int **t, int **s0, int NumRow, int NumCol, 
			    int *rowSum, int *colSum, int *partialColSum, 
			    int *s0RowSum, int *s0WhichCol, int *s0WhichRow, 
			    int *initialPermute, int *initialColPermute, 
			    double *pCurrentSample)
{
  int i, j, k, m, n, diff, numSplit, partialSumOfRow, tmpMax, temp, temp2, temp3, mm, jj, kk, nn, s0Here, sumOne, minOne, maxOne;
  int *permute, *colPermute, *currentS0RowSum, *currentRowSum, 
    *currentS0WhichCol, *currentS0WhichRow, *rowSegment, *leftoverCol, 
    *space, *oneIndex, *numChoice, *oneWay, *tmp;
  int **currentS0,  **splitPoint; 
  double currentSample, prob;

	currentS0=malloc(NumRow*sizeof(int *));
	currentRowSum=malloc(NumRow*sizeof(int));	
	currentS0RowSum=malloc(NumRow*sizeof(int));
	rowSegment=malloc(NumRow*sizeof(int));
	leftoverCol=malloc(NumRow*sizeof(int));
	oneIndex=malloc(NumRow*sizeof(int));
	tmp=malloc(NumRow*sizeof(int));
	splitPoint=malloc((NumRow+1)*sizeof(int *));
	splitPoint[NumRow]=malloc(2*sizeof(int));
	oneWay=malloc(NumRow*sizeof(int));
	numChoice=malloc(NumRow*sizeof(int));
	space=malloc(NumRow*sizeof(int));
	permute=malloc(NumRow*sizeof(int));
	colPermute=malloc(NumCol*sizeof(int));
	currentS0WhichCol=malloc(NumRow*sizeof(int));
	currentS0WhichRow=malloc(NumCol*sizeof(int));
        for (i=0; i<NumRow; i++) {  
	  splitPoint[i]=malloc(2*sizeof(int));
	  currentS0[i]=malloc(NumCol*sizeof(int));
        }	

	      for (m=0; m< NumRow; m++) {
		permute[m]=initialPermute[m];
	      }
	      for (m=0; m< NumCol; m++) {
		colPermute[m]=initialColPermute[m];
	      }

	      for (j=0;j<NumRow; j++) {
		currentS0RowSum[j]=s0RowSum[j];
		currentRowSum[j]=rowSum[j];
		currentS0WhichCol[j]=s0WhichCol[j];

		for (m=0; m<NumCol; m++) {
		  t[j][m]=0;
		  currentS0[j][m]=s0[j][m];
		}
	      }

	      for (m=0; m<NumCol; m++) {
		currentS0WhichRow[m]=s0WhichRow[m];
	      }     

	      currentSample=0;	      

	      for (j=0; j<NumCol-1; j++) {
		if (colSum[j]>0) {

		  diff=0;
		  numSplit=0;
		  partialSumOfRow=0;

		  for (k=0; k<NumRow-1; k++) {

		    partialSumOfRow = partialSumOfRow + currentRowSum[k]; 

		    if (currentRowSum[k]>0) { 

		      tmpMax=999999;
		      temp=partialColSum[j+1];;

		      if (tmpMax>temp) tmpMax=temp;

		      for (mm=j+1; mm<NumCol; mm++) {
			for (nn=0; nn<=k; nn++) {
			  temp=temp+1-currentS0[nn][mm];
			}
			temp=temp-colSum[mm];
			if (tmpMax>temp) tmpMax=temp;
		      }

		      diff = tmpMax - partialSumOfRow;

		      if (diff<0) {

			if (currentS0WhichRow[j]>k && currentS0WhichRow[j]<NumRow)
			  temp2=1;
			else temp2=0;

			if ((colSum[j]-(NumRow-k-1-temp2) < (-diff)) && (!((numSplit>0) && ((-diff)<=splitPoint[numSplit-1][1])))) {
			  numSplit++;

			  if ((numSplit>1) && currentS0WhichRow[j] >
			      splitPoint[numSplit-2][0] && currentS0WhichRow[j]
			      <= k )
			    temp2=1;
			  else temp2=0;		

			  while ((numSplit>1) && ((-diff) >=
						  k-splitPoint[numSplit-2][0] -temp2+ splitPoint[numSplit-2][1])) {
			    numSplit--;
			    if ((numSplit>1) && currentS0WhichRow[j] >
				splitPoint[numSplit-2][0] &&
				currentS0WhichRow[j] <= k )
			      temp2=1;
			    else temp2=0;		

			  }
			  splitPoint[numSplit-1][0]=k;
			  splitPoint[numSplit-1][1]=-diff;
			}
		      }
		    }
		  }

		  splitPoint[numSplit][0]=NumRow-1;
		  splitPoint[numSplit][1]=colSum[j];

		  // the case with no split
		  if (numSplit==0) {	

		    if (currentS0WhichRow[j] == NumRow) {
		      for (mm=0; mm<NumRow; mm++) {
			leftoverCol[mm] = NumCol-j-currentS0RowSum[mm];
		      }
		      prob=Drafting(currentRowSum,0,NumRow-1,colSum[j],oneIndex, leftoverCol, NumRow);
		      currentSample=currentSample-log(prob);  
		    }
		    else {
		      for (mm=0; mm<currentS0WhichRow[j]; mm++) {
			rowSegment[mm]=currentRowSum[mm];
			leftoverCol[mm] = NumCol-j-currentS0RowSum[mm];
		      }
		      for (mm=currentS0WhichRow[j]+1; mm<NumRow; mm++) {
			rowSegment[mm-1]=currentRowSum[mm];
			leftoverCol[mm-1] = NumCol-j-currentS0RowSum[mm];
		      }
		      prob=Drafting(rowSegment,0,NumRow-2,colSum[j],oneIndex, leftoverCol, NumRow);
		      currentSample=currentSample-log(prob);  
		      for (mm=0; mm<colSum[j]; mm++) {
			if (oneIndex[mm] >= currentS0WhichRow[j]) oneIndex[mm]++;
		      }
		    }
		  }

		  // allocate ones at split points.
		  else {
		    s0Here=numSplit+1;
		    sumOne=0;
		    minOne=splitPoint[0][1];

		    if (currentS0WhichRow[j] <= splitPoint[0][0]) {
		      temp2=1;
		      s0Here=0;
		    }
		    else temp2=0;

		    space[0]=splitPoint[0][0]+1-temp2;
		    maxOne=(space[0]<colSum[j])? space[0]: colSum[j];
		    numChoice[0]=maxOne-minOne+1;
		    oneWay[0]=RandomInteger(minOne, maxOne);
		    sumOne=sumOne+oneWay[0];
		    for (m=1; m<numSplit; m++) {
		      minOne=((splitPoint[m][1]-sumOne)>0)?(splitPoint[m][1]-sumOne) : 0;
		      if (currentS0WhichRow[j] > splitPoint[m-1][0] &&
			  currentS0WhichRow[j] <= splitPoint[m][0]) {
			temp2=1; s0Here=m;}
		      else temp2=0;

		      space[m]=splitPoint[m][0]-splitPoint[m-1][0]-temp2;
		      maxOne=(space[m]<(colSum[j]-sumOne)) ? space[m]: (colSum[j]-sumOne);
		      numChoice[m]=maxOne-minOne+1;
		      oneWay[m] = RandomInteger(minOne,maxOne);
		      sumOne=sumOne+oneWay[m];
		    }
		    numChoice[numSplit]=1;
		    oneWay[numSplit]=colSum[j]-sumOne;

		    if (currentS0WhichRow[j] > splitPoint[numSplit-1][0] &&
			currentS0WhichRow[j] <= NumRow-1) {
		      temp2=1; s0Here=numSplit;}
		    else temp2=0;

		    space[numSplit]=NumRow-splitPoint[numSplit-1][0]-1-temp2;

		    // start sampling
		    sumOne=0;
		    if (oneWay[0]>0) {
		      if (currentS0WhichRow[j]<=splitPoint[0][0]) {
			for (jj=0; jj<currentS0WhichRow[j]; jj++) {					rowSegment[jj]=currentRowSum[jj];
			 leftoverCol[jj] =
			    NumCol-j-currentS0RowSum[jj];
			}		      

			for (jj=currentS0WhichRow[j]+1;
			   jj<=splitPoint[0][0]; jj++) {	
			  rowSegment[jj-1]=currentRowSum[jj];
			  leftoverCol[jj-1] =
			    NumCol-j-currentS0RowSum[jj];
			}

			prob=Drafting(rowSegment,0,splitPoint[0][0]-1,
				      oneWay[0], tmp, leftoverCol, NumRow);
			for (k=0; k<oneWay[0]; k++) {
			  if (tmp[k]>=currentS0WhichRow[j]) 
			    oneIndex[sumOne]=tmp[k]+1;
			  else oneIndex[sumOne]=tmp[k];
			  sumOne++;
			}
		      }
		      else {
			for (jj=0; jj<=splitPoint[0][0]; jj++) {					rowSegment[jj]=currentRowSum[jj];
			 leftoverCol[jj] =
			    NumCol-j-currentS0RowSum[jj];
			}		
			prob=Drafting(rowSegment,0,splitPoint[0][0],
				      oneWay[0], tmp, leftoverCol, NumRow);
			for (k=0; k<oneWay[0]; k++) {
			  oneIndex[sumOne]=tmp[k];
			  sumOne++;
			}
		      }
		    }
		    else {prob=1;}
		    currentSample=currentSample+log(numChoice[0])-log(prob);

		    // sample sections before s0
		    temp3=(s0Here < numSplit+1)? s0Here : (numSplit+1);
		    for (m=1; m<temp3; m++) {
		      if (oneWay[m]>0) {
			for (jj=splitPoint[m-1][0]+1;
			     jj<=splitPoint[m][0]; jj++) {
			  rowSegment[jj]=currentRowSum[jj];
			  leftoverCol[jj] =
			    NumCol-j-currentS0RowSum[jj];
			}
			prob=Drafting(rowSegment,splitPoint[m-1][0]+1,splitPoint[m][0],oneWay[m],tmp, leftoverCol, NumRow);
		      }
		      else {prob=1;}
		      currentSample=currentSample+log(numChoice[m])-log(prob);

		      for (k=0; k<oneWay[m]; k++) {
			oneIndex[sumOne]=tmp[k];
			sumOne++;
		      }
		    }

		    // sample s0 section
		    if (s0Here>0 && s0Here<=numSplit) {
		      m=s0Here;
		      if (oneWay[m]>0) {
			for (jj=splitPoint[m-1][0]+1;
			     jj<=splitPoint[m][0]; jj++) {
			  if (jj < currentS0WhichRow[j]) {
			    rowSegment[jj]=currentRowSum[jj];
			    leftoverCol[jj] =
			      NumCol-j-currentS0RowSum[jj];
			  }
			  else if (jj > currentS0WhichRow[j]) {
			    rowSegment[jj-1]=currentRowSum[jj];
			    leftoverCol[jj-1] =
			      NumCol-j-currentS0RowSum[jj];
			  }
			}
			prob=Drafting(rowSegment,splitPoint[m-1][0]+1,splitPoint[m][0]-1,oneWay[m],tmp, leftoverCol, NumRow);
		      }
		      else {prob=1;}
		      currentSample=currentSample+log(numChoice[m])-log(prob);

		      for (k=0; k<oneWay[m]; k++) {
			if (tmp[k]>=currentS0WhichRow[j]) 
			  oneIndex[sumOne]=tmp[k]+1;
			else oneIndex[sumOne]=tmp[k];
			sumOne++;
		      }
		    }

		    //sample sections after s0
		    for (m=s0Here+1; m<=numSplit; m++) {
		      if (oneWay[m]>0) {
			for (jj=splitPoint[m-1][0]+1;
			     jj<=splitPoint[m][0]; jj++) {
			  rowSegment[jj]=currentRowSum[jj];
			  leftoverCol[jj] =
			    NumCol-j-currentS0RowSum[jj];
			}
			prob=Drafting(rowSegment,splitPoint[m-1][0]+1,splitPoint[m][0],oneWay[m],tmp, leftoverCol, NumRow);
		      }
		      else {prob=1;}
		      currentSample=currentSample+log(numChoice[m])-log(prob);

		      for (k=0; k<oneWay[m]; k++) {
			oneIndex[sumOne]=tmp[k];
			sumOne++;
		      }
		    }
		  }

		  //update after sampling
		  for (m=0; m<colSum[j]; m++) {
		    currentRowSum[oneIndex[m]]--;
		    t[oneIndex[m]][j]=1;
		  }

		  if (currentS0WhichRow[j]<NumRow) {
		    currentS0RowSum[currentS0WhichRow[j]]--;
		    currentS0WhichCol[currentS0WhichRow[j]]=NumCol;
		  }

		  for (m=0; m<NumRow; m++) {
		    k=m;
		    for (n=m+1; n<NumRow; n++) {
		      if ((currentRowSum[n]>currentRowSum[k]) || (currentRowSum[n]==currentRowSum[k] && currentS0WhichCol[n]<currentS0WhichCol[k]))
			k=n;
		    }
		    if (k != m) {
		      temp=currentRowSum[m];
		      currentRowSum[m]=currentRowSum[k];
		      currentRowSum[k]=temp;

		      temp=currentS0RowSum[m];
		      currentS0RowSum[m]=currentS0RowSum[k];
		      currentS0RowSum[k]=temp;

		      if (currentS0WhichCol[m]<NumCol) currentS0WhichRow[currentS0WhichCol[m]]=k;
		      if (currentS0WhichCol[k]<NumCol) currentS0WhichRow[currentS0WhichCol[k]]=m;

		      temp=currentS0WhichCol[m];
		      currentS0WhichCol[m]=currentS0WhichCol[k];
		      currentS0WhichCol[k]=temp;

		      temp=permute[m];
		      permute[m]=permute[k];
		      permute[k]=temp;	

		      for (kk=0; kk<=j; kk++) {
			temp=t[m][kk];
			t[m][kk]=t[k][kk];
			t[k][kk]=temp;
		      }

		      for (kk=j+1; kk<NumCol; kk++) {
			temp=currentS0[m][kk];
			currentS0[m][kk]=currentS0[k][kk];
			currentS0[k][kk]=temp;
		      }
		    }
		  }
		}
		else {
		  for (kk=0; kk<NumRow; kk++) 
		    t[kk][j]=0;
		}
	      }

	      for (k=0; k<NumRow; k++) {
		t[k][NumCol-1]=currentRowSum[k];
	      }

	      currentSample=exp(currentSample);

	      for (m=0; m<NumRow-1; m++) {
		while (permute[m] != m) {
		  for (kk=0; kk<NumCol; kk++) {
		    temp=t[m][kk];
		    t[m][kk]=t[permute[m]][kk];
		    t[permute[m]][kk]=temp;
		  }
		  temp=permute[m];
		  permute[m]=permute[permute[m]];
		  permute[temp]=temp;
		}
	      }				  

	      for (m=0; m<NumCol-1; m++) {
		while (colPermute[m] != m) {
		  for (kk=0; kk<NumRow; kk++) {
		    temp=t[kk][m];
		    t[kk][m]=t[kk][colPermute[m]];
		    t[kk][colPermute[m]]=temp;
		  }
		  temp=colPermute[m];
		  colPermute[m]=colPermute[colPermute[m]];
		  colPermute[temp]=temp;		  
		}
	      }
	      *pCurrentSample=currentSample;

	      for (i=0; i<NumRow; i++) {  
		free(splitPoint[i]);
		free(currentS0[i]);
	      }	
	      free(currentS0);
	      free(currentRowSum);
	      free(currentS0RowSum);
	      free(rowSegment);
	      free(leftoverCol);
	      free(oneIndex);
	      free(tmp);
	      free(splitPoint);
	      free(splitPoint[NumRow]);
	      free(oneWay);
	      free(numChoice);
	      free(space);
	      free(permute);
	      free(colPermute);
	      free(currentS0WhichCol);
	      free(currentS0WhichRow);
}	            

static void UpdatePValue (double currentSample, int **t, 
			  int NumCol, double reciprocal0, 
			  double *pSampleMean, 
		          double *pSampleVar, double *pPValue) 
{
	      (*pSampleMean) = (*pSampleMean)+currentSample;
	      (*pSampleVar)=(*pSampleVar)+currentSample*currentSample;

	      if (Reciprocal(t, NumCol) >= reciprocal0)
		//    if (Reciprocal(t, NumCol) >= 18)	
		(*pPValue)=(*pPValue)+currentSample;	      
}

static void PrintOutput (FILE *outFile, int N, double sampleMean, 
			 double sampleVar, double pValue)
{	    
            double cvsquare;

	    cvsquare=N*sampleVar/sampleMean/sampleMean-1;	    
	    pValue=pValue/sampleMean;
	    sampleMean=sampleMean/N;

            fprintf(outFile, "%20g %20g %20g\n", sampleMean, pValue, cvsquare);
	    printf("%20g %20g %20g\n", sampleMean, pValue, cvsquare);
}

static double  NChooseM(int n, int m)
{
	int i;
	double value;

	value=1.0;

	if ((m<=0) || (n<=0) || (n<m)) return(value);
	for (i=1; i<=m; i++) {
		value=value*(n-i+1)/i;
	}
	return (value);

}

static double Mean(double x[], int n)
{
	int i;
	double sum=0;

	for(i=0; i<n; i++) sum = sum+x[i];
	return (sum/n);
}

static int RandomSample2(double weight[], int n)
{
	int low, up, mid;
	double randomNum;

	randomNum=RandomReal(0, weight[n-1]);
	low=0;
	up=n-1;
	mid=n/2;

	while (1) {
		if (randomNum>=weight[mid]) {
			low=mid+1;
			mid=(mid+up)/2;
			if (low==up) {
				mid=up;
				break;
			}
		}
		else if ((mid>0) && (randomNum<weight[mid-1])) {
			up=mid-1;
			mid=(mid+low)/2;
		}
		else break;
	}
	return (mid);
}

static double Drafting(int w[], int start, int end, int n, int a[],
		       int s[], int NumRow)
{
        double *p, *newP;
	int i,j;
	double tmp;
	double prob;
	double *currentWeight;
	double r1, r2;
	double *weight;

	p=malloc(NumRow*sizeof(double));
	newP=malloc(NumRow*sizeof(double));
	currentWeight=malloc(NumRow*sizeof(double));
	weight=malloc(NumRow*sizeof(double));

	if (n==0) return (1);

	if (end-start+1==n) {
	  for (i=start; i<=end; i++) {
	    a[i-start]=i;
	  }
	  return(1);
	}

	for (i=start; i<=end; i++) {
		if (w[i]==0) {
			weight[i]=0;
		}
		else {
			weight[i]=exp(log(1.0*w[i]/(s[i]-w[i])));
		}
	}

	for (i=start; i<=end; i++) {
		tmp=weight[i];
		weight[i]=weight[end];
		p[i]=tmp*R(n-1, weight, start, end-1, NumRow);
		weight[i]=tmp;
		currentWeight[i]=weight[i];
	}

	a[0]=ChooseOneInteger(p, start, end, NumRow);
	currentWeight[a[0]]=0;
	j=0;

	while (n>(j+1)) {
		for (i=start; i<=end; i++) {
		  if (currentWeight[i]<=0.00000000001)

				newP[i]=0;
		  else if (abs(weight[i]-weight[a[j]])<0.0000000001) {

				currentWeight[i]=0;
				r1=R(n-j-2, currentWeight, start, end, NumRow);
				currentWeight[i]=weight[i];
				r2=R(n-j-1, currentWeight, start, end, NumRow);
				newP[i]=weight[i]*r1/(n-j-1)/r2;
			}
			else {
				newP[i]=(weight[a[j]]*p[i]-weight[i]*p[a[j]])/(n-j-1)/(weight[a[j]]-weight[i])/p[a[j]];
			}
		}
		j++;

		a[j]=ChooseOneInteger(newP, start, end, NumRow);
		currentWeight[a[j]]=0;

		for (i=start; i<=end; i++) {
			p[i]=newP[i];
		}
	}

	prob=1;
	for (i=0; i<n; i++) {
		prob=prob*weight[a[i]];
	}
	prob=prob/R(n, weight, start, end, NumRow);

	free(p);
	free(newP);
	free(currentWeight);
	free(weight);

	return(prob);
}

static int ChooseOneInteger(double weight[], int start, int end, int NumRow)
{
	int i;
	double *totalWeight, temp;

	totalWeight=malloc(NumRow*sizeof(double));

	totalWeight[start]=weight[start];
	for (i=start+1; i<=end; i++) {
		totalWeight[i]=totalWeight[i-1]+weight[i];
	}

	temp=RandomReal(0,totalWeight[end]);

	i=start;
	while (totalWeight[i]<=temp) {
		i++;
	}
	free(totalWeight);
	return (i);
}

static double R(int k, double w[], int start, int end, int NumRow)
{
	int i, j;
	double *t; 
	double *r;

	t=malloc((NumRow+1)*sizeof(double));
	r=malloc(NumRow*sizeof(double));

	for (i=1; i<=k; i++) {
		t[i]=T(i, w, start, end);
	}

	r[0]=1;
	for (j=1; j<=k; j++) {
		r[j]=0;
		for (i=1; i<=j; i++) {
			if (i%2==0) 
				r[j]=r[j]-t[i]*r[j-i];
			else 
				r[j]=r[j]+t[i]*r[j-i];
		}
		r[j]=r[j]/j;
	}
	free(t);
	return (r[k]);
}

static double T(int i, double w[], int start, int end) 
{
	double sum;
	int j, k;
	double tmp;

	sum=0;
	for (j=start; j<=end; j++) {
	  tmp=1.0;
	  for (k=0; k<i; k++) {
	    tmp=tmp*w[j];
	  }
	  sum = sum+tmp;
	}
	return (sum);
}

static int RandomInteger(int low, int high)
{
  int k;
  double d;

  d=(double) rand()/((double) RAND_MAX+1);
  k=(int) (d*(high-low+1));
  return (low+k);
}

static double RandomReal(double low, double high)
{
       double d;

       d=(double) rand()/((double) RAND_MAX+1);
       return (low+d*(high-low));
}