#include<iostream>#include<fstream>#include <cmath>#include<cstdlib>#include<time.h>using namespace std;const double r=3.442619855899;  //N = 128//const double r=3.21365765342241;//N = 64;//const double r =2.96130012126409; //N =32;const double rr = 1/r;   const double m1 = 2147483648.0, m2 =4294967296.;const double mm2 = 1/m2;const double v=9.91256303526217e-3;  //N =128//const double v =  0.02002445560588; //N = 64;//const double v =  0.04075874443221; //N =32const int N =128;inline unsigned long IUniform();inline double Uniform();inline double Gaussian();double Tail();void generate_table();void fsolve(int i,double& s2,double &k0);static long hh;    static unsigned long ii, k[N];  //ii is the last 8 bits of hh static double w[N],f[N],xx[N],d[N],x0[N],f0[N];int main(){	ofstream out("gau2.txt");	generate_table();    int Cout  =1000000;	int t1 = clock();	while(Cout--)	{		out<<Gaussian()<<endl;		//Gaussian();	}	out.close();	cout<<clock()-t1<<endl;	return 0;}inline unsigned long IUniform()   //該算法來源于網絡{	static unsigned long s2,ss2 = 86947731;	s2=ss2; ss2^=(ss2<<13); ss2^=(ss2>>17); ss2^=(ss2<<5);	return s2+ss2;}inline double Uniform(){	return (0.5 + (signed)IUniform()*mm2);}inline double Gaussian(){	hh=(signed)IUniform();	ii=hh&(N-1);	return abs(hh)<=k[ii] ? hh*w[ii] : Tail();}double Tail(){	static double x, y,z,zz;	for(;;)	{		x=hh*w[ii];      // ii==0, tail		if(ii==0)		{ 			do{				x=-log(Uniform())*rr;				y=-log(Uniform());			}while(y+y<x*x);			return (hh>0)? r+x : -r-x;		}		// ii>0, wedges 		//if( f[ii]+Uniform()*(f[ii-1]-f[ii]) < exp(-0.5*x*x) )			//return x;				//new rejection methods		z = Uniform();		//z = fabs(z)*m2<hh? z:hh/m2;		//if(f[ii]+z*(f[ii-1]-f[ii]) < f0[ii]+d[ii]*(x-x0[ii]))		//	return x;		//else		if( f[ii]+z*(f[ii-1]-f[ii]) < exp(-0.5*x*x))			return x;				hh=IUniform();		ii=hh&(N-1);		if(abs(hh)<k[ii])			return (hh*w[ii]);	}}//void generate_table()//(unsigned long jsrseed)//{ //	double q;//	int i;//	q=v/exp(-0.5*r*r);//	k[0]=((r/q)*m1);   //tail//	k[1]=0;	//	w[0]=q/m1;//	w[127]=r/m1;//	f[0]=1.;//	f[127]=exp(-0.5*r*r);//	double xi = r,xj = r;//	for(i=126;i>=1;i--)//	{//		 xi=sqrt(-2.*log(v/xi+exp(-0.5*xi*xi)));//		 k[i+1]=((xi/xj)*m1);//		 xj = xi;//		 f[i]=exp(-0.5*xi*xi);//		 w[i]=xi/m1;//	}//}void generate_table(){ 	double q;	int i;	q=v/exp(-0.5*r*r);	k[0]=((r/q)*m1);   //tail	k[1]=0;		w[0]=q/m1;	w[127]=r/m1;	f[0]=1.;	f[N-1]=exp(-0.5*r*r);	xx[N-1] = r;xx[N-2] = r;xx[0]=0;	for(i=N-2;i>=1;i--)	{		 xx[i]=sqrt(-2.*log(v/xx[i+1]+exp(-0.5*xx[i+1]*xx[i+1])));		 k[i+1]=((xx[i]/xx[i+1])*m1);		 f[i]=exp(-0.5*xx[i]*xx[i]);		 w[i]=xx[i]/m1;	}	double a,b,so;	for(i=1;i<=N-1;i++)	{		if(xx[i]<=0.5)		{			//a = f[i-1]-f[i];			//b = xx[i]*f[i]*(xx[i]-xx[i-1]);			d[i] = (f[i]-f[i-1])/(xx[i]-xx[i-1]);			x0[i]=xx[i];			f0[i]=exp(-0.5*x0[i]*x0[i]);		}		else if(xx[i-1]>=0.5)		{			//b = f[i-1]-f[i];			//so = fsolve(i);			//a =exp(-so*so*0.5)-(f[i]-f[i-1])/(xx[i]-xx[i-1])*(so-xx[i-1])-f[i];			fsolve(i,x0[i],d[i]);			f0[i]=exp(-0.5*x0[i]*x0[i]);		}		else		{			d[i] = 0.0; //不優化			f0[i] =-2.0;		}		d[0]= 0.0;f0[0] =1.0;	}}void fsolve(int i,double& s2,double& k0){	k0 = (f[i]-f[i-1])/(xx[i]-xx[i-1]);	double s1 = (xx[i]+xx[i-1])/2.0;	s2= 0;	double tol = 1e-10;	while(fabs(s2-s1)>=tol)    //Newton	{		s2 = s1-(s1+k0*exp(s1*s1/2.0))/(1.0-s1*s1);		s1 = s2;	}}