#--------------------------------CUT HERE------------------------------------- #! /bin/sh # # This is a shell archive. Save this into a file, edit it # and delete all lines above this comment. Then give this # file to sh by executing the command 'sh file'. The files # will be extracted into the current directory owned by # you with default permissions. echo 'x - README' sed 's/^X//' << '________This_Is_The_END________' > README X XThis file contains instructions for using the software we have written Xto compute the approximations and bounds described in our paper: X XHuffer, F.W. and Lin, C. T. (1997). X"Approximating the distribution of the scan statistic using moments of Xthe number of clumps", Journal of the American Statistical XAssociation, Vol 92, pp. 1466--1475. X XIn what follows, we refer by number to the tables in this paper. XAlso, we use the variable names and acronyms from this paper. X XQuestions concerning the software can be sent to the authors: X huffer@stat.fsu.edu, chien@math.tku.edu.tw X XThis file gives a brief explanation of how to use our programs. The Xcomments in the files maplecode and moment_four.c contain more Xdetailed information on many points. X XThese programs may be freely distributed for non-commercial purposes. X XSome related programs are posted at the web address: X http://stat.fsu.edu/~huffer X X=============================================================== X X===================== XUsing the C programs X==================== X XThe C program "moment_four.c" produces Maple expressions for the first Xfour moments of the number of m:d clumps Y. On our local system we Xcompile this program using the command X X cc moment_four.c -lm X Xwhich creates an executable file named a.out. (The program must be Xcompiled using ANSI standard C, not the old C. You may need to add an Xoption to the above command to ensure that ANSI C is used.) You now Xexecute this program with two command line arguments. The first Xargument is the clump size m, an integer from 3 to 10. The second Xargument is the name of the output file. The output file will contain Xthe expressions for the moments. For example, executing X X a.out 7 fourmom.7 X Xcreates the output file fourmom.7. The output files MUST be given Xnames of the form fourmom.m because the maple procedure readmom Xmentioned below expects the output files to follow this naming Xconvention. When you execute readmom(m), it will look for a file Xnamed fourmom.m. X XThe C program moment_two.c is similar to moment_four.c except that Xit only gives expressions for the first two moments. It can handle Xvalues of m up to 22. X XFor a given value of the clump size m, before you can use the Maple Xprocedure readmom(m), you must first use the C program to create the Xfile fourmom.m. X X===================== XUsing the Maple code X===================== X XWe assume that all the files we need are in the current Xdirectory. Minor changes must be made in the code if the files Xare somewhere else. X XThe file maplecode contains the maple commands for computing Xall of the approximations and bounds. X XAfter starting Maple, type this READ command: X X read `maplecode`; X XThis must be done before any of the following. X X------------------------ XCreating Tables 1 and 3. X------------------------ X XTo compute the results in our Tables 1 and 3 you first choose Xvalues for m, d, n, and then execute the following commands X(with m, d, n replaced by numerical values). X X readmom(m); [Note: produces no output] X setval(d,n); [Note: produces no output] X XNow to get the results in Table 1 execute: X X mc2(); lp4(0); cpg4(1); klb(); bounds(); X XTo get the results in Table 3 execute: X X cpg2(5); cpg4(5); gcp(5); X XNotes: X X1. In order to use any of the procedures lp4, cpg4 or bounds, you MUST Xfirst run the readmom and setval procedures. X X2. Executing readmom(m) reads in the expressions for the first four Xmoments for the given clump size m. These expressions must already Xexist and be stored in a file named fourmom.m (where m is an integer Xbetween 3 and 10). To create these expressions you run the C program Xmoment_four.c. X X3. The procedures mc2 and cpg2 require only the first two moments Xwhich you can compute using the fsmval procedure (described below). Xgcp is more self-contained and does not even require fsmval. X X4. The argument k in cpg2(k), cpg4(k), gcp(k) produces approximations Xfor P(Y >= j) for j=1,2,...,k. X X5. Executing lp4(1) gives additional output (the lambda's) not Xsupplied by lp4(0). X X6. If you want to do more calculations with different values of d and Xn, but the SAME value of m, you need only re-execute setval(d,n) with Xthe new values. If you change the value of m, you must re-execute XBOTH readmom(m) and setval(d,n) even if you are not changing d or n. X X7. WARNING: The variables n, m, d, w, m1, m2, m3, m4, c1, c2, c3, c4, Xmom1, mom2, mom3, mom4 are used as global variables by the programs. XThe user should not use these variable names for other purposes. X X8. The file maplecode also contains procedures ap1, ap2, ap3, ap4 for Xcomputing the approximations with those names from Huffer and Lin X(1995). (This is the technical report version of our paper and is Xavailable from the authors.) X X---------------- XCreating Table 2 X---------------- X XTo compute the results in our Table 2, you choose values Xfor m, d, n, and then execute the following commands X(with m, d, n replaced by numerical values). X X fsmval(m,d,n); [Note: produces no output] X cpg2(1); mc2(); X XNotes: X X1. The procedure fsmval computes the first two moments of Y. It Xcan be used for large values of m, but becomes progressively Xslower as m increases. X X2. fsmval does NOT require running the C program before it is Xused (whereas readmom and setval do require it). X X3. cpg2(k) will produce approximate values for P( Y >= j), Xj=1,...,k. X ________This_Is_The_END________ if test `wc -l < README` -ne 154; then echo 'shar: README was damaged during transit (should have had 154 lines)' fi echo 'x - moment_four.c' sed 's/^X//' << '________This_Is_The_END________' > moment_four.c X/* X moment_four.c 8-5-1997 X ---------------------- X This program prints out the expressions of the first X four moments for different value of m. X X Execute the program as the followings. X gcc -o 'execute file' moment_four.c -lm X 'execute file' m 'output file' X*/ X X#include X#include X#include X Xtypedef struct cnode{ X int j, /* the values of b(j,k). */ X k; X double a; /* coefficients of the weight. */ X struct cnode *N,*P; X}LIST; X Xtypedef struct lnode{ X int **mat; /* DE matrix */ X int **se; /* starting and ending index matrix */ X int row, /* rows of the matrix. */ X num, /* 1's in side the matrix */ X bk; /* the number of block inside the matrix. */ X LIST *C; X struct lnode *N,*P; X}LINK; X Xtypedef struct tnode{ X int *e; X int row; /* dimension of the matrix. */ X LIST *C; X struct tnode *L,*R; X}TREE; X Xtypedef struct{ /* comparing the blocks inside the matrix */ X int **t; /* the submatirx of the block */ X int len; /* length of the block */ X}com; X Xtypedef struct pnode{ X int s,t; /* the values in R(s,t) */ X LIST *C; X struct pnode *L,*R; X}PP; X XFILE *fp; XLIST *W; XTREE *S; XPP *R; Xint clump; Xint **New_Array(int,int); XLINK *Lalloc(void); XTREE *Talloc(void); XLIST *Calloc(void); X Xmain(argc, argv) Xint argc; Xchar *argv[]; X{ X double ca,tca; X int **nmat; X int m,I,M,i,j,k,ii,jj,incre; X int colnum,cj,ck; X LINK *Head; X LIST *c; X double acoeff(int,int); X LIST *listfree(LIST *); X LINK *translate(LINK *,int **,double,int,int,int,int); X LINK *decom(LINK *); X LIST *bcoef(LIST *); X PP *Pfree(PP *); X TREE *Delete(TREE *); X void print_tree(TREE *); X void polyprint(PP *); X X X X if (argc != 3) X { X fprintf(stderr, "Usage: %s m outfilename\n", argv[0]); X return 1; X } X X X if((fp=fopen(argv[2],"a+"))==NULL){ /* data file not found */ X fprintf(stderr, "Cannot open data file %s\n", argv[2]); X exit(1); X } X X m = atoi(argv[1]); X X clump=m-1; X X X fprintf(fp,"m:=%d;\n",m); X X for(M=4;M>=1;M--){ X Head=NULL; X S=NULL; X R=NULL; X W=NULL; X if(M==1) X fprintf(fp,"m1"); X else if(M==2) X fprintf(fp,"m2"); X else if(M==3) X fprintf(fp,"m3"); X else if(M==4) X fprintf(fp,"m4"); X fprintf(fp,":='\n"); X for(I=1;I<=M;I++){ X tca=acoeff(M,I); X if(I==1){ X cj=m-1; X ck=1; X ca=tca; X colnum=clump; X nmat=New_Array(I,colnum); X for(i=0;iN){ X if(c->a>0) X fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k); X else if(c->a<0) X fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k); X } X fprintf(fp,")\n"); X X W=listfree(W); X X polyprint(R); X X R=Pfree(R); X fprintf(fp,"';\n"); X } X fclose(fp); X} X X X X/* New_Array: make an allocation for a matrix. */ Xint **New_Array(int row, int col) X{ X int **A,i,j; X X A=(int **)malloc(row*sizeof(int*)); X if(A==NULL) X return NULL; X for(i=0;i a) X return 0.0; X else if (a<=0 || b<=0) X return 0.0; X else X return b*acoeff(a-1,b)+b*acoeff(a-1,b-1); X} X X XLINK *translate(LINK *rt,int **d,double ca,int cj,int ck,int row,int col) X{ X int asize,i,j,k,ii,jj; X int **a,**b,*delcl,*indices,**temp; X int star,location,jump,count,ncol; X LINK *NEW; X LIST *cnew; X double Comb(double,double); X int **block(LINK *); X LINK *createlink(LINK *,LINK *); X LIST *combine(LIST *,LIST *); X void permute(int,int,int *,int,int **,int *); X X X cnew=Calloc(); X if(cnew==NULL){ X fprintf(fp,"no room for temp(cplist)\n"); X exit(0); X } X cnew->j=cj; X cnew->k=ck; X cnew->a=ca; X cnew->N=cnew->P=NULL; X if(W==NULL) X W=cnew; X else X W=combine(W,cnew); X X delcl=(int *)malloc(col*sizeof(int)); X if(delcl==NULL){ X fprintf(fp,"No room for delcl(translate)\n"); X exit(0); X } X indices=(int *)malloc(row*sizeof(int)); X if(indices==NULL){ X fprintf(fp,"No room for indices(translate)\n"); X exit(0); X } X for(k=1;k<=row;k++){ X ca *= -1.0; X location=0; X star=k; X asize=Comb(row,star); X /* This is achievable job since the maximum row is 4.*/ X a=New_Array(asize,row); X if(a==NULL){ X fprintf(fp,"No room for a(translate)\n"); X exit(0); X } X for(i=0;imat=b; /* cp the DE matrix to node */ X NEW->row=k; X NEW->num=0; /* count 1's inside the matrix */ X for(i=0;irow;i++) X NEW->num += NEW->mat[i][1]; X NEW->bk=0; /* initialize the block be 0 */ X NEW->N=NEW->P=NULL; X NEW->se=block(NEW); /* examine the blocks inside the matrix */ X NEW->C=Calloc(); X if(NEW->C==NULL){ X fprintf(fp,"No room for NEW->C(translate)\n"); X exit(0); X } X NEW->C->j=cj; X NEW->C->k=ck; X NEW->C->a=ca; X NEW->C->N=NEW->C->P=NULL; X rt=createlink(rt,NEW); X } X for(i=0;i= m) X { X for (i = 0; i < n; i++) X { X int j; X int value = 0; /* first assumes the value as 1 */ X X /* if i is not in the indIndices, then it's 0 */ X for (j = 0; !value && j < whichIndex; j++) X { X value = i == indIndices[j]; X } X a[*location][i]=value; X if(i==n-1) X (*location)++; X } X return; X } X for (i = whichIndex == 0 ? 0 : indIndices[whichIndex - 1] + 1; X i <= n - (m - whichIndex); X i++) X { X indIndices[whichIndex] = i; X permute(n, m, indIndices, whichIndex + 1,a,location); X } X return; X} X Xdouble Fact(double n) X{ X if(n==0.0) X return 1.0; X else X return n*Fact(n-1.0); X} X Xdouble Comb(double a, double b) X{ X int i,ans=1; X double Fact(double); X X if(a < b) X return 0.0; X else{ X for(i=b;i>=1.0;i--,a--) X ans=ans*a; X ans=ans/Fact(b); X return ans; X } X} X X/* createlink: add the node to the linklist */ XLINK *createlink(LINK *rt,LINK *NEW) X{ X int i,j,match; X LINK *cur,*prev; X LINK *addchecking(LINK *,LINK *,LINK *,LINK *); X LIST *combine(LIST *,LIST *); X TREE *addtree(TREE *,LINK *); X X X/* check if the new matrix is simple matrix or not */ X if(NEW->bk==0){ X /*simple matrix call the shorthand fun. */ X S=addtree(S,NEW); X return rt; X }else{ X if(rt==NULL){ X return NEW; X }else{ X for(prev=NULL,cur=rt;cur!=NULL;prev=cur,cur=cur->N){ X /* larger size of matrix first be decomposed */ X if(NEW->rowrow) X continue; X else if(NEW->row>cur->row) X return addchecking(rt,prev,cur,NEW); X else if(NEW->row==cur->row){ X /* more 1's inside the matrix first be decomposed */ X if(NEW->numnum) X continue; X else if(NEW->num>cur->num) X return addchecking(rt,prev,cur,NEW); X else if(NEW->num==cur->num){ X /* larger size of block is decomposed first */ X if(NEW->se[NEW->bk-1][2]se[cur->bk-1][2]) X continue; X else if(NEW->se[NEW->bk-1][2]>cur->se[cur->bk-1][2]) X return addchecking(rt,prev,cur,NEW); X else if(NEW->se[NEW->bk-1][2]==cur->se[cur->bk-1][2]){ X /*check if the identical term */ X for(i=0,match=1;irow && match==1;i++){ X for(j=0;j<3 && match==1;j++){ X if(NEW->mat[i][j]!=cur->mat[i][j]) X match=0; X } X } X if(match==0) /* Not identical, move to next node to check. */ X continue; X else{ /* Identical, add coeffs to LIST */ X cur->C=combine(cur->C,NEW->C); X for(i=0;ibk;i++) /* free the node */ X free(NEW->se[i]); X free(NEW->se); X for(i=0;irow;i++) X free(NEW->mat[i]); X free(NEW->mat); X free(NEW); X return rt; X } X } X } X } X } X /* end of for-loop */ X /*connect the new node to the end of list */ X prev->N=NEW; X NEW->P=prev; X return rt; X } /* end of else-when rt!=NULL */ X } X} X X/* addchecking: add the node into the appropriate place */ XLINK *addchecking(LINK *rt,LINK *prev,LINK *cur,LINK *NEW) X{ X if(prev==NULL){ X NEW->N=cur; X cur->P=NEW; X return NEW; X }else{ X NEW->N=cur; X NEW->P=prev; X prev->N=NEW; X cur->P=NEW; X return rt; X } X} X X/* combine: combine the two list */ XLIST *combine(LIST *old,LIST *new) X{ X int jp; X LIST *p,*c,*temp; X LIST *listcheck(LIST *,LIST *,LIST *,LIST *); X X /* old has at least one node in the LIST */ X while(new!=NULL){ X temp=new; X new=new->N; X if(new!=NULL) X new->P=NULL; X temp->N=NULL; X for(p=NULL,c=old,jp=0;c!=NULL;p=c,c=c->N){ X if(temp->kk) X continue; X else if(temp->k>c->k){ X old=listcheck(old,p,c,temp); X jp=1; X break; X }else if(temp->k==c->k){ X if(temp->jj) X continue; X else if(temp->j>c->j){ X old=listcheck(old,p,c,temp); X jp=1; X break; X }else if(temp->j==c->j){ X c->a += temp->a; X jp=1; X free(temp); X break; X } X X } X } X if(jp==0){ X p->N=temp; X temp->P=p; X } X } X return old; X} X X/* listcheck: add the node into the appropriate place */ XLIST *listcheck(LIST *rt,LIST *prev,LIST *cur,LIST *NEW) X{ X if(prev==NULL){ X NEW->N=cur; X cur->P=NEW; X return NEW; X }else{ X NEW->N=cur; X NEW->P=prev; X prev->N=NEW; X cur->P=NEW; X return rt; X } X} X X/* block: check the blocks inside the matrix and return the number of block and matrix */ Xint **block(LINK *b) X{ X X int i,j,k,l,kk,ll,end; X int *bk=&b->bk; /* initialize to be 0 */ X int **temp; X com *c; X int Num_Cmp1( const void * , const void * ); X int Num_Cmp2( const void * , const void * ); X X if(b->row==1) X return NULL; X /* record the matrix where is not overlapped with other rows */ X temp=New_Array(b->row,3); X if(temp==NULL){ X fprintf(fp,"No room for temp(block)\n"); X exit(0); X } X /* X Good for b->row > 2 X When b->mat[i][2]=0 and b->mat[i+1][2]!=0 it starts the block X b->mat[i][2]!=0 and b->mat[i+1][2]=0 it ends the block X b->se[i][0] records the starting row X b->se[i][1] records the ending row X b->se[i][2] records the number of rows inside the block X */ X for(k=0,i=end=1;irow;i++){ X if(b->mat[i][2]==0 && b->mat[i-1][2]==0){ X for(ll=0;ll<3;ll++) /* cp the nonoverlapped row */ X temp[k][ll]=b->mat[i-1][ll]; X k++; X }else if(b->mat[i][2]!=0 && b->mat[i-1][2]!=0){ /* still inside the block */ X continue; X }else if(b->mat[i][2]==0 && b->mat[i-1][2]!=0){ /* it ends the block */ X b->se[*bk][1]=i-1; /* the ending row */ X b->se[*bk][2]=i-b->se[*bk][0]; /* the number of rows inside this block */ X (*bk)++; /* increment to next block */ X end=1; /* index for ending the block */ X }else if(b->mat[i][2]!=0 && b->mat[i-1][2]==0){ /* it starts the block */ X if(*bk==0){ X b->se=New_Array((int)b->row/2,3); /* The maximum blocks will be row/2 */ X if(b->se==NULL){ X fprintf(fp,"No room for b->se(block)\n"); X exit(0); X } X } X b->se[*bk][0]=i-1; /* the staring row of the first block */ X end=0; /* index for starting the block */ X } X if(i==b->row-1 && b->mat[i][2]==0){ /* examine the last row */ X for(ll=0;ll<3;ll++) X temp[k][ll]=b->mat[i][ll]; X k++; X } X } X if(end==0){ /* completes the record of the block */ X b->se[*bk][1]=b->row-1; X b->se[*bk][2]=b->row-b->se[*bk][0]; X (*bk)++; X } X /* if it is simple, sort(increasing) the number of 1's among rows */ X if(*bk==0 && k==b->row){ X qsort(b->mat,b->row,sizeof(int *),Num_Cmp1); X for(i=0;irow;i++){ X if(i==0) X b->mat[i][0]=0; X else X b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1]; X } X b->se=NULL; X }else{ /* There is at least one block inside the matrix */ X for(i=*bk;i<(int)b->row/2;i++) /* Free the un-used part */ X free(b->se[i]); X /* if k!=0 then put those unoverlapped rows to the front and sort them X Also compare the size of the blocks and sort them at the back */ X if(k!=0){ X if(b->se[0][0]!=k){ X for(i=0,kk=k;i<*bk;i++){ X l=kk; X for(j=b->se[i][0];j<=b->se[i][1];j++){ X for(ll=0;ll<3;ll++) X temp[kk][ll]=b->mat[j][ll]; X kk++; X } X b->se[i][0]=l; X b->se[i][1]=kk-1; X } X for(i=0;irow;i++){ X for(ll=0;ll<3;ll++) X b->mat[i][ll]=temp[i][ll]; X } X qsort(b->mat,b->se[0][0],sizeof(int *),Num_Cmp1); X for(i=0;ise[0][0];i++){ X if(i==0) X b->mat[i][0]=0; X else X b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1]; X } X } X } X if(*bk>=2){ X c=(com *)malloc(*bk*sizeof(com)); X if(c==NULL){ X fprintf(fp,"No c space(block)\n"); X exit(0); X } X for(i=0;i<*bk;i++){ X c[i].t=New_Array(b->se[i][2],3); X if(c[i].t==NULL){ X fprintf(fp,"No c.t space(block)\n"); X exit(0); X } X for(j=0,kk=b->se[i][0];kk<=b->se[i][1];j++,kk++){ X for(ll=0;ll<3;ll++) X c[i].t[j][ll]=b->mat[kk][ll]; X } X c[i].len=b->se[i][2]; X } X qsort(c, *bk, sizeof(com), Num_Cmp2); X for(i=0,j=b->se[0][0];i<*bk;i++){ X for(kk=0;kkmat[j][ll]=c[i].t[kk][ll]; X } X /* rearrange the values inside b->se */ X b->se[i][2]=c[i].len; X if(i!=0) X b->se[i][0]=b->se[i-1][1]+1; X else X b->se[0][0]=k; X b->se[i][1]=b->se[i][0]+b->se[i][2]-1; X for(kk=0;kkmat[i][0] */ X X for(i=0,j=b->se[0][0];i<*bk;i++){ X for(kk=0;kkse[i][2];kk++,j++){ X if(kk==0){ X if(i==0 && k==0) X b->mat[0][0]=0; X else X b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1]; X }else X b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1]+b->mat[j][2]-b->mat[j][1]; X } X } X } X for(i=0;irow;i++) X free(temp[i]); X free(temp); X return b->se; X} X Xint Num_Cmp1( const void * S1, const void * S2) X{ X int ** T1 = (int **) S1, X ** T2 = (int **) S2; X X if(*(*T1+1) < *(*T2+1)) X return -1; X else if (*(*T1+1) > *(*T2+1)) X return 1; X else X return 0; X} X Xint Num_Cmp2( const void * S1, const void * S2) X{ X com * T1 = (com *) S1, X * T2 = (com *) S2; X X if(T1->len < T2->len) X return -1; X else if (T1->len > T2->len) X return 1; X else X return 0; X} X X X XTREE *addtree(TREE *p,LINK *I) X{ X int i,match; X LIST *combine(LIST *,LIST *); X X if(p==NULL){ X p=Talloc(); X if(p==NULL){ X fprintf(fp,"No room for p(addtree)\n"); X exit(0); X } X p->e=(int *)malloc(I->row*sizeof(int)); X if(p->e==NULL){ X fprintf(fp,"No room for p->e(addtree)\n"); X exit(0); X } X for(i=0;irow;i++) X p->e[i]=I->mat[i][1]; X p->row=I->row; X p->C=I->C; X p->L=p->R=NULL; X for(i=0;irow;i++) X free(I->mat[i]); X free(I->mat); X /* I->se is NULL, no need to free it */ X free(I); X }else if(I->row>p->row){ /* if row is greater, goes to Left Tree */ X p->L=addtree(p->L,I); X }else if(I->rowrow){ /* if row is lesser, goes to Right Tree */ X p->R=addtree(p->R,I); X }else if(I->row==p->row){ X for(i=match=0;irow;i++){ X if(I->mat[i][1]==p->e[i]){ X continue; X }else if(I->mat[i][1]>p->e[i]){ X match=1; X break; X }else if(I->mat[i][1]e[i]){ X match=2; X break; X } X } X if(match==1){ /* Not identical, insert to Left Tree */ X p->L=addtree(p->L,I); X }else if(match==2){ X p->R=addtree(p->R,I); X }else if(match==0){ X /* combine two link list of coefficients */ X p->C=combine(p->C,I->C); X for(i=0;irow;i++) X free(I->mat[i]); X free(I->mat); X free(I); X } X } X return p; X} X X/* decom: decompose the matrices inside the link list */ XLINK *decom(LINK *p) X{ X int i,j,col,len,tg; X int *loc; X LIST *cr; X int *coef(LINK *,int *,int *,int *); X LINK *organ(LINK *,int,int,double); X LINK *Lfree(LINK *); X X while(p!=NULL){ X /* We only look at the last big block and X put corresponding columns to coeffs inside loc */ X col=p->mat[p->row-1][0]+p->mat[p->row-1][1]; X loc=(int *)malloc(col*sizeof(int)); X if(loc==NULL){ X fprintf(fp,"No room for loc(decom)\n"); X exit(0); X } X /* initialize the values of array loc */ X for(i=0;iN!=NULL){ X p=p->N; X p->P=Lfree(p->P); X }else{ X p=Lfree(p); X } X } X return p; X} X X/* Lfree: To free the node inside the link list */ XLINK *Lfree(LINK *p) X{ X int i; X LIST *cr; X X for(i=0;ibk;i++) X free(p->se[i]); X free(p->se); X for(i=0;irow;i++) X free(p->mat[i]); X free(p->mat); X cr=p->C; X while(cr!=NULL){ X if(cr->N!=NULL){ X cr=cr->N; X free(cr->P); X cr->P=NULL; X }else{ X free(cr); X cr=NULL; X } X } X free(p); X return NULL; X} X X X/* coef: find the location of coefficients and return the taget */ Xint *coef(LINK *p,int *c,int *tar,int *l) X{ X int mr,k,sum,i,j; X int *fill(LINK *,int,int *,int *); X X mr=k=p->row-1; /* start to mark the set at the last row of the block */ X for(sum=0;k>=p->se[p->bk-1][0];k--){ X sum=sum+p->mat[k][2]; X if(sum>=p->mat[mr][1]){ X c=fill(p,mr,c,l); X /* if equal continue to mark set X not equal means stop and tg=0 and return */ X if(sum==p->mat[mr][1]){ X mr=k-1; /* the next marking row and initialize sum=0 */ X sum=0; X }else if (sum>p->mat[mr][1]){ /* break here, no need to continue mark set */ X *tar=0; /* target is 0 vector */ X return c; X } X } X /* if m* is the first row of the block */ X if(mr==p->se[p->bk-1][0]){ X *tar=1; /* target is (1,0,...0) */ X return c; X } X } X /* if m* is not the first row of the block */ X c=fill(p,mr,c,l); X *tar=0; /* target is 0 vector */ X return c; X} X X/* fill: find the ending and starting columns of current marking row */ Xint *fill(LINK *p,int mr,int *c,int *l) X{ X if(mr==p->row-1){ X c[(*l)++]=p->mat[mr][0]+p->mat[mr][1]-1; /* the last ending column */ X } X c[(*l)++]=p->mat[mr][0]; /* the starting column of current one */ X c[(*l)++]=p->mat[mr][0]-1; /* the ending column of next one */ X return c; X} X X/* organ: organize the new produced matrix X (also delet the redundant rows) and to the link list */ XLINK *organ(LINK *p,int c,int tg,double r) X{ X int i,j,k,kk,l,jp; X int **d,*delrw; X LINK *NEW; X LIST *cur; X LIST *cplist(LIST *,LIST *,double); X X /* d is the replaced matrix accroding to tg */ X d=New_Array(p->row,3); X if(d==NULL){ X fprintf(fp,"No room for d(organ)\n"); X exit(0); X } X for(i=0;irow;i++) X for(j=0;j<3;j++) X d[i][j]=p->mat[i][j]; X /* depending on tg to change the values inside matrix X c is the current replaced column X Spiritly we move the new replaced column to the first column X */ X kk=p->se[p->bk-1][0]; X for(i=kk;irow;i++){ X /* when the replaced column is not inside the leading 0's */ X if(c>d[i][0]-1){ X /* when the replaced column is inside the lasting 0's */ X if(c>d[i][0]+d[i][1]-1){ X if(tg==0) X d[i][0] += 1; X if(tg==1){ X /* add 1 to the first row of the block and 0 in the rest */ X if(i==kk) X d[i][1] += 1; X else if(i!=kk) X d[i][0] += 1; X } X }else if(c<=d[i][0]+d[i][1]-1){ X /* when the replaced column is inside the 1's X when tg=1 the first row not change, X others, the number of 1's reduce 1 X 0's add 1 */ X if(tg==0 || (tg==1 && i!=kk)){ X d[i][0] += 1; X d[i][1] -= 1; X } X } X } X X if(i==kk||(i!=kk && d[i][0]==d[i-1][0]+d[i-1][1])) X d[i][2]=0; X else if(i!=kk && d[i][0]se[p->bk-1][2]*sizeof(int)); X if(delrw==NULL){ X fprintf(fp,"No room for delrw(organ)\n"); X exit(0); X } X for(i=0;ise[p->bk-1][2];i++) X delrw[i]=0; X for(i=1,k=kk+1;krow;k++,i++){ X if(delrw[i]==1) X continue; X if(d[k-1][0]==d[k][0]){ X if(d[k-1][1]<=d[k][1]) X delrw[i]=1; X }else if(d[k-1][0]row-1 (the last row to be deleted) no change needed; X For other k's, if d[k+1][0]>=d[k-1][0]+d[k-1][1] (nonoverlapped) X then d[k+1][2]=0; X else d[k+1][2] += d[k][2] X */ X for(i=0,l=k=kk;krow;i++,k++){ X if(delrw[i]==1){ X if(k==kk) X d[k+1][2]=0; X else if(krow-1){ X if(d[k+1][0] >= d[k-1][0]+d[k-1][1]) X d[k+1][2]=0; X else X d[k+1][2] += d[k][2]; X } X continue; X } X for(j=0;j<3;j++) X d[l][j]=d[k][j]; X l++; X } X /* rearrange the values in d[k][0] */ X if(kk!=0) X d[kk][0]=d[kk-1][0]+d[kk-1][1]; X else X d[kk][0]=0; X for(i=kk+1 ;imat=New_Array(l,3); X if(NEW->mat==NULL){ X fprintf(fp,"No room for NEW->mat(organ)\n"); X exit(0); X } X X for(i=0;imat[i][j]=d[i][j]; X } X } X X NEW->row=l; X NEW->num=0; X for(i=0;irow;i++) X NEW->num += NEW->mat[i][1]; X /* Here is wasting time to do it, Might need to fix it later X Still hard to say which is better */ X NEW->bk=0; /* initialize the block be 0 */ X NEW->N=NEW->P=NULL; X NEW->se=block(NEW); /* examine the blocks inside the matrix */ X /* copy down the list of p->C into NEW->C */ X NEW->C=NULL; X X NEW->C=cplist(NEW->C,p->C,r); X X p=createlink(p,NEW); X X return p; X X} X X/* cplist: copy the list of old to new */ XLIST *cplist(LIST *new,LIST *old,double sign) X{ X LIST *c,*temp,*tail; X X for(c=old;c!=NULL;c=c->N){ X temp=Calloc(); X if(temp==NULL){ X fprintf(fp,"no room for temp(cplist)\n"); X exit(0); X } X temp->j=c->j; X temp->k=c->k; X /* the odd term with coef 1 and the even term with -1 */ X if(sign==1.0) X temp->a=c->a; X else if(sign==0.0) X temp->a=c->a*(-1.0); X temp->N=temp->P=NULL; X X if(new==NULL){ /* the starting node in the list */ X new=temp; X }else{ /* add the node to the last of the list */ X tail=new; X while(tail->N!=NULL) X tail=tail->N; X temp->P=tail; X tail->N=temp; X } X X } X return new; X} X X/* X polynomials: transfer the matrix notation into modified polynomials Xp(interaction)=\sum_{0<=k<=l}{k \choose k1,k2,...,km}R(a,k) Xand R(a,k)={n \choose k}(1-sum ai)^(n-k)\prod ai^ki X*/ Xvoid polynomials(TREE *p) X{ X int i,j,k,ktotal; X double ans; X int *b; X LIST *c,*new; X double Fact(double); X LIST *cplist(LIST *,LIST *,double); X PP *addpoly(PP *,LIST *,int,int); X X /* the array to store the values of k's inside the multinomials. */ X b=(int *)malloc(p->row*sizeof(int)); X for(i=0;irow;i++) X b[i]=0; X /* from 0,...,0 to the values-1's inside p->e[i] X start from the first one(add 1) up to its max (value-1) and set to 0 X then the second one (add 1) the move the first one(add 1) to its max X and then set to 0; then move up 1 in the second one till it reached X its max and set to 0 on the second one and move to the third and X repeat the previous work again then move to the fourth and so on. */ X i=0; X while(1){ X for(k=ktotal=0;krow;k++) X ktotal += b[k]; X /* calculate {ktotal \choose b[0] ... b[p->row-1]} */ X for(ans=Fact(ktotal),k=0;krow;k++) X ans /= Fact(b[k]); X new=NULL; X new=cplist(new,p->C,1.0); X for(c=new;c!=NULL;c=c->N){ X c->a *= ans; X } X X R=addpoly(R,new,ktotal,p->row); X X if(b[i]+1e[i]) /* the first one increase step by step */ X b[i] += 1; X else{ X b[i]=0; X for(j=i+1;jrow;j++){ X b[j] += 1; X if(b[j]e[j]) /* not reach its max */ X break; X else /* reach the max set back to 0 and move to next one */ X b[j]=0; X } X if(j==p->row) /* increment reaches the last component */ X break; X i=0; /* back to the first one component */ X } X } X} X X XPP *palloc(void) X{ X return (PP *)malloc(sizeof(PP)); X} X XPP *addpoly(PP *o,LIST *d,int s,int t) X{ X X X LIST *combine(LIST *,LIST *); X PP *palloc(void); X LIST *c; X X if(o==NULL){ X o=palloc(); X if(o==NULL){ X fprintf(fp,"no room for o(addpoly)\n"); X exit(0); X } X o->C=NULL; X o->C=d; X o->s=s; X o->t=t; X o->L=o->R=NULL; X }else if(t>o->t){ X o->L=addpoly(o->L,d,s,t); X }else if(tt){ X o->R=addpoly(o->R,d,s,t); X }else if(t==o->t){ X if(s>o->s){ X o->L=addpoly(o->L,d,s,t); X }else if(ss){ X o->R=addpoly(o->R,d,s,t); X }else if(s==o->s){ X o->C=combine(o->C,d); X } X } X return o; X} X X Xvoid polyprint(PP *p) X{ X X LIST *c; X LIST *bcoef(LIST *); X X if(p!=NULL){ X polyprint(p->R); X X /* Now b(j,k)=b(j+1,k)+b(j+1,k-1) */ X X p->C=bcoef(p->C); X fprintf(fp,"+(\n"); X for(c=p->C;c!=NULL;c=c->N){ X if(c->a>0) X fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k); X else if(c->a<0) X fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k); X } X fprintf(fp,")*"); X fprintf(fp,"R(%d,%d)\n",p->s,p->t); X polyprint(p->L); X } X} X X/* X bcoef: use b(j,k)=b(j+1,k)+b(j+1,k-1) X find the largest value of j amongst k and when k=0 jump out X the loop. (the largest j is always in the first node among the X same k. Traverse the list and find the last node and execute X the fact and delete the node from the list and add two more X new nodes into the list. Continue this process until it reaches X the first node and move to next k. X*/ XLIST *bcoef(LIST *p) X{ X int jj,kk,num; X LIST *tail,*head,*end,*temp1,*temp2,*temp3; X LIST *combine(LIST *,LIST *); X X for(head=p,tail=p->N; tail!=NULL;){ X num=0; X jj=head->j; X kk=head->k; X if(kk==0) X break; X while(tail!=NULL){ X if(tail->k==kk){ X num++; X if(tail->N==NULL) X end=tail; X tail=tail->N; X }else{ X end=tail->P; X break; X } X } X if(num==0){ X head=end->N; X if(head!=NULL) X tail=head->N; X else X tail=NULL; X }else{ X tail=end; X while(head->j!=tail->j){ X temp1=Calloc(); X if(temp1==NULL){ X fprintf(fp,"No room for temp1(polyprint)\n"); X exit(0); X } X temp1->j=tail->j+1; X temp1->k=tail->k; X temp1->a=tail->a; X temp1->P=NULL; X temp2=Calloc(); X if(temp2==NULL){ X fprintf(fp,"No room for temp2(polyprint)\n"); X exit(0); X } X temp2->j=tail->j+1; X temp2->k=tail->k-1; X temp2->a=tail->a; X temp2->P=temp1; X temp1->N=temp2; X temp2->N=NULL; X combine(p,temp1); X temp3=tail->P; X temp3->N=tail->N; X tail->N->P=temp3; X free(tail); X tail=temp3; X } X if(head->N!=NULL){ X head=head->N; X tail=head->N; X }else{ X tail=NULL; X } X } X } X return p; X} X X/* X Delete the item in the given tree X Delete the item when there is no child. X*/ XTREE *Delete(TREE *p) X{ X int i; X TREE *temp; X TREE *Tfree(TREE *); X X if (p!= NULL){ X if(p->R==NULL && p->L==NULL){ X p=Tfree(p); X return p; X }else if(p->R==NULL){ X temp=p; X p=p->L; X p=Delete(p); X temp=Tfree(temp); X return temp; X }else if(p->L==NULL){ X temp=p; X p=p->R; X p=Delete(p); X temp=Tfree(temp); X return temp; X }else{ X temp=p; X p->L=Delete(p->L); X p->R=Delete(p->R); X temp=Tfree(temp); X return temp; X } X } else { X return p; X } X X} X X XTREE *Tfree(TREE *p) X{ X X LIST *listfree(LIST *); X X X free(p->e); X p->C=listfree(p->C); X free( p); X return NULL; X} X XLIST *listfree(LIST *p) X{ X LIST *c; X X c=p; X while(c!=NULL){ X if(c->N!=NULL){ X c=c->N; X free(c->P); X c->P=NULL; X }else{ X free(c); X c=NULL; X } X } X return NULL; X} X X Xvoid print_tree(TREE *p) X{ X int i; X LIST *c; X void polynomials(TREE *); X X if(p!= NULL){ X print_tree(p->R); X polynomials(p); X print_tree(p->L); X } X} X X XPP *Pfree(PP *p) X{ X X if(p!= NULL){ X p->R=Pfree(p->R); X p->L=Pfree(p->L); X if(p->R==NULL && p->L==NULL){ X p->C=listfree(p->C); X free(p); X } X } X return p; X} X X X X ________This_Is_The_END________ if test `wc -l < moment_four.c` -ne 1541; then echo 'shar: moment_four.c was damaged during transit (should have had 1541 lines)' fi echo 'x - moment_two.c' sed 's/^X//' << '________This_Is_The_END________' > moment_two.c X/* X moment_two.c 8-5-1997 X X Only print out the first two moments. X------------------------------------------ X X Execute the program as the followings. X gcc -o 'execute file' moment_two.c -lm X 'execute file' m 'output file' X*/ X X#include X#include X#include X Xtypedef struct cnode{ X int j, /* the values of b(j,k). */ X k; X double a; /* coefficients of the weight. */ X struct cnode *N,*P; X}LIST; X Xtypedef struct lnode{ X int **mat; /* DE matrix */ X int **se; /* starting and ending index matrix */ X int row, /* rows of the matrix. */ X num, /* 1's in side the matrix */ X bk; /* the number of block inside the matrix. */ X LIST *C; X struct lnode *N,*P; X}LINK; X Xtypedef struct tnode{ X int *e; X int row; /* dimension of the matrix. */ X LIST *C; X struct tnode *L,*R; X}TREE; X Xtypedef struct{ /* comparing the blocks inside the matrix */ X int **t; /* the submatirx of the block */ X int len; /* length of the block */ X}com; X Xtypedef struct pnode{ X int s,t; /* the values in R(s,t) */ X LIST *C; X struct pnode *L,*R; X}PP; X XFILE *fp; XLIST *W; XTREE *S; XPP *R; Xint clump; Xint **New_Array(int,int); XLINK *Lalloc(void); XTREE *Talloc(void); XLIST *Calloc(void); X Xmain(argc, argv) Xint argc; Xchar *argv[]; X{ X double ca,tca; X int **nmat; X int m,I,M,i,j,k,ii,jj,incre; X int colnum,cj,ck; X LINK *Head; X LIST *c; X double acoeff(int,int); X LIST *listfree(LIST *); X LINK *translate(LINK *,int **,double,int,int,int,int); X LINK *decom(LINK *); X LIST *bcoef(LIST *); X PP *Pfree(PP *); X TREE *Delete(TREE *); X void print_tree(TREE *); X void polyprint(PP *); X X X X if (argc != 3) X { X fprintf(stderr, "Usage: %s m outfilename\n", argv[0]); X return 1; X } X X X if((fp=fopen(argv[2],"a+"))==NULL){ /* data file not found */ X fprintf(stderr, "Cannot open data file %s\n", argv[2]); X exit(1); X } X X m = atoi(argv[1]); X X clump=m-1; X X X fprintf(fp,"m:=%d;\n",m); X X for(M=2;M>=1;M--){ X Head=NULL; X S=NULL; X R=NULL; X W=NULL; X if(M==1) X fprintf(fp,"m1"); X else if(M==2) X fprintf(fp,"m2"); X fprintf(fp,":='\n"); X for(I=1;I<=M;I++){ X tca=acoeff(M,I); X if(I==1){ X cj=m-1; X ck=1; X ca=tca; X colnum=clump; X nmat=New_Array(I,colnum); X for(i=0;iN){ X if(c->a>0) X fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k); X else if(c->a<0) X fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k); X } X fprintf(fp,")\n"); X X W=listfree(W); X polyprint(R); X X R=Pfree(R); X fprintf(fp,"';\n"); X } X fclose(fp); X} X X X X/* New_Array: make an allocation for a matrix. */ Xint **New_Array(int row, int col) X{ X int **A,i,j; X X A=(int **)malloc(row*sizeof(int*)); X if(A==NULL) X return NULL; X for(i=0;i a) X return 0.0; X else if (a<=0 || b<=0) X return 0.0; X else X return b*acoeff(a-1,b)+b*acoeff(a-1,b-1); X} X X XLINK *translate(LINK *rt,int **d,double ca,int cj,int ck,int row,int col) X{ X int asize,i,j,k,ii,jj; X int **a,**b,*delcl,*indices,**temp; X int star,location,jump,count,ncol; X LINK *NEW; X LIST *cnew; X double Comb(double,double); X int **block(LINK *); X LINK *createlink(LINK *,LINK *); X LIST *combine(LIST *,LIST *); X void permute(int,int,int *,int,int **,int *); X X X cnew=Calloc(); X if(cnew==NULL){ X fprintf(fp,"no room for temp(cplist)\n"); X exit(0); X } X cnew->j=cj; X cnew->k=ck; X cnew->a=ca; X cnew->N=cnew->P=NULL; X if(W==NULL) X W=cnew; X else X W=combine(W,cnew); X X delcl=(int *)malloc(col*sizeof(int)); X if(delcl==NULL){ X fprintf(fp,"No room for delcl(translate)\n"); X exit(0); X } X indices=(int *)malloc(row*sizeof(int)); X if(indices==NULL){ X fprintf(fp,"No room for indices(translate)\n"); X exit(0); X } X for(k=1;k<=row;k++){ X ca *= -1.0; X location=0; X star=k; X asize=Comb(row,star); X /* This is achievable job since the maximum row is 4.*/ X a=New_Array(asize,row); X if(a==NULL){ X fprintf(fp,"No room for a(translate)\n"); X exit(0); X } X for(i=0;imat=b; /* cp the DE matrix to node */ X NEW->row=k; X NEW->num=0; /* count 1's inside the matrix */ X for(i=0;irow;i++) X NEW->num += NEW->mat[i][1]; X NEW->bk=0; /* initialize the block be 0 */ X NEW->N=NEW->P=NULL; X NEW->se=block(NEW); /* examine the blocks inside the matrix */ X NEW->C=Calloc(); X if(NEW->C==NULL){ X fprintf(fp,"No room for NEW->C(translate)\n"); X exit(0); X } X NEW->C->j=cj; X NEW->C->k=ck; X NEW->C->a=ca; X NEW->C->N=NEW->C->P=NULL; X rt=createlink(rt,NEW); X } X for(i=0;i= m) X { X for (i = 0; i < n; i++) X { X int j; X int value = 0; /* first assumes the value as 1 */ X X /* if i is not in the indIndices, then it's 0 */ X for (j = 0; !value && j < whichIndex; j++) X { X value = i == indIndices[j]; X } X a[*location][i]=value; X if(i==n-1) X (*location)++; X } X return; X } X for (i = whichIndex == 0 ? 0 : indIndices[whichIndex - 1] + 1; X i <= n - (m - whichIndex); X i++) X { X indIndices[whichIndex] = i; X permute(n, m, indIndices, whichIndex + 1,a,location); X } X return; X} X Xdouble Fact(double n) X{ X if(n==0.0) X return 1.0; X else X return n*Fact(n-1.0); X} X Xdouble Comb(double a, double b) X{ X int i,ans=1; X double Fact(double); X X if(a < b) X return 0.0; X else{ X for(i=b;i>=1.0;i--,a--) X ans=ans*a; X ans=ans/Fact(b); X return ans; X } X} X X/* createlink: add the node to the linklist */ XLINK *createlink(LINK *rt,LINK *NEW) X{ X int i,j,match; X LINK *cur,*prev; X LINK *addchecking(LINK *,LINK *,LINK *,LINK *); X LIST *combine(LIST *,LIST *); X TREE *addtree(TREE *,LINK *); X X X/* check if the new matrix is simple matrix or not */ X if(NEW->bk==0){ X /*simple matrix call the shorthand fun. */ X S=addtree(S,NEW); X return rt; X }else{ X if(rt==NULL){ X return NEW; X }else{ X for(prev=NULL,cur=rt;cur!=NULL;prev=cur,cur=cur->N){ X /* larger size of matrix first be decomposed */ X if(NEW->rowrow) X continue; X else if(NEW->row>cur->row) X return addchecking(rt,prev,cur,NEW); X else if(NEW->row==cur->row){ X /* more 1's inside the matrix first be decomposed */ X if(NEW->numnum) X continue; X else if(NEW->num>cur->num) X return addchecking(rt,prev,cur,NEW); X else if(NEW->num==cur->num){ X /* larger size of block is decomposed first */ X if(NEW->se[NEW->bk-1][2]se[cur->bk-1][2]) X continue; X else if(NEW->se[NEW->bk-1][2]>cur->se[cur->bk-1][2]) X return addchecking(rt,prev,cur,NEW); X else if(NEW->se[NEW->bk-1][2]==cur->se[cur->bk-1][2]){ X /*check if the identical term */ X for(i=0,match=1;irow && match==1;i++){ X for(j=0;j<3 && match==1;j++){ X if(NEW->mat[i][j]!=cur->mat[i][j]) X match=0; X } X } X if(match==0) /* Not identical, move to next node to check. */ X continue; X else{ /* Identical, add coeffs to LIST */ X cur->C=combine(cur->C,NEW->C); X for(i=0;ibk;i++) /* free the node */ X free(NEW->se[i]); X free(NEW->se); X for(i=0;irow;i++) X free(NEW->mat[i]); X free(NEW->mat); X free(NEW); X return rt; X } X } X } X } X } X /* end of for-loop */ X /*connect the new node to the end of list */ X prev->N=NEW; X NEW->P=prev; X return rt; X } /* end of else-when rt!=NULL */ X } X} X X/* addchecking: add the node into the appropriate place */ XLINK *addchecking(LINK *rt,LINK *prev,LINK *cur,LINK *NEW) X{ X if(prev==NULL){ X NEW->N=cur; X cur->P=NEW; X return NEW; X }else{ X NEW->N=cur; X NEW->P=prev; X prev->N=NEW; X cur->P=NEW; X return rt; X } X} X X/* combine: combine the two list */ XLIST *combine(LIST *old,LIST *new) X{ X int jp; X LIST *p,*c,*temp; X LIST *listcheck(LIST *,LIST *,LIST *,LIST *); X X /* old has at least one node in the LIST */ X while(new!=NULL){ X temp=new; X new=new->N; X if(new!=NULL) X new->P=NULL; X temp->N=NULL; X for(p=NULL,c=old,jp=0;c!=NULL;p=c,c=c->N){ X if(temp->kk) X continue; X else if(temp->k>c->k){ X old=listcheck(old,p,c,temp); X jp=1; X break; X }else if(temp->k==c->k){ X if(temp->jj) X continue; X else if(temp->j>c->j){ X old=listcheck(old,p,c,temp); X jp=1; X break; X }else if(temp->j==c->j){ X c->a += temp->a; X jp=1; X free(temp); X break; X } X X } X } X if(jp==0){ X p->N=temp; X temp->P=p; X } X } X return old; X} X X/* listcheck: add the node into the appropriate place */ XLIST *listcheck(LIST *rt,LIST *prev,LIST *cur,LIST *NEW) X{ X if(prev==NULL){ X NEW->N=cur; X cur->P=NEW; X return NEW; X }else{ X NEW->N=cur; X NEW->P=prev; X prev->N=NEW; X cur->P=NEW; X return rt; X } X} X X/* block: check the blocks inside the matrix and return the number of block and matrix */ Xint **block(LINK *b) X{ X X int i,j,k,l,kk,ll,end; X int *bk=&b->bk; /* initialize to be 0 */ X int **temp; X com *c; X int Num_Cmp1( const void * , const void * ); X int Num_Cmp2( const void * , const void * ); X X if(b->row==1) X return NULL; X /* record the matrix where is not overlapped with other rows */ X temp=New_Array(b->row,3); X if(temp==NULL){ X fprintf(fp,"No room for temp(block)\n"); X exit(0); X } X /* X Good for b->row > 2 X When b->mat[i][2]=0 and b->mat[i+1][2]!=0 it starts the block X b->mat[i][2]!=0 and b->mat[i+1][2]=0 it ends the block X b->se[i][0] records the starting row X b->se[i][1] records the ending row X b->se[i][2] records the number of rows inside the block X */ X for(k=0,i=end=1;irow;i++){ X if(b->mat[i][2]==0 && b->mat[i-1][2]==0){ X for(ll=0;ll<3;ll++) /* cp the nonoverlapped row */ X temp[k][ll]=b->mat[i-1][ll]; X k++; X }else if(b->mat[i][2]!=0 && b->mat[i-1][2]!=0){ /* still inside the block */ X continue; X }else if(b->mat[i][2]==0 && b->mat[i-1][2]!=0){ /* it ends the block */ X b->se[*bk][1]=i-1; /* the ending row */ X b->se[*bk][2]=i-b->se[*bk][0]; /* the number of rows inside this block */ X (*bk)++; /* increment to next block */ X end=1; /* index for ending the block */ X }else if(b->mat[i][2]!=0 && b->mat[i-1][2]==0){ /* it starts the block */ X if(*bk==0){ X b->se=New_Array((int)b->row/2,3); /* The maximum blocks will be row/2 */ X if(b->se==NULL){ X fprintf(fp,"No room for b->se(block)\n"); X exit(0); X } X } X b->se[*bk][0]=i-1; /* the staring row of the first block */ X end=0; /* index for starting the block */ X } X if(i==b->row-1 && b->mat[i][2]==0){ /* examine the last row */ X for(ll=0;ll<3;ll++) X temp[k][ll]=b->mat[i][ll]; X k++; X } X } X if(end==0){ /* completes the record of the block */ X b->se[*bk][1]=b->row-1; X b->se[*bk][2]=b->row-b->se[*bk][0]; X (*bk)++; X } X /* if it is simple, sort(increasing) the number of 1's among rows */ X if(*bk==0 && k==b->row){ X qsort(b->mat,b->row,sizeof(int *),Num_Cmp1); X for(i=0;irow;i++){ X if(i==0) X b->mat[i][0]=0; X else X b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1]; X } X b->se=NULL; X }else{ /* There is at least one block inside the matrix */ X for(i=*bk;i<(int)b->row/2;i++) /* Free the un-used part */ X free(b->se[i]); X /* if k!=0 then put those unoverlapped rows to the front and sort them X Also compare the size of the blocks and sort them at the back */ X if(k!=0){ X if(b->se[0][0]!=k){ X for(i=0,kk=k;i<*bk;i++){ X l=kk; X for(j=b->se[i][0];j<=b->se[i][1];j++){ X for(ll=0;ll<3;ll++) X temp[kk][ll]=b->mat[j][ll]; X kk++; X } X b->se[i][0]=l; X b->se[i][1]=kk-1; X } X for(i=0;irow;i++){ X for(ll=0;ll<3;ll++) X b->mat[i][ll]=temp[i][ll]; X } X qsort(b->mat,b->se[0][0],sizeof(int *),Num_Cmp1); X for(i=0;ise[0][0];i++){ X if(i==0) X b->mat[i][0]=0; X else X b->mat[i][0]=b->mat[i-1][0]+b->mat[i-1][1]; X } X } X } X if(*bk>=2){ X c=(com *)malloc(*bk*sizeof(com)); X if(c==NULL){ X fprintf(fp,"No c space(block)\n"); X exit(0); X } X for(i=0;i<*bk;i++){ X c[i].t=New_Array(b->se[i][2],3); X if(c[i].t==NULL){ X fprintf(fp,"No c.t space(block)\n"); X exit(0); X } X for(j=0,kk=b->se[i][0];kk<=b->se[i][1];j++,kk++){ X for(ll=0;ll<3;ll++) X c[i].t[j][ll]=b->mat[kk][ll]; X } X c[i].len=b->se[i][2]; X } X qsort(c, *bk, sizeof(com), Num_Cmp2); X for(i=0,j=b->se[0][0];i<*bk;i++){ X for(kk=0;kkmat[j][ll]=c[i].t[kk][ll]; X } X /* rearrange the values inside b->se */ X b->se[i][2]=c[i].len; X if(i!=0) X b->se[i][0]=b->se[i-1][1]+1; X else X b->se[0][0]=k; X b->se[i][1]=b->se[i][0]+b->se[i][2]-1; X for(kk=0;kkmat[i][0] */ X X for(i=0,j=b->se[0][0];i<*bk;i++){ X for(kk=0;kkse[i][2];kk++,j++){ X if(kk==0){ X if(i==0 && k==0) X b->mat[0][0]=0; X else X b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1]; X }else X b->mat[j][0]=b->mat[j-1][0]+b->mat[j-1][1]+b->mat[j][2]-b->mat[j][1]; X } X } X } X for(i=0;irow;i++) X free(temp[i]); X free(temp); X return b->se; X} X Xint Num_Cmp1( const void * S1, const void * S2) X{ X int ** T1 = (int **) S1, X ** T2 = (int **) S2; X X if(*(*T1+1) < *(*T2+1)) X return -1; X else if (*(*T1+1) > *(*T2+1)) X return 1; X else X return 0; X} X Xint Num_Cmp2( const void * S1, const void * S2) X{ X com * T1 = (com *) S1, X * T2 = (com *) S2; X X if(T1->len < T2->len) X return -1; X else if (T1->len > T2->len) X return 1; X else X return 0; X} X X X XTREE *addtree(TREE *p,LINK *I) X{ X int i,match; X LIST *combine(LIST *,LIST *); X X if(p==NULL){ X p=Talloc(); X if(p==NULL){ X fprintf(fp,"No room for p(addtree)\n"); X exit(0); X } X p->e=(int *)malloc(I->row*sizeof(int)); X if(p->e==NULL){ X fprintf(fp,"No room for p->e(addtree)\n"); X exit(0); X } X for(i=0;irow;i++) X p->e[i]=I->mat[i][1]; X p->row=I->row; X p->C=I->C; X p->L=p->R=NULL; X for(i=0;irow;i++) X free(I->mat[i]); X free(I->mat); X /* I->se is NULL, no need to free it */ X free(I); X }else if(I->row>p->row){ /* if row is greater, goes to Left Tree */ X p->L=addtree(p->L,I); X }else if(I->rowrow){ /* if row is lesser, goes to Right Tree */ X p->R=addtree(p->R,I); X }else if(I->row==p->row){ X for(i=match=0;irow;i++){ X if(I->mat[i][1]==p->e[i]){ X continue; X }else if(I->mat[i][1]>p->e[i]){ X match=1; X break; X }else if(I->mat[i][1]e[i]){ X match=2; X break; X } X } X if(match==1){ /* Not identical, insert to Left Tree */ X p->L=addtree(p->L,I); X }else if(match==2){ X p->R=addtree(p->R,I); X }else if(match==0){ X /* combine two link list of coefficients */ X p->C=combine(p->C,I->C); X for(i=0;irow;i++) X free(I->mat[i]); X free(I->mat); X free(I); X } X } X return p; X} X X/* decom: decompose the matrices inside the link list */ XLINK *decom(LINK *p) X{ X int i,j,col,len,tg; X int *loc; X LIST *cr; X int *coef(LINK *,int *,int *,int *); X LINK *organ(LINK *,int,int,double); X LINK *Lfree(LINK *); X X while(p!=NULL){ X /* We only look at the last big block and X put corresponding columns to coeffs inside loc */ X col=p->mat[p->row-1][0]+p->mat[p->row-1][1]; X loc=(int *)malloc(col*sizeof(int)); X if(loc==NULL){ X fprintf(fp,"No room for loc(decom)\n"); X exit(0); X } X /* initialize the values of array loc */ X for(i=0;iN!=NULL){ X p=p->N; X p->P=Lfree(p->P); X }else{ X p=Lfree(p); X } X } X return p; X} X X/* Lfree: To free the node inside the link list */ XLINK *Lfree(LINK *p) X{ X int i; X LIST *cr; X X for(i=0;ibk;i++) X free(p->se[i]); X free(p->se); X for(i=0;irow;i++) X free(p->mat[i]); X free(p->mat); X cr=p->C; X while(cr!=NULL){ X if(cr->N!=NULL){ X cr=cr->N; X free(cr->P); X cr->P=NULL; X }else{ X free(cr); X cr=NULL; X } X } X free(p); X return NULL; X} X X X/* coef: find the location of coefficients and return the taget */ Xint *coef(LINK *p,int *c,int *tar,int *l) X{ X int mr,k,sum,i,j; X int *fill(LINK *,int,int *,int *); X X mr=k=p->row-1; /* start to mark the set at the last row of the block */ X for(sum=0;k>=p->se[p->bk-1][0];k--){ X sum=sum+p->mat[k][2]; X if(sum>=p->mat[mr][1]){ X c=fill(p,mr,c,l); X /* if equal continue to mark set X not equal means stop and tg=0 and return */ X if(sum==p->mat[mr][1]){ X mr=k-1; /* the next marking row and initialize sum=0 */ X sum=0; X }else if (sum>p->mat[mr][1]){ /* break here, no need to continue mark set */ X *tar=0; /* target is 0 vector */ X return c; X } X } X /* if m* is the first row of the block */ X if(mr==p->se[p->bk-1][0]){ X *tar=1; /* target is (1,0,...0) */ X return c; X } X } X /* if m* is not the first row of the block */ X c=fill(p,mr,c,l); X *tar=0; /* target is 0 vector */ X return c; X} X X/* fill: find the ending and starting columns of current marking row */ Xint *fill(LINK *p,int mr,int *c,int *l) X{ X if(mr==p->row-1){ X c[(*l)++]=p->mat[mr][0]+p->mat[mr][1]-1; /* the last ending column */ X } X c[(*l)++]=p->mat[mr][0]; /* the starting column of current one */ X c[(*l)++]=p->mat[mr][0]-1; /* the ending column of next one */ X return c; X} X X/* organ: organize the new produced matrix X (also delet the redundant rows) and to the link list */ XLINK *organ(LINK *p,int c,int tg,double r) X{ X int i,j,k,kk,l,jp; X int **d,*delrw; X LINK *NEW; X LIST *cur; X LIST *cplist(LIST *,LIST *,double); X X /* d is the replaced matrix accroding to tg */ X d=New_Array(p->row,3); X if(d==NULL){ X fprintf(fp,"No room for d(organ)\n"); X exit(0); X } X for(i=0;irow;i++) X for(j=0;j<3;j++) X d[i][j]=p->mat[i][j]; X /* depending on tg to change the values inside matrix X c is the current replaced column X Spiritly we move the new replaced column to the first column X */ X kk=p->se[p->bk-1][0]; X for(i=kk;irow;i++){ X /* when the replaced column is not inside the leading 0's */ X if(c>d[i][0]-1){ X /* when the replaced column is inside the lasting 0's */ X if(c>d[i][0]+d[i][1]-1){ X if(tg==0) X d[i][0] += 1; X if(tg==1){ X /* add 1 to the first row of the block and 0 in the rest */ X if(i==kk) X d[i][1] += 1; X else if(i!=kk) X d[i][0] += 1; X } X }else if(c<=d[i][0]+d[i][1]-1){ X /* when the replaced column is inside the 1's X when tg=1 the first row not change, X others, the number of 1's reduce 1 X 0's add 1 */ X if(tg==0 || (tg==1 && i!=kk)){ X d[i][0] += 1; X d[i][1] -= 1; X } X } X } X X if(i==kk||(i!=kk && d[i][0]==d[i-1][0]+d[i-1][1])) X d[i][2]=0; X else if(i!=kk && d[i][0]se[p->bk-1][2]*sizeof(int)); X if(delrw==NULL){ X fprintf(fp,"No room for delrw(organ)\n"); X exit(0); X } X for(i=0;ise[p->bk-1][2];i++) X delrw[i]=0; X for(i=1,k=kk+1;krow;k++,i++){ X if(delrw[i]==1) X continue; X if(d[k-1][0]==d[k][0]){ X if(d[k-1][1]<=d[k][1]) X delrw[i]=1; X }else if(d[k-1][0]row-1 (the last row to be deleted) no change needed; X For other k's, if d[k+1][0]>=d[k-1][0]+d[k-1][1] (nonoverlapped) X then d[k+1][2]=0; X else d[k+1][2] += d[k][2] X */ X for(i=0,l=k=kk;krow;i++,k++){ X if(delrw[i]==1){ X if(k==kk) X d[k+1][2]=0; X else if(krow-1){ X if(d[k+1][0] >= d[k-1][0]+d[k-1][1]) X d[k+1][2]=0; X else X d[k+1][2] += d[k][2]; X } X continue; X } X for(j=0;j<3;j++) X d[l][j]=d[k][j]; X l++; X } X /* rearrange the values in d[k][0] */ X if(kk!=0) X d[kk][0]=d[kk-1][0]+d[kk-1][1]; X else X d[kk][0]=0; X for(i=kk+1 ;imat=New_Array(l,3); X if(NEW->mat==NULL){ X fprintf(fp,"No room for NEW->mat(organ)\n"); X exit(0); X } X X for(i=0;imat[i][j]=d[i][j]; X } X } X X NEW->row=l; X NEW->num=0; X for(i=0;irow;i++) X NEW->num += NEW->mat[i][1]; X /* Here is wasting time to do it, Might need to fix it later X Still hard to say which is better */ X NEW->bk=0; /* initialize the block be 0 */ X NEW->N=NEW->P=NULL; X NEW->se=block(NEW); /* examine the blocks inside the matrix */ X /* copy down the list of p->C into NEW->C */ X NEW->C=NULL; X X NEW->C=cplist(NEW->C,p->C,r); X X p=createlink(p,NEW); X X return p; X X} X X/* cplist: copy the list of old to new */ XLIST *cplist(LIST *new,LIST *old,double sign) X{ X LIST *c,*temp,*tail; X X for(c=old;c!=NULL;c=c->N){ X temp=Calloc(); X if(temp==NULL){ X fprintf(fp,"no room for temp(cplist)\n"); X exit(0); X } X temp->j=c->j; X temp->k=c->k; X /* the odd term with coef 1 and the even term with -1 */ X if(sign==1.0) X temp->a=c->a; X else if(sign==0.0) X temp->a=c->a*(-1.0); X temp->N=temp->P=NULL; X X if(new==NULL){ /* the starting node in the list */ X new=temp; X }else{ /* add the node to the last of the list */ X tail=new; X while(tail->N!=NULL) X tail=tail->N; X temp->P=tail; X tail->N=temp; X } X X } X return new; X} X X/* X polynomials: transfer the matrix notation into modified polynomials Xp(interaction)=\sum_{0<=k<=l}{k \choose k1,k2,...,km}R(a,k) Xand R(a,k)={n \choose k}(1-sum ai)^(n-k)\prod ai^ki X*/ Xvoid polynomials(TREE *p) X{ X int i,j,k,ktotal; X double ans; X int *b; X LIST *c,*new; X double Fact(double); X LIST *cplist(LIST *,LIST *,double); X PP *addpoly(PP *,LIST *,double,int,int); X X /* the array to store the values of k's inside the multinomials. */ X b=(int *)malloc(p->row*sizeof(int)); X for(i=0;irow;i++) X b[i]=0; X /* from 0,...,0 to the values-1's inside p->e[i] X start from the first one(add 1) up to its max (value-1) and set to 0 X then the second one (add 1) the move the first one(add 1) to its max X and then set to 0; then move up 1 in the second one till it reached X its max and set to 0 on the second one and move to the third and X repeat the previous work again then move to the fourth and so on. */ X i=0; X while(1){ X for(k=ktotal=0;krow;k++) X ktotal += b[k]; X /* calculate {ktotal \choose b[0] ... b[p->row-1]} */ X for(ans=Fact(ktotal),k=0;krow;k++) X ans /= Fact(b[k]); X new=NULL; X new=cplist(new,p->C,1.0); X for(c=new;c!=NULL;c=c->N){ X c->a *= ans; X } X X R=addpoly(R,new,ans,ktotal,p->row); X X if(b[i]+1e[i]) /* the first one increase step by step */ X b[i] += 1; X else{ X b[i]=0; X for(j=i+1;jrow;j++){ X b[j] += 1; X if(b[j]e[j]) /* not reach its max */ X break; X else /* reach the max set back to 0 and move to next one */ X b[j]=0; X } X if(j==p->row) /* increment reaches the last component */ X break; X i=0; /* back to the first one component */ X } X } X} X X XPP *palloc(void) X{ X return (PP *)malloc(sizeof(PP)); X} X XPP *addpoly(PP *o,LIST *d,double kcom,int s,int t) X{ X X X LIST *combine(LIST *,LIST *); X PP *palloc(void); X LIST *c; X X if(o==NULL){ X o=palloc(); X if(o==NULL){ X fprintf(fp,"no room for o(addpoly)\n"); X exit(0); X } X o->C=NULL; X o->C=d; X o->s=s; X o->t=t; X o->L=o->R=NULL; X }else if(t>o->t){ X o->L=addpoly(o->L,d,kcom,s,t); X }else if(tt){ X o->R=addpoly(o->R,d,kcom,s,t); X }else if(t==o->t){ X if(s>o->s){ X o->L=addpoly(o->L,d,kcom,s,t); X }else if(ss){ X o->R=addpoly(o->R,d,kcom,s,t); X }else if(s==o->s){ X o->C=combine(o->C,d); X } X } X return o; X} X X Xvoid polyprint(PP *p) X{ X X LIST *c; X LIST *bcoef(LIST *); X X if(p!=NULL){ X polyprint(p->R); X X /* Now b(j,k)=b(j+1,k)+b(j+1,k-1) */ X X p->C=bcoef(p->C); X fprintf(fp,"+(\n"); X for(c=p->C;c!=NULL;c=c->N){ X if(c->a>0) X fprintf(fp,"+%.0f*b(%d,%d)\n",c->a,c->j,c->k); X else if(c->a<0) X fprintf(fp,"%.0f*b(%d,%d)\n",c->a,c->j,c->k); X } X fprintf(fp,")*"); X fprintf(fp,"R(%d,%d)\n",p->s,p->t); X polyprint(p->L); X } X} X X/* X bcoef: use b(j,k)=b(j+1,k)+b(j+1,k-1) X find the largest value of j amongst k and when k=0 jump out X the loop. (the largest j is always in the first node among the X same k. Traverse the list and find the last node and execute X the fact and delete the node from the list and add two more X new nodes into the list. Continue this process until it reaches X the first node and move to next k. X*/ XLIST *bcoef(LIST *p) X{ X int jj,kk,num; X LIST *tail,*head,*end,*temp1,*temp2,*temp3; X LIST *combine(LIST *,LIST *); X X for(head=p,tail=p->N; tail!=NULL;){ X num=0; X jj=head->j; X kk=head->k; X if(kk==0) X break; X while(tail!=NULL){ X if(tail->k==kk){ X num++; X if(tail->N==NULL) X end=tail; X tail=tail->N; X }else{ X end=tail->P; X break; X } X } X if(num==0){ X head=end->N; X if(head!=NULL) X tail=head->N; X else X tail=NULL; X }else{ X tail=end; X while(head->j!=tail->j){ X temp1=Calloc(); X if(temp1==NULL){ X fprintf(fp,"No room for temp1(polyprint)\n"); X exit(0); X } X temp1->j=tail->j+1; X temp1->k=tail->k; X temp1->a=tail->a; X temp1->P=NULL; X temp2=Calloc(); X if(temp2==NULL){ X fprintf(fp,"No room for temp2(polyprint)\n"); X exit(0); X } X temp2->j=tail->j+1; X temp2->k=tail->k-1; X temp2->a=tail->a; X temp2->P=temp1; X temp1->N=temp2; X temp2->N=NULL; X combine(p,temp1); X temp3=tail->P; X temp3->N=tail->N; X tail->N->P=temp3; X free(tail); X tail=temp3; X } X if(head->N!=NULL){ X head=head->N; X tail=head->N; X }else{ X tail=NULL; X } X } X } X return p; X} X X X/* X Delete the item in the given tree X Delete the item when there is no child. X*/ XTREE *Delete(TREE *p) X{ X int i; X TREE *temp; X TREE *Tfree(TREE *); X X if (p!= NULL){ X if(p->R==NULL && p->L==NULL){ X p=Tfree(p); X return p; X }else if(p->R==NULL){ X temp=p; X p=p->L; X p=Delete(p); X temp=Tfree(temp); X return temp; X }else if(p->L==NULL){ X temp=p; X p=p->R; X p=Delete(p); X temp=Tfree(temp); X return temp; X }else{ X temp=p; X p->L=Delete(p->L); X p->R=Delete(p->R); X temp=Tfree(temp); X return temp; X } X } else { X return p; X } X X} X X XTREE *Tfree(TREE *p) X{ X X LIST *listfree(LIST *); X X X free(p->e); X p->C=listfree(p->C); X free( p); X return NULL; X} X XLIST *listfree(LIST *p) X{ X LIST *c; X X c=p; X while(c!=NULL){ X if(c->N!=NULL){ X c=c->N; X free(c->P); X c->P=NULL; X }else{ X free(c); X c=NULL; X } X } X return NULL; X} X X Xvoid print_tree(TREE *p) X{ X int i; X LIST *c; X void polynomials(TREE *); X X if(p!= NULL){ X print_tree(p->R); X polynomials(p); X print_tree(p->L); X } X} X X XPP *Pfree(PP *p) X{ X X if(p!= NULL){ X p->R=Pfree(p->R); X p->L=Pfree(p->L); X if(p->R==NULL && p->L==NULL){ X p->C=listfree(p->C); X free(p); X } X } X return p; X} X X X X ________This_Is_The_END________ if test `wc -l < moment_two.c` -ne 1452; then echo 'shar: moment_two.c was damaged during transit (should have had 1452 lines)' fi echo 'x - maplecode' sed 's/^X//' << '________This_Is_The_END________' > maplecode X X################################################### X# Date: July 31, 1996 X# By: Fred W. Huffer, Florida State Univ. X# Chien-Tai Lin, Tamkang Univ. X# X# Content: Maple code for computing approximations X# and bounds given in the JASA article by X# Huffer and Lin (1997), pp 1466-1475. X####################################################### X X XDigits:=20; X Xreadmom:=proc(mval) X###################################################################### X# This procedure reads in the expressions for the first four moments X# which are named m1, m2, m3, m4. X# Example: To get the expressions for m=7, execute readmom(7). X###################################################################### Xglobal m,mmmmm,m1,m2,m3,m4; Xm:=mval; Xmmmmm:=mval; Xread `fourmom.`.mval; X################################################################# X# Read the output files from your own directory inside the codes. X################################################################# Xend; X X############################################## X Xsetval:=proc(dval,nval) X########################################################################## X# This computes the first four moments and first four cumulants X# for d=dval and n=nval. X# These are named mom1, mom2, mom3, mom4, c1, c2, c3, c4. X# You must execute the procedure readmom before using setval X# in order to specify the value of m and read in the expressions X# for the moments. However, as long as you stick with the same X# value of m, you only need to execute readmom once. You can then run X# setval as many times as you like. X# You may wish to set Digits before using setval. X########################################################################## Xglobal m,mmmmm,d,ddddd,n,nnnnn,w,wwwww, X c1,c2,c3,c4, mom1,mom2,mom3,mom4; Xd:=dval; Xn:=nval; Xddddd:=dval; Xnnnnn:=nval; Xwwwww:=nnnnn-mmmmm+1; Xw:=wwwww; Xmom1:=evalf(m1); Xmom2:=evalf(m2); Xmom3:=evalf(m3); Xmom4:=evalf(m4); Xc1:=mom1; Xc2:=mom2-mom1^2; Xc3:=mom3-2*mom1*c2-mom2*mom1; Xc4:=mom4-3*mom1*c3-3*mom2*c2-mom3*mom1; XNULL; Xend; X X############################ X Xb:=proc(j,k) X Heaviside(n-j-k)*binomial(n-j,k); Xend; X XR:=proc(s,t) X if 1-t*d<=0 then X 0; X else X binomial(n,s)*d^s*(1-t*d)^(n-s); X fi; Xend; X X############################################################## X# Evaluate the first two moments from G and F for larger m. X# Faster Routine to get the first two moments X############################################################## X Xacc2:=proc(a,dim) Xlocal i,j,tot; Xfor j from 1 to dim do Xtot:=0; Xfor i from 0 to j-1 do Xtot:=tot+a[i,j]; Xa[i,j]:=a[i,j-1]+tot; Xod; Xa[j,j]:=a[j-1,j-1]+2*tot+a[j,j] Xod; Xeval(a); Xend; X X Xacc1:=proc(a,dim) Xlocal j; Xfor j from 1 to dim do Xa[j]:=a[j-1]+a[j]; Xod; Xeval(a); Xend; X X#################### X# The first moment X#################### X Xonemom:=proc(mm,dd,nn) Xlocal y,i,g,ans; Xy:=1-dd; Xif y<=0 then RETURN(`Failed.`); fi; Xg:=array(0..(mm-2),[seq(evalf(binomial(nn,i)*dd^i*y^(nn-i)),i=0..(mm-2))]); Xacc1(g,mm-2); Xans:=(nn-mm+1)*(1-g[mm-2]); Xans; Xend; X X#################### X# The second moment X#################### X Xtwomom:=proc(mm,dd,nn) Xlocal f0,i,j,y,z,g,f,ans; Xy:=1-dd; Xz:=1-2*dd; Xif y <= 0 or z <=0 then RETURN(`Failed.`); fi; Xg:=array(0..(mm-2),[seq(evalf(binomial(nn,i)*dd^i*y^(nn-i)),i=0..(mm-2))]); Xf:=array(0..(mm-2),0..(mm-2)); Xfor j from 0 to (mm-2) do for i from 0 to j do Xf[i,j] := evalf(binomial(nn,i)*binomial(nn-i,j)*d^(i+j)*z^(nn-i-j)); Xod; Xod; Xans:=0; Xacc1(g,mm-2); Xans:=(nn-mm+1)*(1-g[mm-2])+(nn-mm+1)*(nn-mm)*(1-2*g[mm-2]); Xacc1(g,mm-3); Xacc1(g,mm-3); Xans:=ans+4*(nn-mm)*g[mm-3]; Xacc1(g,mm-4); Xacc1(g,mm-4); Xans:=ans-12*g[mm-4]; Xacc2(f,mm-2); Xans:=ans+(nn-2*mm+3)*(nn-2*mm+2)*f[mm-2,mm-2]; Xacc2(f,mm-3); Xans:=ans-2*(nn-2*mm+3)*f[mm-3,mm-3]; Xacc2(f,mm-4); Xans:=ans+2*f[mm-4,mm-4]; Xans; Xend; X X########################################################### X# This procedure computes the first two moments and thus X# the first two cumulants by faster routine in order to X# get the result for larger m=mval,d=dval,n=nval. X# Run mc2() and cpg2(1) procedures after this procedure can X# obtain the values of Table 2. X############################################################ Xfsmval:=proc(mval,dval,nval) X global n,m,d,w,mom1,mom2,c1,c2; X n:=nval; X d:=dval; X m:=mval; X w:=n-m+1; X if m<=3 then ERROR(`Try bigger m`); fi; X mom1:=onemom(m,d,n); X if type(mom1,string) then ERROR(`Try another d.`); fi; X if n>=2*(m-1) then X mom2:=twomom(m,d,n); X if type(mom2,string) then ERROR(`Try another d.`); fi; X else X ERROR(`Failed. n<2(m-1). Try another n or m.`); X fi; X c1:=mom1; X c2:=mom2-mom1^2; X NULL; Xend; X Xmax_val_i:=30; X X###################################################### X# BOUNDS: max_val_i can be changed if you like. X###################################################### Xbounds:=proc() Xlocal i,j,p,ps,lb,ub,y; Xi:='i'; Xj:='j'; Xp:=collect(simplify(expand(1-(y-i)*(y-i-1)*(y-j)*(y-j-1)/(i*(i+1)*j*(j+1)))),y); Xps:=unapply(subs(y^4=mom4,y^3=mom3,y^2=mom2,y=mom1,p),i,j); Xlb:=max(seq(seq(ps(i,j),j=(i+2)..max_val_i+2),i=1..max_val_i)); Xi:='i'; Xp:=collect(simplify(expand(1-(y-1)*(y-w)*(y-i)*(y-i-1)/(w*i*(i+1)))),y); Xps:=unapply(subs(y^4=mom4,y^3=mom3,y^2=mom2,y=mom1,p),i); Xub:=min(seq(ps(i),i=2..min(max_val_i,w-2))); Xprintf(`LB=%.9f\n`,lb); Xprintf(`UB=%.9f\n`,ub); Xend; X X########################################### X# MC2 : 0 < pi/s < 1 and (w-1)^2-4*c > 0 X########################################### Xmc2 := proc() Xlocal MC2,pi,c,s; Xpi := mom1/w ; Xc :=(c2 - w*pi*(1-pi))/(2*pi*(1-pi)); Xs := (1/2)*(w + 1 - ( (w-1)^2 - 4*c)^(1/2)); XMC2 := 1 - (1-pi)*(1-pi/s)^(w-1); Xif pi/s > 1 or pi/s <=0 then ERROR(`Failed on 0 < pi/s < 1.`); fi; Xif (w-1)^2-4*c <=0 then ERROR(`Failed on (w-1)^2-4*c > 0`); fi; Xprintf(`MC2=%.9f\n`,MC2); Xend; X X############################### X# CPG2 : c2-m1>0 and 0<=r<1 X############################### Xcpg2:=proc(jmax) Xlocal r,i,j; Xif c2-mom1 <0 then ERROR(`Failed on c2-m1 >= 0`); fi; Xr:=((c2/c1)-1)/((c2/c1)+1); Xif r>1 or r< 0 then ERROR(`Failed on 0<= r < 1`); fi; Xp.0 := exp(- 2*mom1/(1 + c2/mom1)); Xl.1:=c1*(1-r)^2; Xfor j from 2 to jmax do l.j:=l.1*r^(j-1); od; Xfor j from 1 to jmax do X p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,m))],`+`)/j; Xod; Xprintf(`===CPG2===\n`); Xpge.0:=1; Xfor j from 1 to jmax do X pge.j := pge.(j-1)-p.(j-1); X printf(`%d,%.9f\n`,j,pge.j); Xod; Xend; X X################################################################# X# CPG4 approximation X# all lambdas should be positive and 0<=x<1 (r in the paper) . X################################################################## Xcpg4:= proc(jmax) Xlocal f,v,F,a,xi,i,j,soln,g1,g2,g3,g4,r; Xf:=x/(1-x); Xwith(linalg): Xv:=vector(4); Xv[1]:=simplify(x*diff(f,x)); Xfor j from 2 to 4 do v[j] :=simplify(x*diff(v[j-1],x)) od; XF:=matrix(4,3,[1,2,v[1],1,4,v[2],1,8,v[3],1,16,v[4]]); Xa:=matrix(3,1,[a0,a1,a2]); Xxi:=multiply(F,a); Xg1 := xi[1,1]=c1: Xg2 := xi[2,1]=c2: Xg3 := xi[3,1]=c3: Xg4 := xi[4,1]=c4: Xsoln:= fsolve({g1,g2,g3,g4},{a0,a1,a2,x}); Xl.1:=evalf(subs(soln,a0+a2*x)); Xif l.1 < 0 then ERROR(`Failed on lambda 1.`); fi; Xl.2:=evalf(subs(soln,a1+a2*x^2)); Xif l.2 < 0 then ERROR(`Failed on lambda 2.`); fi; Xl.3:=evalf(subs(soln,a2*x^3)); Xif l.3 < 0 then ERROR(`Failed on lambda 3.`); fi; Xr:=evalf(subs(soln,x)); Xif r < 0 or r > 1 then ERROR(`Failed on r.`); fi; Xp.0:=exp(-evalf(subs(soln,a0+a1+a2*x/(1-x)))); Xfor j from 4 to jmax do l.j:=l.3*r^(j-3); od; Xfor j from 1 to jmax do X p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,m))],`+`)/j; Xod; Xprintf(`===CPG4===\n`); Xpge.0:=1; Xfor j from 1 to jmax do X pge.j := pge.(j-1)-p.(j-1); X printf(`%d,%.9f\n`,j,pge.j); Xod; Xend; X X######################### X Xlpmin:=proc(q,deg,output) X########################################################################## X# This procedure is used by lp4. It computes the minimum probability p X# assuming a compound Poisson distribution for Y in which lambda[i]=L[i] X# = zero for i>q and the nonzero lambda's satisfy constraints X# of order deg. X# deg=0 assumes only that L[i]>=0 for all i. X# deg=1 assumes in addition that L[i]>=L[i+1] for all i. X# (in other words, the first differences (L[i]-L[i+1]) are >= 0.) X# deg=2 assumes in addition that the second differences are >=0. X########################################################################## Xglobal c1,c2,c3,c4; Xlocal kmom,i,j,e,constraint,func,L,A,T,tri,B,solmin,pmin; Xkmom:=4; XA:=array(1..(kmom+1),1..q,[seq([seq(j^i,j=1..q)],i=0..kmom)]); Xe:=array(1..q,1..1); Xif deg=0 then XT:=&*(); Xelse X# Create an upper triangular matrix of 1's (Is there a better way?) Xtri:=array(1..q,1..q,[seq([seq((1+signum(j-i+.5))/2,j=1..q)],i=1..q)]); XT:=evalm(tri^deg); Xfi; XB:=evalm(A &* T); Xconstraint:={seq(c.i=convert( [seq(B[i+1,j]*e[j,1],j=1..q)] ,`+`),i=1..kmom)}; Xfunc:=convert([seq(B[1,j]*e[j,1],j=1..q)],`+`); Xsolmin:=simplex[minimize](func,constraint,NONNEGATIVE); Xif nops(solmin)=0 then XRETURN(`Failed.`); Xfi; Xpmin:=1-evalf(subs(solmin,exp(-func))): Xif output=1 then Xassign(eval(solmin,1)); XL:=evalm(T &* e); Xprint(`Lambda's are`,L); Xfi; Xpmin; Xend; X X########### X### LP4 X########### Xlp4:=proc(output) X####################################################################### X# This procedure computes the LP4 approximation described in the paper. X# If output=1, the lambda's will be printed. X# It prints out the values of p, deg, and q. X# See further comments in the procedure lpmin. X####################################################################### Xlocal deg,p,qmax,qqqqq,qend; Xdeg := 2; Xqqqqq:=mmmmm+1; Xqmax:=max(2*mmmmm,10); Xwhile deg >= 0 do X p := lpmin(qqqqq,deg,output); X qend:=qqqqq; X if type(p,string) and qend <=qmax then X qqqqq:=qqqqq+5; X elif qend>qmax then X qqqqq:=mmmmm+1; X deg:=deg-1; X else X break X fi; Xod; Xif deg<0 then Xprintf(`No solution exists. Try a value of q > %d\n`,qend); Xelse X printf(`lp4=%.9f,deg=%d,q=%d\n`,p,deg,qqqqq); Xfi; Xend; X X########### X# GCP X########### Xgcp := proc(jmax) Xlocal pi,EY,mm,EI1I2,p,j,i; Xmm:=m-1: Xpi:=1-evalf(Sum('R(j,1)','j'=0..(mm-1))); XEY:=(n-mm)*pi; XEI1I2:=1-2*evalf(Sum('R(mm-2*j-1,1)','j'=0..floor((mm-1)/2))) X +(-1)^(mm+1)*R(0,2); Xp:=EI1I2/pi; Xfor j from 2 to mm do l.j := (n-mm)*pi*(1-p)^2*p^(j-1); od; Xl.1:=EY-evalf(Sum('j*l.j','j'=2..mm)); Xp.0:=exp(-evalf(Sum('l.j','j'=1..mm))); Xfor j from 1 to jmax do X p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,mm))],`+`)/j; Xod; Xprintf(`===GCP===\n`); Xpge.0:=1; Xfor j from 1 to jmax do X pge.j := pge.(j-1)-p.(j-1); X printf(`%d,%.9f\n`,j,pge.j); Xod; Xend; X X############### X# KLB X############### Xklb:=proc() Xlocal r,s2,lb; Xs2:=(mom2-mom1)/2; Xr:=floor(2*s2/mom1)+1; Xlb:=2*(r*mom1-s2)/(r*(r+1)); Xprintf(`klb=%.9f\n`,lb); Xend; X X############### X#### AP1-4 X############### Xap:= proc() X global c1,c2,c3,c4; X local i,f,A,B,C,l; X for i to 4 do X f := (j,k) -> k^j: X A := matrix(i,i,f); X B := inverse(A); X if i=1 then X C := multiply(B,[c1]); X elif i=2 then X C := multiply(B,[c1,c2]); X elif i=3 then X C := multiply(B,[c1,c2,c3]); X else X C := multiply(B,[c1,c2,c3,c4]); X fi; X print(C); X for l to i do X if C[l] < 0 then X l := 6; X fi; X od; X if(l<=5) then X print(`AP`,i,1-exp(-sum('C[k]','k'=1..i))); X else X print(`AP`,i); X fi; X od; Xend; X X X X X X ________This_Is_The_END________ if test `wc -l < maplecode` -ne 421; then echo 'shar: maplecode was damaged during transit (should have had 421 lines)' fi