#include <math.h>
#include<stdio.h>
#include<stdlib.h>
#define ITMAX 200
#define EPS 3.0e-7
#define FPMIN 1.0e-30

/*simple function*/
double dmin(double,double);
int imin(int,int);
int imax(int,int);
double power(double,int);
int ipower (int x, int i); ///x^i
/*==========================================*/
/* the subroutines to locate the memory of vector, matrix and array*/ 
double * dvector(int ,int );
int * ivector(int,int);
char * cvector(int, int);
double ** matrix( int, int, int,int);
char ** cmatrix( int, int, int,int);
int ** imatrix( int, int, int,int);
int ** ihmatrix( int, int, int *H);
double *** array(int,int, int, int, int,int);
int *** iarray(int,int, int, int, int,int);
char *** carray(int,int, int, int, int,int);

/*the subroutines to free the memory of vector, matrix and array*/
void free_ivector(int *,int,int);
void free_dvector(double *,int,int);
void free_matrix(double **,int , int , int , int );
void free_cvector(char *,int,int);
void free_cmatrix(char **,int , int , int , int );
void free_ihmatrix(int **, int, int);
void free_imatrix(int **, int,int,int,int);
void free_array(double ***,int,int,int , int , int , int );
void free_iarray(int ***,int,int,int , int , int , int );
void free_carray(char ***,int,int,int , int , int , int );
void nrerror(char *);
/*subroutines to generate random number===============================*/
double rand1(); /*uniform random number in interval (0,1].*/
double nrandom(double mean,double stdev ); /*normal random number */
double chisq(int dgfd);/*random number of chi_square distribution with dgree freedom dgfd*/
void multinomial(int notrial, double *theta, int *alpha, int length);
/*the random number of multinomial distribution output in alpha[] 
(alpha1,alpha2,..alpha_k)--M(n,theta1,..theta_k) length=k; notrial=n;*/
void   Dirichlet(double *theta, int *alpha, int length);
/*output of the random number of Dirichlet distribution in the vector theta[].  
(theta[1],theta[2],..theta[k])--Dirichlet(alpha1,..alpha_k) length=k*/
double tgamma(int,double);/*/random number of gamma distribution*/

/*===subroutines to calculate the distribution or the probability P(X<x)========*/
double pnormal(double x);  /*/distribution of standard normal =P(N(0,1)<=x)*/
double pchisq( int n,double x);/*/distribution of Chi-square distribution P(x^2<=x)*/

/* subroutines to sort the number========================================*/
void randomsort(int length, int* results );
void quicksort(double *vec, int leng);


/*==================auxiliary sunfunctions====================*/
double dennorm(double x,double mean ,double var);/*/ value of density of normal distribution                                                 //at point x*/ 
double gammp(double a,double x);     
void tred2(double **a, int n, double d[], double e[]);/*/ used by eigen()*/
double SIGN (double a,double b);/*/used by eigen();*/
void eigen( double **z, double d[], int n);
 /*/ a subrouting to calculate the eigen value output in vector d;
	//z as a input is a symmetric matrix; output a as a orthogonal matric
	//=the ith column is the eigen vector corresponding to eigen value d[i]*/
void tranmatrix(double **A,double **B,int n1,int n2);/*/input B(n1xn2), out put A(n2xn1)the teanport matrix of B*/
void   timematrix(double **D,double **A,double **B,int n1, int n2, int n);/*/D=A(n1xn)*B(nxn2)*/
double invertpdmatrix(double **oMatrix,double **dMatrix,int ndim);/*/dmatrix=inverse of (omatrix)-input matrix*/

/* special subfunctions====================================*/
void Arrayindex(int index, int *H, int length, int *numall);
 /* transfer an vector H to a integer index where
    length is the length of the vector
	numall is a vector, numall[i] is the total number of different values for H[i]
	when we consider H as a haplotype, numall[i] is the number of alleles at marker i
 */
int INdex(int *H, int length, int *numallele);
int skipline(FILE *fp, int num);

