/* some basic random number generation routines */
    
/* normally distributed random deviates */
#include <math.h>
#include "/data/xiaoxiak/thesis/header/nr.h"            
float gasdev(long *idum)
{
        float ran1(long *idum);
        static int iset=0;
        static double gset;
        double fac,rsq,v1,v2;
            
        if (*idum < 0) iset=0;
        if  (iset == 0) {
                do {
                        v1=2.0*ran1(idum)-1.0;
                        v2=2.0*ran1(idum)-1.0;
                        rsq=v1*v1+v2*v2;
                } while (rsq >= 1.0 || rsq == 0.0);
                fac=sqrt(-2.0*log(rsq)/rsq);
                gset=v1*fac;
                iset=1;
                return v2*fac;
        } else {
                iset=0;
                return gset;
        }
}


void multnorm(double *mean, double **var,int n,long *setSeed)
{
/* this draws a multivariate normal vector of dimension n
   and returns it in mean-so the mean is lost and lower half of var
   destroyed by choldc routine */
   int i,j;
   double *p,*mm;
   double sum;
   p=dvector(1,n);
   mm=dvector(1,n);
   choldc(var,n,p);
   for(i=1;i<=n;i++) mm[i]=gasdev(setSeed);
   for(i=1;i<=n;i++){
      sum=0.0;
      for(j=1;j<i;j++) sum += var[i][j]*mm[j];
      mean[i] += sum+p[i]*mm[i];
   }
   free_dvector(mm,1,n);
   free_dvector(p,1,n);
}

double gamdev(double ia, long *idum)
{
   double am,e,s,v1,v2,x,y,p,q;
   if(ia<=0.0) nrerror("error in gamdev: parameter non-positive");
   if(ia==1.0) x=-log(ran1(idum));
   if(ia<1.0){
      p=exp(1.0)/(ia+exp(1.0));
      do{
         if(ran1(idum)<p){
            x=pow(ran1(idum),1.0/ia);
            q=exp(-x);
         }
         else{
            x=1.0-log(ran1(idum));
            q=pow(x,ia-1.0);
         }
      } while(ran1(idum)>=q);
   }
   else{
      do{
         do{
            do{
               v1=ran1(idum);
               v2=2.0*ran1(idum)-1.0;
            } while (v1*v1+v2*v2 > 1.0);
            y=v2/v1;
            am=ia-1.0;
            s=sqrt(2.0*am+1.0); 
            x=s*y+am;
         } while(x <= 0);
         e=(1.0+y*y)*exp(am*log(x/am)-s*y);
      } while(ran1(idum) > e);
   }
   return x;
}


#define IA 16807                
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
         
float ran1(long *idum)
{
        int j;
        long k;
        static long iy=0;
        static long iv[NTAB];
        double temp;

        if (*idum <= 0 || !iy) {
                if (-(*idum) < 1) *idum=1;
                else *idum = -(*idum);
                for (j=NTAB+7;j>=0;j--) {
                        k=(*idum)/IQ;
                        *idum=IA*(*idum-k*IQ)-IR*k;
                        if (*idum < 0) *idum += IM;
                        if (j < NTAB) iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ;
        *idum=IA*(*idum-k*IQ)-IR*k;
        if (*idum < 0) *idum += IM;
        j=iy/NDIV;
        iy=iv[j];  
        iv[j] = *idum;
        if ((temp=AM*iy) > RNMX) return RNMX;
        else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX



/* random deviate by L'Ecuyer long period plus shuffle, uniform dist. */
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran2(long *idum)
{
	int j;
	long k;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if (*idum <= 0) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		idum2=(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ1;
			*idum=IA1*(*idum-k*IQ1)-k*IR1;
			if (*idum < 0) *idum += IM1;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ1;
	*idum=IA1*(*idum-k*IQ1)-k*IR1;
	if (*idum < 0) *idum += IM1;
	k=idum2/IQ2;
	idum2=IA2*(idum2-k*IQ2)-k*IR2;
	if (idum2 < 0) idum2 += IM2;
	j=iy/NDIV;
	iy=iv[j]-idum2;
	iv[j] = *idum;
	if (iy < 1) iy += IMM1;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX