#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, ¤tSample); 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)); }