/*=============================================================
==================================================================*/
void Dirichlet(double *theta, int *alpha, int length)
{
	int i;
	double * H,temp;
	H=dvector(1,length);
	temp=0;
	for(i=1;i<=length;i++)
	{
		H[i]=chisq(2*alpha[i]);
		temp+=H[i];
	}
    for(i=1;i<=length;i++)
	{
		theta[i]=H[i]/temp;
	}
	free_dvector(H,1,length);

/*/the random number of Dirichlet distribution 
//(theta1,theta2,..theta_k)--Dirichlet(alpha1,..alpha_k)
//length=k
*/
}
/*/===================*/
void multinomial(int notrial, double *theta, int *alpha, int length)
{ 
	int i,j;
	double *H,temp;
	H=dvector(1,length);
	temp=0;
	for(i=1;i<=length;i++)
	{   
		alpha[i]=0;
		temp+=theta[i];
		H[i]=temp;
	}
	for(i=1;i<=notrial;i++)
	{  temp=rand1();
	    for(j=1;j<=length;j++)
		{ 
		   if((j==1)&&(temp<H[j]))alpha[j]++;
		   if((j>1)&&(temp>H[j-1])&&(temp<=H[j]))alpha[j]++;
		}
	}
free_dvector(H,1,length);



}








/*/======================*/
int imin(int i,int j)
{
	if(i<=j)return(i);
	else return(j);
}
/*/=============================*/
int imax(int i,int j)
{
	if(i>=j)return(i);
	else return(j);
}
/*/=============================*/
double dennorm(double x, double mean, double var)
{
	double a,pi;
    pi=3.1415926535;
	a=1/sqrt(2*pi*var)*exp(-(x-mean)*(x-mean)/2/var);
	return a;


}
/*/=======================================*/
double power(double x,int i)
{
	int j;
	double y;
	y=1;
	for(j=1;j<=i;j++)
	{
		y=y*x;
	}
	return y;
}
/*/=====================================*/
int ipower(int x,int i)
{
	int j;
	int y;
	y=1;
	for(j=1;j<=i;j++)
	{
		y=y*x;
	}
	return y;
}
/*/==================================*/
  double **matrix(int nrl, int nrh, int ncl, int nch)
/* int nrl,nrh,ncl,nch;*/

  {      

        int i;
        double **m;

        m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m -= nrl;

        for(i=nrl;i<=nrh;i++) {
                m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= ncl;
        }
        return m;
}
/*/===================================================*/
   int **ihmatrix(int nrl, int nrh, int*H)
/* int nrl,nrh,ncl,nch;*/

  {      

        int i,index;
        int **m;

        m=(int **) malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m -= nrl;
        index=0;
        for(i=nrl;i<=nrh;i++) {
			index++;
                m[i]=(int *) malloc((unsigned) (H[index])*sizeof(int));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= 1;
        }
        return m;
}

/*/===================================================*/

void free_matrix(double **m,int nrl, int nrh, int ncl, int nch)
/* double **m;
int nrl,nrh,ncl,nch;*/
{
        int i;

        for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
        free((char*) (m+nrl));
}

