#include <math.h>

/*constants for ran1(int *)*/
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 124456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

/*constants for ran11(long *)*/
#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)

/*constants for ran2(long *idum)*/
#define IM1 2147483563
#define IM2 2147483399
#define AM2 (1.0/IM1)
#define IMM1 (IM1-1)
#define IAA1 40014
#define IAA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB2 32
#define NDIV2 (1+IMM1/NTAB2)
#define EPS2 1.2e-7
#define RNMX2 (1.0-EPS2)

/*constants for ran4*/
#define NITER 4

void psdes(unsigned long *lword,unsigned long *irword);

void psdes(unsigned long *lword, unsigned long *irword)
{
	unsigned long i,ia, ib, iswap, itmph=0,itmpl=0;
	static unsigned long c1[NITER]=
		{0xbaa96887L,0x1e17d32cL,0x03bcdc3cL,0x0f33d1b2L};
	static unsigned long c2[NITER]=
		{0x4b0f3b58L,0xe874f0c3L,0x6955c5a6L,0x55a7ca46L};
		
	for(i=0;i<NITER;i++)
		{
		ia=(iswap=(*irword)) ^ c1[i];
		itmpl = ia & 0xffff;
		itmph = ia >> 16;
		ib=itmpl*itmpl+ ~(itmph*itmph);
		*irword=(*lword) ^ (((ia = (ib >> 16) |
			((ib & 0xffff) << 16)) ^ c2[i]) + itmpl*itmph);
		*lword=iswap;
		}
}

float ran4(long *idum)
{
	unsigned long irword,itemp,lword;
	static long idums = 0;

	static unsigned long jflone = 0x3f800000;
	static unsigned long jflmsk = 0x007fffff;
	
	if(*idum < 0)
		{
		idums = -(*idum);
		*idum=1;
		}
	irword=(*idum);
	lword=idums;
	psdes(&lword,&irword);
	itemp=jflone | (jflmsk & irword);
	++(*idum);
	return (*(float *)&itemp)-1.0;
}

float ran2(long *idum)
{
	int j;
	long k;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB2];
	float temp;
	
	if (*idum <= 0)
		{
		if (-(*idum)<1)
			*idum=1;
		else
			*idum= -(*idum);
		idum2=(*idum);
		for(j=NTAB2+7;j>=0;j--)
			{
			k=(*idum)/IQ1;
			*idum=IAA1*(*idum-k*IQ1)-k*IR1;
			if(*idum<0)
				*idum += IM1;
			if(j<NTAB2)
				iv[j] = *idum;
			}
		iy=iv[0];
		}
	k=(*idum)/IQ1;
	*idum=IAA1*(*idum-k*IQ1)-k*IR1;
	if(*idum<0)
		*idum += IM1;
	k=idum2/IQ2;
	idum2=IAA2*(idum2-k*IQ2)-k*IR2;
	if(idum2<0)
		idum2 += IM2;
	j=iy/NDIV2;
	iy=iv[j]-idum2;
	iv[j] = *idum;
	if(iy<1)
		iy += IMM1;
	if ((temp=AM2*iy)> RNMX2)
		return RNMX2;
	else
		return temp;
}

float ran1(int *idum)
{
	static long ix1,ix2,ix3;
	static float r[98];
	float temp;
	static int iff=0;
	int j;
	
	extern void nrerror(char []);

	if (*idum < 0 || iff == 0)
		{
		iff=1;
		ix1=(IC1-(*idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++)
			{
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
			}
		*idum=1;
		}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		nrerror("ran1: This cannot happen.");
	temp=r[j];
	r[j]=(ix1+ix2*RM2)*RM1;
	return temp;
}


float ran11(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float 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;
}

    
float gasdev(long *idum)
{
	static int iset=0;
	static double gset;
	float fac,r,v1,v2;
	float rangas();
    
	if  (iset == 0) {
		do {
			v1=2.0*ran11(idum)-1.0;
			v2=2.0*ran11(idum)-1.0;
			r=v1*v1+v2*v2;
		} while (r >= 1.0);
		fac=sqrt(-2.0*log(r)/r);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}