#include <math.h>




#include <math.h>
typedef double Number;

/* tests use of exponential deviates  */
#include<iostream.h>
#include<math.h>
int main() {
	Number gasdev(long*);
	long n,i;
	n = -1;
	for(i=1;i<20;i++){
		cout << gasdev(&n) << endl;
		n = 2*n;
	}
	return 0;
}

/**************************************************************/
/*  Returns a normally distributed deviate with zero mean and 
	unit variance, using random1 as the source of uniform
	deviates.
*/

Number gasdev(long *idum)
{
	Number random(long *idum);
	static int iset=0;
	static Number gset;
	Number fac,rsq,v1,v2;

	if  (iset == 0) {					/* We don't have an extra deviate handy, so */
		do {
			v1=2.0*random(idum)-1.0;    /* pick two uniform numbers in the square */
			v2=2.0*random(idum)-1.0;    /* extending from -1 to +1 in each direction */
			rsq=v1*v1+v2*v2;			/* See if they are in the unit circle */
		} while (rsq >= 1.0 || rsq == 0.0);  /* If they are not, try again */
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;							/* Set flag  */
		return v2*fac;
	} else {							/* We have an extra devate handy, */
		iset=0;							/* so unset the flag */
		return gset;					/* and return it. */
	}
}

/*     End of expdev routine   */
/**************************************************************/
/**************************************************************/
/*  routine random   */                                           
/*  Generates random numbers between 0.0 and 1.0. Call with idum
	a negative number to initialize; thereafter, do not alter
	idum between successive deviates in a sequence. RNMX 
	approximates the largest Numbering value that is less than 1.\
	The period of this routine is ~ 10^8
*/
#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)

Number random(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	Number temp;

	if (*idum <= 0 || !iy) {			/* initialize  */
		if (-(*idum) < 1) *idum=1;		/* Be sure to prevent idum=0 */
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {		/* load the shuffle table after 8 warm-ups) */
			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;						/* Start here when not initializing  */
	*idum=IA*(*idum-k*IQ)-IR*k;			/* Compute idum=(IA*idum) % IM without */
	if (*idum < 0) *idum += IM;			/* overflows by Schrage's method */
	j=iy/NDIV;							/* Will be in the range 0..NTAB-1. */
	iy=iv[j];							/* Output previously stored value and */
	iv[j] = *idum;						/* refill the suffle table  */ 
	if ((temp=AM*iy) > RNMX) return RNMX; /* Because users don't expect endpoint values */
	else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

/*     End of random routine   */
/**************************************************************/