// A program to calculate by random selection, the coefficients of
// addition and subtraction for Gaussian variables.

#include <stdio.h>
#include <math.h>

FILE *randomizer;

#define mask (8192-1) 
unsigned char mixArray[mask+1];

unsigned char RandomByte(void){
	static unsigned long long seed = -1;
	static int p0=0, p1=0;
	unsigned char result, swap;
	
	p1 = (p1 + mixArray[p1]+2) & mask;
	result = mixArray[p1];

	seed += (seed>>16);
	p0 = seed & mask;
	swap  =mixArray[ p0 ];

	mixArray[p0] = result;
	mixArray[p1] = swap;
	
	seed -= ((unsigned long long)swap)<<48;
	return result;
}

static void InitRandom(void) {
	int i,o;
	unsigned char swapme[mask+1];
	unsigned char a,b;

	for(i=0; i<sizeof(mixArray); ++i) mixArray[i]=i&0xff;
	fread( (unsigned char*)(void*)&swapme, sizeof(swapme), 1, randomizer );
	for(i=0; i<sizeof(mixArray); ++i) {
		o = swapme[i];
		a = mixArray[i]; b=mixArray[ o ];
		mixArray[o]=a; mixArray[i]=b;
	}
	for(i=0; i<1000000; ++i) RandomByte();
}


double Random01(){
	static double adjust = 1.0/(1ULL<<5*8);
	static unsigned long long result=0;
	static double homeing = 0.5; 
	double res;
	unsigned char* p = (unsigned char*)(void*)&result;	
	int i;
	
	for( i=0; i<5; ++i) *(p++) = RandomByte();
	res=(double)result * adjust;
	if (RandomByte()==127) { // Attempt to remove long term biases.
		if (homeing>0.5) 
			while (res>0.5) res=Random01();
		else
			while (res<0.5) res=Random01();	
	}
	homeing = homeing * 0.999 + res*0.001;
	return res;
}

double RejectGaussian(){
	static unsigned char sign=0;
	double u1,u2,z;
	// Attempt to find a random variable of a Gaussian fit, quickly...
	// By pure rejection...
	while (1){
		u1=Random01(); // pick a number
		u2=Random01(); // any number...

		if (exp( -0.5*64.*u1*u1 )>u2) {
			z =  u1*8;
			break;
		}
	}

	return z; // unsigned return value...
	
	sign >>=1;
	if (sign<2) sign=RandomByte();
	return (sign&1) ? -z:z;
}

double BGaussian(){
	// Use the "Box" method to generate gaussian distribution
	double z;
	z=sqrt( -2. * log( 1.-Random01() ) )*sin( 2*M_PI*Random01() );	
if (z<0) return -z;
	return z;
}

double NBGaussian(){
	// Use the No-Sine "Box" method to generate gaussian distribution
	// Under calculates for some reason...
	// Add a scalar to scale the population standard deviation to 1.0
	// as a drift correction.
	static double lzz=1.0;
	double xs1, xs2, z, w, drift=1.0;
	do {
		xs1=2.0*Random01()-1; // Cancel any net mean value...
		xs2=Random01();
		w = xs1*xs1 + xs2*xs2;
	} while (w >= 1.0 );

	z = xs1 * sqrt( (-2.0 * log( w ) )/w ) * drift;
	w = z*z;
	drift += (w + lzz)>2.0 ? -1e12:1e12; // Too big scales down, ...
	lzz = w;
	if (z<0) return -z;
	return z;
}

#define Gaussian  BGaussian
int main(void) {
	double a=0.;
	double tmp;
	double i=0;
	int j=0;
	double frac=0;
	long long even=0;
	char display[48]="1.", next[48];
	
	randomizer=fopen("/dev/urandom","rb");
	InitRandom();

	while (1) {
		tmp = Gaussian();
		a+=tmp;
		i+=1.0;
		if (tmp>0) frac+=1.0;

		if ( (++j & 0xfffff)==0 ){
			double res = sqrt(a / i);
			char *p=display, *q=next;
		       	sprintf( next,"%f", res );

			while (*p && *q && *p==*q) {++p;++q;}
			if (!(++even&0x7)) { *p++=*q; } // occasionally...
			*p=0;
			printf("%s    %8.6f  %4.2e  \r",display,frac/i,i);
			fflush(stdout);
		}
	}
	fclose(randomizer);
}