/*/=================================================*/
double ***array(int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
/* allocate a double 3tensor with range t[nrl...nrh][ncl...nch][ndl...ndh]*/
{
	int i,j,nrow=nrh-nrl+1, ncol = nch -ncl +1, ndep = ndh-ndl+1;
	double ***t;

	/* allocate pointers to pointers to rows*/
	t = ((double ***) malloc((size_t)(nrow+1)*sizeof(double **)));
	if (!t) nrerror("allocation failure in 1 in d3tensor()");
	t +=1;
	t -=nrl;

	/* allocate pointers to rows and set pointers to them*/
	t[nrl] = (double **) malloc((size_t)((nrow*ncol+1)*sizeof(double *)));
	if (!t[nrl]) nrerror("allocation failure 2 in d3tensor()");
	t[nrl] +=1;
	t[nrl] -=ncl;

	/* allocate rows and set pointers to them*/
    t[nrl][ncl] =(double *) malloc((size_t)((nrow*ncol*ndep+1)*sizeof(double)));
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in d3tensor()");
	t[nrl][ncl] +=1;
	t[nrl][ncl] -= ncl;

	for (j = ncl+1; j<=nch;j++)
		t[nrl][j] = t[nrl][j-1] + ndep;
		for(i=nrl+1;i<=nrh;i++){
			t[i]=t[i-1]+ncol;
			t[i][ncl]=t[i-1][ncl]+ncol*ndep;
			for(j=ncl+1;j<=nch;j++) t[i][j] = t[i][j-1]+ndep;
		}

		/* return pointer to array of pointers to rows*/
		return t;
	}

void free_array(double ***t, int nrl,int nrh,int ncl, int nch, int ndl, int ndh)
/* free a matrix allocated by d3tensor()*/
{
free((char *) (t[nrl][ncl]+ndl-1));
free((char *) (t[nrl]+ncl-1));
free((char *) (t+nrl -1));
}
/*/======================================*/
  char **cmatrix(int nrl, int nrh, int ncl, int nch)
/* int nrl,nrh,ncl,nch;*/

  {      

        int i;
        char **m;

        m=(char **) malloc((unsigned) (nrh-nrl+1)*sizeof(char*));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m -= nrl;

        for(i=nrl;i<=nrh;i++) {
                m[i]=(char *) malloc((unsigned) (nch-ncl+1)*sizeof(char));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= ncl;
        }
        return m;
}

/*/===================================================*/

void free_cmatrix(char **m,int nrl, int nrh, int ncl, int nch)
/* double **m;
int nrl,nrh,ncl,nch;*/
{
        int i;

        for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
        free((char*) (m+nrl));
}

/*/=================================================*/
char ***carray(int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
/* allocate a double 3tensor with range t[nrl...nrh][ncl...nch][ndl...ndh]*/
{
	int i,j,nrow=nrh-nrl+1, ncol = nch -ncl +1, ndep = ndh-ndl+1;
	char ***t;

	/* allocate pointers to pointers to rows*/
	t = ((char ***) malloc((size_t)(nrow+1)*sizeof(char **)));
	if (!t) nrerror("allocation failure in 1 in d3tensor()");
	t +=1;
	t -=nrl;

	/* allocate pointers to rows and set pointers to them*/
	t[nrl] = (char **) malloc((size_t)((nrow*ncol+1)*sizeof(char *)));
	if (!t[nrl]) nrerror("allocation failure 2 in d3tensor()");
	t[nrl] +=1;
	t[nrl] -=ncl;

	/* allocate rows and set pointers to them*/
    t[nrl][ncl] =(char *) malloc((size_t)((nrow*ncol*ndep+1)*sizeof(char)));
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in d3tensor()");
	t[nrl][ncl] +=1;
	t[nrl][ncl] -= ncl;

	for (j = ncl+1; j<=nch;j++)
		t[nrl][j] = t[nrl][j-1] + ndep;
		for(i=nrl+1;i<=nrh;i++){
			t[i]=t[i-1]+ncol;
			t[i][ncl]=t[i-1][ncl]+ncol*ndep;
			for(j=ncl+1;j<=nch;j++) t[i][j] = t[i][j-1]+ndep;
		}

		/* return pointer to array of pointers to rows*/
		return t;
	}

void free_carray(char ***t, int nrl,int nrh,int ncl, int nch, int ndl, int ndh)
/* free a matrix allocated by d3tensor()*/
{
free((char *) (t[nrl][ncl]+ndl-1));
free((char *) (t[nrl]+ncl-1));
free((char *) (t+nrl -1));
}
/*/====================================== */


void nrerror(char error_text[])
{
        void exit(int);

        fprintf(stderr,"Numerical Recipes run-time error...\n");
        fprintf(stderr,"%s\n",error_text);
        fprintf(stderr,"...now exiting to system...\n");
        exit(1);
}
/*/======================================================
//=================================================*/
int ***iarray(int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
/* allocate a double 3tensor with range t[nrl...nrh][ncl...nch][ndl...ndh]*/
{
	int i,j,nrow=nrh-nrl+1, ncol = nch -ncl +1, ndep = ndh-ndl+1;
	int ***t;

	/* allocate pointers to pointers to rows*/
	t = ((int ***) malloc((size_t)(nrow+1)*sizeof(int **)));
	if (!t) nrerror("allocation failure in 1 in d3tensor()");
	t +=1;
	t -=nrl;

	/* allocate pointers to rows and set pointers to them*/
	t[nrl] = (int **) malloc((size_t)((nrow*ncol+1)*sizeof(int *)));
	if (!t[nrl]) nrerror("allocation failure 2 in d3tensor()");
	t[nrl] +=1;
	t[nrl] -=ncl;

	/* allocate rows and set pointers to them*/
    t[nrl][ncl] =(int *) malloc((size_t)((nrow*ncol*ndep+1)*sizeof(int)));
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in d3tensor()");
	t[nrl][ncl] +=1;
	t[nrl][ncl] -= ncl;

	for (j = ncl+1; j<=nch;j++)
		t[nrl][j] = t[nrl][j-1] + ndep;
		for(i=nrl+1;i<=nrh;i++){
			t[i]=t[i-1]+ncol;
			t[i][ncl]=t[i-1][ncl]+ncol*ndep;
			for(j=ncl+1;j<=nch;j++) t[i][j] = t[i][j-1]+ndep;
		}

		/* return pointer to array of pointers to rows*/
		return t;
	}

void free_iarray(int ***t, int nrl,int nrh,int ncl, int nch, int ndl, int ndh)
/* free a matrix allocated by d3tensor()*/
{
free((char *) (t[nrl][ncl]+ndl-1));
free((char *) (t[nrl]+ncl-1));
free((char *) (t+nrl -1));
}
/*/======================================================*/

double *dvector(int n1,int n2)
{double *v;
v=(double*)malloc((size_t)((n2-n1+1)*sizeof(double)));
if(!v) nrerror("allocation failure in dvector()");
return v-n1;
}
/*/==================================================*/
int * ivector(int n1,int n2)
{int *v;
v=(int*)malloc((size_t)((n2-n1+1)*sizeof(int)));
if(!v) nrerror("allocation failure in dvector()");
return v-n1;
}
char *cvector(int n1,int n2)
{char *v;
v=(char*)malloc((size_t)((n2-n1+1)*sizeof(char)));
if(!v) nrerror("allocation failure in dvector()");
return v-n1;
}

/*/=================================================*/
void free_dvector(double *v, int n1,int n2)
{free((char*)(v+n1));
}
void free_cvector(char *v, int n1,int n2)
{free((char*)(v+n1));
}

void free_ivector(int *v, int n1,int n2)
{free((char*)( v+n1));
}

/*/==================================================*/
 int **imatrix(int nrl, int nrh, int ncl, int nch)
/* int nrl,nrh,ncl,nch;*/

{
        int i;
        int **m;

        m=(int **) malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m -= nrl;

        for(i=nrl;i<=nrh;i++) {
                m[i]=(int *) malloc((unsigned) (nch-ncl+1)*sizeof(int));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= ncl;
        }
        return m;
}

/*/===================================================*/

void free_imatrix(int **m,int nrl, int nrh, int ncl, int nch)
/* int **m;
int nrl,nrh,ncl,nch;*/
{
        int i;

        for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
        free((char*) (m+nrl));
}
/*/=============================*/
void free_ihmatrix(int **m,int nrl, int nrh)
/* int **m;
int nrl,nrh,ncl,nch;*/
{
        int i;

        for(i=nrh;i>=nrl;i--) free((char*) (m[i]+1));
        free((char*) (m+nrl));
}
/*//===============================================================
int random2(int m)
  
 { 
	
	float x; int z;
    x=rand1();
    z=(int)(m*x);
    return(z+1);
 }
*/

/*/===================================================*/
/*/subroutine to sort a vector;*/
void quicksort(double *vec, int leng)
{
	double a, *b;
	int i, k, begin, end;
	b=dvector(0,leng);
	if(leng<=1)return;
	if(leng==2)
	{
		if(vec[0]>vec[1])
		{
			a=vec[0]; 
			vec[0]=vec[1];
			vec[1]=a;
		}
		return;
	}
	k=(int)(rand1()*leng);
	for(i=0;i<leng;i++) b[i]=vec[i];
	begin=0;end=0;
    for(i=0;i<leng;i++)
	{
		if((i!=k)&&(b[i]>=b[k]))
		{
			vec[leng-1-end]=b[i];
			end++;
		}
		if((i!=k)&&(b[i]<b[k]))
		{
			vec[begin]=b[i];
			begin++;
		}
	}
	vec[begin]=b[k];
	if(begin==0)
		quicksort(vec+1,leng-1);
	else if(end==0)
		quicksort(vec,leng-1);
	else
	{
		quicksort(vec,begin);
		quicksort(vec+begin+1,end);
	}
	free_dvector(b,0,leng);
}



/*============================================================*/
/*=====================================================*/
void randomsort(int length,int *result)
{
	int i,j, k,*H,length1;
    length1=length;
	H=ivector(1,length);
	for(i=1;i<=length;i++)H[i]=i;
	for(i=1; i<=length;i++)
	{ j=(int)(rand1()*length1)+1;
	  result[i]=H[j];
	  for(k=j;k<length1;k++)
	  {
		  H[k]=H[k+1];
	  }
	  length1--;
	}
free_ivector(H,1,length);
}
/*/subroutine to sort a vector from minmum to maxmum;*/
void iquicksort(int *vec, int leng)
{
	int a, *b;
	int i, k, begin, end;
	b=ivector(0,leng);
	if(leng<=1)return;
	if(leng==2)
	{
		if(vec[0]>vec[1])
		{
			a=vec[0]; 
			vec[0]=vec[1];
			vec[1]=a;
		}
		return;
	}
	k=(int)(rand1()*leng);
	for(i=0;i<leng;i++) b[i]=vec[i];
	begin=0;end=0;
    for(i=0;i<leng;i++)
	{
		if((i!=k)&&(b[i]>=b[k]))
		{
			vec[leng-1-end]=b[i];
			end++;
		}
		if((i!=k)&&(b[i]<b[k]))
		{
			vec[begin]=b[i];
			begin++;
		}
	}
	vec[begin]=b[k];
	if(begin==0)
		iquicksort(vec+1,leng-1);
	else if(end==0)
		iquicksort(vec,leng-1);
	else
	{
		iquicksort(vec,begin);
		iquicksort(vec+begin+1,end);
	}
	free_ivector(b,0,leng);
}


/*/==============================================
//function to generate the uniformly distributed sample on (0,1) 
//long RX=1,RY=1,RZ=1;
 
 float rand1()
{  extern RX,RY,RZ; 
	
	 float t;
     int t1;
     RX=(171*RX)%30269;
     RY=(172*RY)%30309;
     RZ=(170*RZ)%30323;
     t=(float)(RX/30269.0+RY/30309.0+RZ/30323.0);
     t1=(int)(t);
     return(t-t1);
}

*/
/*/function to generate the uniformly distributed sample on (0,1) 
//long RX=1,RY=1,RZ=1;
 */
 double rand1()
{  extern RX,RY,RZ; 
	
	 double t;
     int t1;
     RX=(171*RX)%30269;
     RY=(172*RY)%30309;
     RZ=(170*RZ)%30323;
     t=(double)(RX/30269.0+RY/30309.0+RZ/30323.0);
     t1=(int)(t);
     return(t-t1);
}


/*============================================================*/
double nrandom(double mean,double stdev)
{
	double v1,v2,rsq,fac,r;
	do
	{
		v1=2.0*rand1()-1.0;
		v2=2.0*rand1()-1.0;
		rsq=v1*v1+v2*v2;
	}
	while (rsq>=1.0||rsq==0.0);
	fac=sqrt(-2.0*log(rsq)/rsq);
	r=v1*fac*stdev+mean;
	return(r);
/*/	r=v2*fac*stdev+mean;
//	return(r);
//  Generate the random number with univariate normal distributionN(mean,stdev);
//	Input
//		mean : the mean of distribution;
//		stdev : the standard deviation of distribution;
//	Output
//		a random number with univariate normal distribution;
//	Note
//		transformation method:
//		    y1=sqrt(-log(x1))cos(2*pi*x2);
//          y1=sqrt(-log(x1))sin(2*pi*x2);
*/
}


double chisq(int ndim)
{	int i;
	double tmp,sum;
	if(ndim<=0)
	{
		printf("The freedom of chi-square must be positive integer!\n");
		exit(1);
	}

	sum=0;
	for(i=0;i<ndim;i++)
	{
		tmp=nrandom(0,1);
		sum+=tmp*tmp;
	}
	return(sum);
/*/	generate the random number with chi-square distribution with freedom ndim;
//	Input
//		ndim : the freedom of chi-square distribution;
//	Return
//		the random number with chi-square distribution;
*/
}

double tgamma(int shape,double scale)
{
	int j;
	double am,e,s,v1,v2,x,y;

	if(shape<1)
		printf("shape parameter must be greater than 1\n");
	if(shape<6)
	{
		x=1.0;
		for(j=0;j<shape;j++)
			x=x*rand1();
		x=-log(x);
	}
	else
	{
		do
		{
			do
			{
				do
				{
					v1=rand1();
					v2=2.0*rand1()-1.0;
				}while(v1*v1+v2*v2>1.0);
				y=v2/v1;
				am=shape-1;
				s=sqrt(2.0*am+1.0);
				x=s*y+am;
			}while(x<=0.0);
			e=(1.0+y*y)*exp(am*log(x/am)-s*y);
		}while(rand1()>e);
	}
	return(x/scale);
/*/ 	Generate a random numbers with GAMA distribution;
//  the mean of disribution is shape/scale and the variance is shape/scale/scale;
//  Input 
//		shape : the shape parameter of gama distribution;
//		scale : the scale parameter of gama distribution;
//	Return
//		a random number with gama distribution;
*/
}



double dmin(double a1,double a2)
{
	if(a1>=a2)return(a2);
		else return(a1);
}

/*/======== in complete gamma function*/
double gammp(double a, double x)
{
    void gser(double *gamser, double a, double x, double *gln);
    void gcf(double *gammcf, double a,double x, double *gln);

       double gamser, gammcf, gln;
	   if(x<0.0||a<=0.0)
	   {
		   printf(" error for the input parameter in gammp function\n");
		   exit(1);
	   }
	   if(x<a+1.0)
	   {
		   gser(&gamser, a, x,&gln);/*/use series representation*/
           return gamser;
	   }
	   else
	   {
		   gcf(&gammcf,a,x,&gln);
		   return 1.0-gammcf;
	   }
}

/*/=======incomplete gamma fuction using series*/
 void gser(double *gamser, double a, double x, double *gln)
 {
	 double gammln(double xx);
	 int n;
	 double sum, del, ap;
	 *gln=gammln(a);
	 if(x<=0.0)
	 {
		 if(x<0.0){printf("error input parameter in routing gser\n"); exit(1);}
		 *gamser=0.0;
		 return;
	 }
	 else
	 {
		 ap=a;
		 del=sum=1.0/a;
		 for(n=1;n<=ITMAX;n++)
		 {
			 ++ap;
			 del *=x/ap;
			 sum+=del;
			 if(fabs(del)<fabs(sum)*EPS)
			 {
		
				 *gamser=sum*exp(-x+a*log(x)-(*gln));
			     return;
			 }
		 }
		 printf("a too large, ITMAX too small in routing gser\n");
		 return;
	 }
 }

 /*/====caculate incomplete gamma function using continue frection representation*/
 void gcf(double *gammcf, double a,double x, double *gln)
 {
	 double gammln(double xx);
      int i;
	  double an, b, c, d, del, h;
	  *gln=gammln(a);
	  b=x+1.0-a;
	  c=1.0/FPMIN;
	  d=1.0/b;
	  h=d;
	  for(i=1;i<=ITMAX;i++)
	  {
		  an=-i*(i-a);
		  b+=2.0;
		  d=an*d+b;
		  if(fabs(d)<FPMIN)d=FPMIN;
		  c=b+an/c;
		  if(fabs(c)<FPMIN)c=FPMIN;
		  d=1.0/d;
		  del=d*c;
		  h*=del;
		  if(fabs(del-1.0)<EPS)break;
	  }
	  if(i>ITMAX) printf("a too large, ITMAX too small in gcf\n");
	  *gammcf=exp(-x+a*log(x)-(*gln))*h;
 }
 /*/===log gamma function*/
 double gammln(double xx)
 {
	 double x,y,tmp, ser;
	 static double cof[6]={76.18009172947146, -86.50532032941677, 24.01409824083091,
		 -1.23173957240155, 0.1208650973866179e-2,-0.5395239384953e-5};
	 int j;
	 y=x=xx;
	 tmp=x+5.5;
	 tmp-=(x+0.5)*log(tmp);
	 ser=1.000000000190015;
     for(j=0;j<=5;j++)ser+=cof[j]/++y;
	 return -tmp+log(2.5066282746310005*ser/x);
 }

 /*/=======probbility of chisquare distribution with degree freedom n;*/
 double pchisq( int n,double x)
 {
	 return gammp((double)(n/2.0),x/2);
 }
 /*/=========distribution of standard normal distribution*/
 double pnormal(double x)
 {
	 double gammp(double a,double xx);
	 if(x>=0){return 0.5+0.5*gammp(0.5,x*x/2.0);}
	 else {return 0.5-0.5*gammp(0.5, x*x/2);}
 }
/*/=============== the following three subrouting for calculate the eigen value and eigenvector*/
 void tred2(double **a, int n, double d[], double e[])
 { 
	 /*/=this is a routing for transform a symmitric matrix to a tridiagonal matrix
	 //====a is*/ 
	int ll,k,j,i;
	double scale, hh, h,g,f;
	for(i=n;i>=2;i--)
	{
		ll=i-1;
		h=scale=0.0;
		if(ll>1)
		{
			for(k=1;k<=ll;k++)
				scale+=fabs(a[i][k]);
				if(scale==0.0) e[i]=a[i][ll];
				else
				{
					for(k=1;k<=ll;k++)
					{
						a[i][k]/=scale;
						h+=a[i][k]*a[i][k];
					}
					f=a[i][ll];
					g=(f>=0.0? -sqrt(h):sqrt(h));
					e[i]=scale*g;
					h-=f*g;
					a[i][ll]=f-g;
					f=0.0;
					for(j=1;j<=ll;j++)
					{
						a[j][i]=a[i][j]/h;
						g=0.0;
						for(k=1;k<=j;k++)
						  g+=a[j][k]*a[i][k];
						for(k=j+1;k<=ll;k++)
							g+=a[k][j]*a[i][k];
						e[j]=g/h;
						f+=e[j]*a[i][j];
					}
					hh=f/(h+h);
					for(j=1;j<=ll;j++)
					{
						f=a[i][j];
						e[j]=g=e[j]-hh*f;
						for(k=1;k<=j;k++)
						a[j][k]-=(f*e[k]+g*a[i][k]);
					}
				}
		}

			else
				e[i]=a[i][ll];
			d[i]=h;
	}
		d[1]=0.0;
		e[1]=0.0;
		for(i=1;i<=n;i++)
		{
			ll=i-1;
			if(d[i])
			{
				for(j=1;j<=ll;j++)
				{
					g=0.0;
					for(k=1;k<=ll;k++)
					   g+=a[i][k]*a[k][j];
					for(k=1;k<=ll;k++)
						a[k][j]-=g*a[k][i];
				}
			}
			d[i]=a[i][i];
			a[i][i]=1.0;
			for(j=1;j<=ll;j++)
				a[j][i]=a[i][j]=0.0;
		}
	}

/*/=====================================*/
void eigen( double **z, double d[], int n)
{ /*/ a subrouting to calculate the eigen value output in vector d;
	//z as a input is a symmetric matrix; output a as a orthogonal matric
	//= the ith column is the eigen vector corresponding to eigen value d[i]
	*/
	double SIGN (double a,double b);

	int m,ll,iter,k,i;
	double s,r,p,g,f,dd,c,b,*e;
	e=dvector(1,n);
    tred2(z,n,d,e);
	for(i=2;i<=n;i++) e[i-1]=e[i];
	e[n]=0.0;
	for(ll=1;ll<=n;ll++)
	{
		iter=0;
		do
		{
			for(m=ll;m<=n-1;m++)
			{
				dd=fabs(d[m])+fabs(d[m+1]);
				if((double)(fabs(e[m])+dd)==dd)break;
			}
			if(m!=1)
			{
				if(iter++==30) nrerror("Too many iterations in tqli");
				g=(d[ll+1]-d[ll])/(2.0*e[ll]);
				r=sqrt(g*g+1.0);
				g=d[m]-d[ll]+e[ll]/(g+SIGN(r,g));
				s=c=1.0;
				p=0.0;
				for(i=m-1;i>=ll;i--)
				{
					f=s*e[i];
					b=c*e[i];
					e[i+1]=(r=sqrt(f*f+g*g));
					if(r==0.0)
					{
						d[i+1]-=p;
						e[m]=0.0;
						break;
					}
					s=f/r;
					c=g/r;
					g=d[i+1]-p;
					r=(d[i]-g)*s+2.0*c*b;
					d[i+1]=g+(p=s*r);
					g=c*r-b;
					for(k=1;k<=n;k++)
					{
						f=z[k][i+1];
						z[k][i+1]=s*z[k][i]+c*f;
						z[k][i]=c*z[k][i]-s*f;
					}
				}
				if(r==0.0&&i>=1)continue;
				d[ll]-=p;
				e[ll]=g;
				e[m]=0.0;
			}
		}while(m!=ll);
	}
free_dvector(e,1,n);
}


double SIGN (double a,double b)
{
	if(b>=0)return fabs(a);
	else return -fabs(a);
}
/*/=====================
//====product of two matrix*/
 void timematrix(double **D,double **A,double **B,int n1, int n2, int n)
 {
	 int i,j,k;
	 for(i=1;i<=n1;i++)
	 {
		 for(j=1;j<=n2;j++)
		 { D[i][j]=0;
			 for(k=1;k<=n;k++)
			 {
				 D[i][j]+=A[i][k]*B[k][j];
			 }
		 }
	 }
 }
 /*/=================*/
void transmatrix(double **A,double **B,int n1,int n2)
{
	int i,j;
	for(i=1;i<=n1;i++)
	{
		for(j=1;j<=n2;j++)
		{
			A[j][i]=B[i][j];
		}
	}
}
/*/=============================*/
double invertpdmatrix(double **oMatrix,double **dMatrix,int ndim)
{
	int i,j,k;
	double dis=1,sum,*p;

	p=dvector(1,ndim);
	for(i=1;i<=ndim;i++)
	{
		for(j=i;j<=ndim;j++)
		{
			sum=oMatrix[i][j];
			for(k=i-1;k>=1;k--)
				sum-=oMatrix[i][k]*oMatrix[j][k];
			if(i==j)
			{
			
				dis*=sum;
				p[i]=sqrt(sum);
			}
			else
				oMatrix[j][i]=sum/p[i];
		}
	}
	for(i=1;i<=ndim;i++)
	{
		oMatrix[i][i]=1.0/p[i];
		for(j=i+1;j<=ndim;j++)
		{
			sum=0.0;
			for(k=i;k<j;k++)
				sum-=oMatrix[j][k]*oMatrix[k][i];
			oMatrix[j][i]=sum/p[j];
		}
	}
	for(i=1;i<=ndim;i++)
		for(j=i+1;j<=ndim;j++)
			oMatrix[i][j]=0;
	for(i=1;i<=ndim;i++)
		for(j=1;j<=ndim;j++)
		{
			dMatrix[i][j]=0;
			for(k=1;k<=ndim;k++)
				dMatrix[i][j]+=oMatrix[k][i]*oMatrix[k][j];
		}
	free_dvector(p,1,ndim);
	return(dis);
/*/  Get the inverse matrix of a positive definite matrix;
//	Input
//		**oMatrix : the input positive definite matrix;
//      ndim : the dimension of a matrix;
//  Output
//      **dMatrix : the inverse matrix;
//  Return
//		the determinant of the original matrix;
//  Note
//		the original matrix is changed;
*/
}
/*/=====================================*/
int skipline(FILE *fp, int num)
{
	int k=0;
	char ch;
	while(k<num)
	{
		for(;;)
		{
			ch=getc(fp);
			if(ch==EOF&&k<num)
				return(-k);
			if(ch=='\n')
			{
				k++;
				break;
			}
		}
	}
	return(num);
/*/  Skip the lines of a text file when read from or write to file.
//  Input
//		fp: the file pointer that can be read;
//		num: the number of lines that be skipped.
//	Return:
//		The function will return a negative integer that represents the total number
//		of this file when it is less than the number of lines that be skipped, otherwise
//		it will return the number of lines that be skipped.
*/
}
/*//========================
/==============================================
subfunction of transit from array index to number index
of haplotype 
=================================================*/
		
int INdex(int *H, int length, int *numallele)
{ long index, scale; int i;
 index=H[1]; scale=1;
  for(i=2;i<=length;i++)
  {scale=scale*numallele[i-1];
    index=index+scale*H[i];
  }
return index;
/*/====================
//return a number represent the array H, length--length of the vector H; 
//numballele is the maxmum number of allele in each marker (in each dimension of H) 
*/
}

/*====================================
function of transiting the number index to arrayindex*/

void Arrayindex(int index, int *H, int length, int *numall)
{
	int scale; int i; int A;
    A=1;
    for(i=1;i<=length;i++)A=A*numall[i];
	scale=A/numall[length];
    for(i=1; i<=length; i++)
	{
		H[length-i+1]=index/scale;
        index=index-H[length-i+1]*scale;
        if(i<length) scale=scale/numall[length-i];
	}
	/*/===================
	// invert transform of INdex( );
*/
}