#include <math.h>
#include<stdio.h>
#include<stdlib.h>
#define ITMAX 200
#define EPS 3.0e-7
#define FPMIN 1.0e-30
#define PI               3.1415926535897
#define CONCRIT                 10.0e-12
#define MAX_TRAP                       2
#define SLOPE_SCALE_FACT             1.0
#define INT_SCALE_FACT               0.5
#define TINY 1.0e-20
static double maxarg1, maxarg2;
#define FMAX(a,b) (maxarg1=(a), maxarg2=(b), (maxarg1) > (maxarg2) ?\
             (maxarg1) : (maxarg2))
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static int iminarg1, iminarg2;
#define IMIN(a, b) (iminarg1=(a), iminarg2=(b), (iminarg1) < (iminarg2) ?\
		    (iminarg1) : (iminarg2))
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;


void Arrayindex(int index, int *H, int length, int *numallele);
int INdex(int *H, int length, int *numallele);

double dmin(double,double);
int imin(int,int);
int imax(int,int);
double power(double,int);
int ipower (int x, int i);
double rand1();
void randomsort(int, int* );
void quicksort(double *vec, int leng);
void iquicksort(int *vec, int leng);
void quicksortsort(double *vec, int *index, int leng);
//befor using quicksorsort, index[i]=i; sort from smallest to largest;
void Lquicksortsort(double *vec, int *index, int leng);
//sort from largest to smallest

void Quicksort(double *arr, int n);
void Quickindex(double *arr, int *indx, int n);
void iQuicksort(int *arr, int n);
void iQuickindex(int *arr, int *indx, int n);
//output indx[i] such that arr[indx[i]] is in ascending order
double ** matrix( int, int, int,int);
char ** cmatrix( int, int, int,int);
char * cvector(int, int);
int ** imatrix( int, int, int,int);
int ** ihmatrix( int, int, int *H);

double *** array(int,int, int, int, int,int);
void free_array(double ***,int,int,int , int , int , int );
int *** iarray(int,int, int, int, int,int);
void free_iarray(int ***,int,int,int , int , int , int );
char *** carray(int,int, int, int, int,int);
void free_carray(char ***,int,int,int , int , int , int );


double * dvector(int ,int );
int * ivector(int,int);
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 nrerror(char *);

void choldc(double **, int, double []);
void cholsl(double **, int, double [], double [], double []);
void lubksb(double **, int, int *, double *);
void ludcmp(double **, int, int *, double *);
double pythag(double, double);
void svdcmp(double **a, int m, int n, double w[], double **v);
/*
Given a matrix a[1..m][1..n], this routine computes its singular value
decomposition, A = U.W.V' The matrix U replaces a on output. The
diagonal matrix of singular values W is output as a vector w[1..n]. The
matrix V (not the transpose V') is output as v[1..n][1..n].

I modified the routine svdcmp to function in double precision.
*/
double nrandom(double mean,double stdev);
//  Generate the random number with univariate normal distributionN(mean,stdev);

double chisq(int ndim );//random number of chi_square distribution with digree fredom ndim
double llgamma(int,double);//random number of gamma distribution
double dennorm(double x,double mean ,double var);// value of density of normal distribution
                                                 //at point x 
double gammp(double a,double 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)
int rpoisson(double lambda,int cutoff);
//generate the poisson random number with mean lambda and largest number=cutoff

void   Dirichlet(double *theta, int *alpha, int length);
//out put of the random number of Dirichlet distribution in the vector theta[].  
//(theta[1],theta[2],..theta[k])--Dirichlet(alpha1,..alpha_k)
//length=k
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 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 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
//vo
int skipline(FILE *fp, int num);
int compvec(int *H1, int *H2, int leng);
void cpvec(int *H1, int *H2, int leng);//copy H2 to H1

//=====================
void cpvec(int *H1, int *H2, int leng)
{
	int i;
	for(i=1;i<=leng;i++)H1[i]=H2[i];
}

