#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 =32
const 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;
}
}