gaussian1.cpp

用ziggurat 算法实现高斯随机数的生成。效率很好。

#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;
	}
}