int compvec(int *H1, int *H2, int leng)
{
	int i;
	for(i=1;i<=leng;i++)
	{
		if(H1[i]!=H2[i])return 0;
	}
	return 1;
}
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( (m[i]+ncl));
        free( (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 cvector()");
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)
	{
		free_dvector(b,0,leng);return;
	}
	if(leng==2)
	{
		if(vec[0]>vec[1])
		{
			a=vec[0]; 
			vec[0]=vec[1];
			vec[1]=a;
		}
        free_dvector(b,0,leng);
		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 quicksortsort(double *vec, int *index, int leng)
{
	double a, *b;
	int i, k, begin,indi, end, *H;
	b=dvector(1,leng);
	H=ivector(1,leng);
	if(leng<=1)
	{ 
		free_dvector(b,1,leng);
	    free_ivector(H,1,leng);
		return;
	}
	if(leng==2)
	{
		if(vec[1]>vec[2])
		{
			a=vec[1]; 
			vec[1]=vec[2];
			vec[2]=a;
			indi=index[1];
			index[1]=index[2];
			index[2]=indi;
		}
		free_dvector(b,1,leng);
	    free_ivector(H,1,leng);
		return;
	}
	k=(int)(rand1()*leng)+1;
	for(i=1;i<=leng;i++) {b[i]=vec[i];H[i]=index[i];}
	begin=1;end=0;
    for(i=1;i<=leng;i++)
	{
		if((i!=k)&&(b[i]>=b[k]))
		{
			vec[leng-end]=b[i];
			index[leng-end]=H[i];
			end++;
		}
		if((i!=k)&&(b[i]<b[k]))
		{
			vec[begin]=b[i];
			index[begin]=H[i];
			begin++;
		}
	}
	vec[begin]=b[k];
	index[begin]=H[k];
	if(begin==1)
		quicksortsort(vec+1,index+1,leng-1);
	else if(end==0)
		quicksortsort(vec,index,leng-1);
	else
	{
		quicksortsort(vec,index,begin-1);
		
		quicksortsort(vec+begin,index+begin,end);
	}
	free_dvector(b,1,leng);
	free_ivector(H,1,leng);
	//=======vec as intput is a double vector, index is the vector with index[i]=i;
	//====as output vec is sort from smallest to largest of input vec
	//==index[i] is the original order of vec[i];
}

/*============================================================*/
//============================================================
void Lquicksortsort(double *vec, int *index, int leng)
{
	double a, *b;
	int i, k, begin,indi, end, *H;
	b=dvector(1,leng);
	H=ivector(1,leng);
	if(leng<=1)
	{
		free_dvector(b,1,leng);
		free_ivector(H,1,leng);
		return;
	}
	if(leng==2)
	{
		if(vec[1]<vec[2])
		{
			a=vec[1]; 
			vec[1]=vec[2];
			vec[2]=a;
			indi=index[1];
			index[1]=index[2];
			index[2]=indi;
		}
		free_dvector(b,1,leng);
		free_ivector(H,1,leng);
		return;
	}
	k=(int)(rand1()*leng)+1;
	for(i=1;i<=leng;i++) {b[i]=vec[i];H[i]=index[i];}
	begin=1;end=0;
    for(i=1;i<=leng;i++)
	{
		if((i!=k)&&(b[i]<=b[k]))
		{
			vec[leng-end]=b[i];
			index[leng-end]=H[i];
			end++;
		}
		if((i!=k)&&(b[i]>b[k]))
		{
			vec[begin]=b[i];
			index[begin]=H[i];
			begin++;
		}
	}
	vec[begin]=b[k];
	index[begin]=H[k];
	if(begin==1)
		Lquicksortsort(vec+1,index+1,leng-1);
	else if(end==0)
		Lquicksortsort(vec,index,leng-1);
	else
	{
		Lquicksortsort(vec,index,begin-1);
		
		Lquicksortsort(vec+begin,index+begin,end);
	}
	free_dvector(b,1,leng);
	free_ivector(H,1,leng);
	//=======vec as intput is a double vector, index is the vector with index[i]=i;
	//====as output vec is sort from smallest to largest of input vec
	//==index[i] is the original order of vec[i];
}

/*============================================================*/
//== subroutine
void Quicksort(double *arr, int n)
{
	int i,ir=n, j,k,l=1, *istack;
	int jstack=0, NSTACK, M;
	double a, temp;

	NSTACK=50;
	M=7;
	istack=ivector(1,NSTACK);
	for(;;)
	{
		if(ir-l<M)
		{
			for(j=l+1;j<=ir;j++)
			{
				a=arr[j];
				for(i=j-1;i>=l;i--)
				{
					if(arr[i]<=a)break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if(jstack==0)break;
			ir=istack[jstack--];
			l=istack[jstack--];
		}
		else
		{
			k=(l+ir)>>1;
			SWAP(arr[k],arr[l+1])
				if(arr[l]>arr[ir])
				{
					SWAP(arr[l],arr[ir])
				}
				if(arr[l+1]>arr[ir])
				{
					SWAP(arr[l+1],arr[ir])
				}
				if(arr[l]>arr[l+1])
				{
					SWAP(arr[l],arr[l+1])
				}
				i=l+1;
				j=ir;
				a=arr[l+1];
				for(;;)
				{
					do i++; while (arr[i]<a);
					do j--; while (arr[j]>a);
					if(j<i)break;
					SWAP(arr[i],arr[j]);
				}
				arr[l+1]=arr[j];
				arr[j]=a;
				jstack+=2;
				if(jstack>NSTACK) printf("NSTACK too small in Quicksort\n");
				if(ir-i+1>=j-l)
				{
					istack[jstack]=ir;
					istack[jstack-1]=i;
					ir=j-1;
				}
				else
				{
					istack[jstack]=j-1;
					istack[jstack-1]=l;
					l=i;
				}
		}
	}
	free_ivector(istack,1,NSTACK);
}
//=========================================================
/*============================================================*/
//== subroutine
void iQuicksort(int *arr, int n)
{
	int i,ir=n, j,k,l=1, *istack;
	int jstack=0, NSTACK, M;
	int a, temp;

	NSTACK=50;
	M=7;
	istack=ivector(1,NSTACK);
	for(;;)
	{
		if(ir-l<M)
		{
			for(j=l+1;j<=ir;j++)
			{
				a=arr[j];
				for(i=j-1;i>=l;i--)
				{
					if(arr[i]<=a)break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if(jstack==0)break;
			ir=istack[jstack--];
			l=istack[jstack--];
		}
		else
		{
			k=(l+ir)>>1;
			SWAP(arr[k],arr[l+1])
				if(arr[l]>arr[ir])
				{
					SWAP(arr[l],arr[ir])
				}
				if(arr[l+1]>arr[ir])
				{
					SWAP(arr[l+1],arr[ir])
				}
				if(arr[l]>arr[l+1])
				{
					SWAP(arr[l],arr[l+1])
				}
				i=l+1;
				j=ir;
				a=arr[l+1];
				for(;;)
				{
					do i++; while (arr[i]<a);
					do j--; while (arr[j]>a);
					if(j<i)break;
					SWAP(arr[i],arr[j]);
				}
				arr[l+1]=arr[j];
				arr[j]=a;
				jstack+=2;
				if(jstack>NSTACK) printf("NSTACK too small in Quicksort\n");
				if(ir-i+1>=j-l)
				{
					istack[jstack]=ir;
					istack[jstack-1]=i;
					ir=j-1;
				}
				else
				{
					istack[jstack]=j-1;
					istack[jstack-1]=l;
					l=i;
				}
		}
	}
	free_ivector(istack,1,NSTACK);
}
//===subroutine to get the index of the sort. output indx[i] such that arr[indx[i]] is in ascending order
void Quickindex(double *arr, int *indx, int n)
{
	int i,indxt,ir=n,temp, j,k,l=1, *istack;
	int jstack=0, NSTACK, M;
	double a;

	NSTACK=50;
	M=7;
	istack=ivector(1,NSTACK);
	for(j=1;j<=n;j++)indx[j]=j;
	for(;;)
	{
		if(ir-l<M)
		{
			for(j=l+1;j<=ir;j++)
			{
				indxt=indx[j];
				a=arr[indxt];
				for(i=j-1;i>=l;i--)
				{
					if(arr[indx[i]]<=a)break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if(jstack==0)break;
			ir=istack[jstack--];
			l=istack[jstack--];
		}
		else
		{
			k=(l+ir)>>1;
			SWAP(indx[k],indx[l+1]);
				if(arr[indx[l]]>arr[indx[ir]])
				{
					SWAP(indx[l],indx[ir])
				}
				if(arr[indx[l+1]]>arr[indx[ir]])
				{
					SWAP(indx[l+1],indx[ir])
				}
				if(arr[indx[l]]>arr[indx[l+1]])
				{
					SWAP(indx[l],indx[l+1])
				}
				i=l+1;
				j=ir;
				indxt=indx[l+1];
				a=arr[indxt];
				for(;;)
				{
					do i++; while (arr[indx[i]]<a);
					do j--; while (arr[indx[j]]>a);
					if(j<i)break;
					SWAP(indx[i],indx[j]);
				}
				indx[l+1]=indx[j];
				indx[j]=indxt;
				jstack+=2;
				if(jstack>NSTACK) printf("NSTACK too small in Quicksort\n");
				if(ir-i+1>=j-l)
				{
					istack[jstack]=ir;
					istack[jstack-1]=i;
					ir=j-1;
				}
				else
				{
					istack[jstack]=j-1;
					istack[jstack-1]=l;
					l=i;
				}
		}
	}
	free_ivector(istack,1,NSTACK);
}
void iQuickindex(int *arr, int *indx, int n)
{
	int i,indxt,ir=n,temp, j,k,l=1, *istack;
	int jstack=0, NSTACK, M;
	int a;

	NSTACK=50;
	M=7;
	istack=ivector(1,NSTACK);
	for(j=1;j<=n;j++)indx[j]=j;
	for(;;)
	{
		if(ir-l<M)
		{
			for(j=l+1;j<=ir;j++)
			{
				indxt=indx[j];
				a=arr[indxt];
				for(i=j-1;i>=l;i--)
				{
					if(arr[indx[i]]<=a)break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if(jstack==0)break;
			ir=istack[jstack--];
			l=istack[jstack--];
		}
		else
		{
			k=(l+ir)>>1;
			SWAP(indx[k],indx[l+1]);
				if(arr[indx[l]]>arr[indx[ir]])
				{
					SWAP(indx[l],indx[ir])
				}
				if(arr[indx[l+1]]>arr[indx[ir]])
				{
					SWAP(indx[l+1],indx[ir])
				}
				if(arr[indx[l]]>arr[indx[l+1]])
				{
					SWAP(indx[l],indx[l+1])
				}
				i=l+1;
				j=ir;
				indxt=indx[l+1];
				a=arr[indxt];
				for(;;)
				{
					do i++; while (arr[indx[i]]<a);
					do j--; while (arr[indx[j]]>a);
					if(j<i)break;
					SWAP(indx[i],indx[j]);
				}
				indx[l+1]=indx[j];
				indx[j]=indxt;
				jstack+=2;
				if(jstack>NSTACK) printf("NSTACK too small in Quicksort\n");
				if(ir-i+1>=j-l)
				{
					istack[jstack]=ir;
					istack[jstack-1]=i;
					ir=j-1;
				}
				else
				{
					istack[jstack]=j-1;
					istack[jstack-1]=l;
					l=i;
				}
		}
	}
	free_ivector(istack,1,NSTACK);
}
/*=====================================================*/
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 llgamma(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);//pythag(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)
			{
				//if(sum<=0)
					//cerr <<" Cholesky Decomposition Failed\n";
				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 *numallele)
//long index;int * H; 
{
	int scale; int i; int A;
    A=1;
    for(i=1;i<=length;i++)A=A*numallele[i];
	scale=A/numallele[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/numallele[length-i];
	}
	//===================
	// invert transform of INdex( );
}
//==========SVD==============
double pythag(double a, double b)
/* pythag computes sqrt(a^2 + b^2) without destructive underflow or
   overflow.
*/
{
  double absa, absb;

  absa = fabs(a);
  absb = fabs(b);
  if (absa > absb)
    return absa*sqrt(1.0+DSQR(absb/absa));
  else
    return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
}

void svdcmp(double **a, int m, int n, double w[], double **v)
{
/*
Given a matrix a[1..m][1..n], this routine computes its singular value
decomposition, A = U.W.V' The matrix U replaces a on output. The
diagonal matrix of singular values W is output as a vector w[1..n]. The
matrix V (not the transpose V') is output as v[1..n][1..n].

I modified the routine svdcmp to function in double precision.
*/

  int flag, i, its, j, jj, k, l, nm;
  double anorm, c, f, g, h, s, scale, x, y, z, *rv1;

  rv1 = dvector(1, n);
  g = scale = anorm = 0.0;

  /* Householder reduction to bidiagonal form */

  for (i = 1; i <= n; i++) {
    l = i + 1;
    rv1[i] = scale*g;
    g = s = scale = 0.0;
    if (i <= m) {
      for (k = i; k <= m; k++) scale += fabs(a[k][i]);
      if (scale) {
        for (k = i; k <= m; k++) {
          a[k][i] /= scale;
          s += a[k][i]*a[k][i];
        }
        f = a[i][i];
        g = -SIGN(sqrt(s),f);
        h = f*g-s;
	a[i][i] = f-g;
        for (j = l; j <= n; j++) {
            for (s = 0.0, k = i; k <= m; k++) s += a[k][i]*a[k][j];
            f = s/h;
            for (k = i; k <= m; k++) a[k][j] += f*a[k][i];
	}
        for (k = i; k <= m; k++) a[k][i] *= scale;
      }
    }
    w[i] = scale*g;
    g = s = scale = 0.0;
    if (i <= m && i != n) {
      for (k = l; k <= n; k++) scale += fabs(a[i][k]);
      if (scale) {
        for (k = l; k <= n; k++) {
          a[i][k] /= scale;
          s += a[i][k]*a[i][k];
        }
        f = a[i][l];
        g = -SIGN(sqrt(s),f);
        h = f*g-s;
        a[i][l] = f-g;
        for (k = l; k <= n; k++) rv1[k] = a[i][k]/h;
        for (j = l; j <= m; j++) {
            for (s = 0.0, k = l; k <= n; k++) s += a[j][k]*a[i][k];
            for (k = l; k <= n; k++) a[j][k] += s*rv1[k];
        }
        for (k = l; k <= n; k++) a[i][k] *= scale;
      }
    }
    anorm = FMAX(anorm, (fabs(w[i])+fabs(rv1[i])));
  }

  /* Accumulation of right-hand transformations. */
  for (i = n; i >= 1; i--) {
    if (i < n) {
      if (g) {
        for (j = l; j <= n; j++)   /* Double division to avoid underflow */
          v[j][i] = (a[i][j]/a[i][l])/g;
        for (j = l; j <= n; j++) {
          for (s = 0.0, k = l; k <= n; k++) s += a[i][k]*v[k][j];
          for (k = l; k <= n; k++) v[k][j] += s*v[k][i];
        }
      }
      for (j = l; j <= n; j++) v[i][j] = v[j][i] = 0.0;
    }
    v[i][i] = 1.0;
    g = rv1[i];
    l = i;
  }

  /* Accumulation of left-hand transformations. */
  for (i = IMIN(m,n); i >= 1; i--) {
    l = i+1;
    g = w[i];
    for (j = l; j <= n; j++) a[i][j] = 0.0;
    if (g) {
      g = 1.0/g;
      for (j = l; j <= n; j++) {
        for (s = 0.0, k = l; k <= m; k++) s += a[k][i]*a[k][j];
          f = (s/a[i][i])*g;
          for (k = i; k <= m; k++) a[k][j] += f*a[k][i];
      }
      for (j = i; j <= m; j++) a[j][i] *= g;
    } else for (j = i; j <= m; j++) a[j][i] = 0.0;
    ++a[i][i];
  }

  /* Diagonalization of the bidagonal form: Loop over singular
     values, and over allowed iterations.
  */
  for (k = n; k >= 1; k--) {
    for (its = 1; its <= 30; its++) {
      flag = 1;
      for (l = k; l >= 1; l--) {        /* Test for splitting */
	nm = l-1;                       /* rv1[1] is always zero. */
        if ((double)(fabs(rv1[l])+anorm) == anorm) {
          flag = 0;
          break;
        }
        if ((double)(fabs(w[nm])+anorm) == anorm) break;
      }
      if (flag) {
	c = 0.0;                        /* Cancellation of rv1[l], if l > 1 */
        s = 1.0;
        for (i = l; i <= k; i++) {
          f = s*rv1[i];
          rv1[i] = c*rv1[i];
          if ((double)(fabs(f)+anorm) == anorm) break;
          g = w[i];
          h = pythag(f, g);
          w[i] = h;
          h = 1.0/h;
          c = g*h;
          s = -f*h;
          for (j = 1; j <= m; j++) {
            y = a[j][nm];
            z = a[j][i];
            a[j][nm] = y*c+z*s;
            a[j][i] = z*c-y*s;
          }
        }
      }
      z = w[k];
      if (l == k) {                 /* Convergence. */
        if (z < 0.0) {              /* Sing. val. is made non-negative */
          w[k] = -z;
          for (j = 1; j <= n; j++) v[j][k] = -v[j][k];
        }
        break;
      }
      //if (its == 30) nrerror("No convergence in 30 svdcmp iterations");
      x = w[l];                     /* Shift from bottom 2-by-2 minor. */
      nm = k-1;
      y = w[nm];
      g = rv1[nm];
      h = rv1[k];
      f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
      g = pythag(f, 1.0);
      f = ((x-z)*(x+z)+h*((y/(f+SIGN(g, f)))-h))/x;

      /* Next QR transformation. */
      c = s = 1.0;
      for (j = l; j <= nm; j++) {
        i = j + 1;
        g = rv1[i];
        y = w[i];
        h = s*g;
        g = c*g;
        z = pythag(f, h);
        rv1[j] = z;
        c = f/z;
        s = h/z;
        f = x*c+g*s;
        g = g*c-x*s;
        h = y*s;
        y *= c;
        for (jj = 1; jj <= n; jj++) {
          x = v[jj][j];
          z = v[jj][i];
          v[jj][j] = x*c+z*s;
          v[jj][i] = z*c-x*s;
        }
        z = pythag(f, h);
        w[j] = z;                  /* Rotation can be arbitrary if z = 0. */
        if (z) {
          z = 1.0/z;
          c = f*z;
          s = h*z;
        }
        f = c*g+s*y;
        x = c*y-s*g;
        for (jj = 1; jj <= m; jj++) {
          y = a[jj][j];
          z = a[jj][i];
          a[jj][j] = y*c+z*s;
          a[jj][i] = z*c-y*s;
        }
      }
      rv1[l] = 0.0;
      rv1[k] = f;
      w[k] = x;
    }
  }
  free_dvector(rv1, 1, n);
}

//generate the poisson random number with mean lambda and largest number=cutoff
int rpoisson(double lambda,int cutoff)
{
	double g,t;
	int i,em;
	em=-1;
	g=exp(-lambda);
	t=1.0;
	for(i=0;i<=cutoff;i++)
	{
		em++;
		t*=rand1();
		if(t<=g)break;
	}
	return em;
}