#include <math.h>
typedef double Number;

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

/**************************************************************/
/*  Returns an exponentially distributed, positive, random 
	deviate of unit mean, using random1 as the source of uniform
	deviates.
*/

Number expdev(long *idum)
{
	Number random(long *idum);
	Number dum;

	do
		dum=random(idum);
	while (dum == 0.0);
	return -log(dum);
}
/*     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   */
/**************************************************************/