zig.cpp

高斯随机数发生器

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <ctime>
#define NN 1024
#define KK 1000000

static unsigned long iz,jz,jsr=(((unsigned)time(NULL))%100000)*(((unsigned)time(NULL))%1000);
#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define UNI (.5 + (signed) SHR3 * .2328306e-9)

#define f(x) exp(-0.5*x*x)

using namespace std;
double R=4.03884984610951;
double S=0.001226324646353;
double tableZ[NN+1][2]={0.0};

void formthetable(double a[][2])
{
	int i;
	a[NN-1][0]=R;
	a[NN-1][1]=f(a[NN-1][0]);
	a[NN][0]=S/a[NN-1][1];
	a[NN][1]=0;
	for(i=NN-1;i>1;i--)
	{
		a[i-1][1]=a[i][1]+S/a[i][0];
		a[i-1][0]=sqrt(-2*log(a[i-1][1]));
	}
}

int temp;
double XY[2],x,y;
double Gauss()
{	
LOOP:
	temp=rand()&(NN-1);
	XY[0]=(signed(SHR3)/2147483648.0)*tableZ[temp+1][0];
	if(XY[0]<tableZ[temp][0]&&XY[0]>-tableZ[temp][0])
		return(XY[0]);
	else
	{
		if(temp==1023)
		{
			do
			{
				x=-log(UNI)/R;
				y=-log(UNI);
			}
			while(y+y<x*x);
			return ((XY[0]>0)?R+x:-R-x);
		}
		else
		{
			XY[1]=tableZ[temp+1][1]+UNI*(tableZ[temp][1]-tableZ[temp+1][1]);
			if(XY[1]<f(XY[0]))
				return(XY[0]);
			else goto LOOP;
		}	
	}
}

void main()
{
	int i,j;
	double result,max=-1.0,a[1000]={0},begin;
	double mean=0.0,square=0.0;
	ofstream output;
	output.open("result1.txt");
	formthetable(tableZ);
	srand((((unsigned)time(NULL))%100)*(((unsigned)time(NULL))%10000));
	//for(i=0;i<NN;i++)
		//cout<<setprecision(20)<<tableZ[i][1]<<endl;
	//begin=clock();
	for(i=0;i<KK;i++)
	{
		//result=Gauss();
		output<<setprecision(20)<<Gauss()<<endl;
//		j=long((result+8)*1000.0/16.0);
//		a[j]+=1.0;
		//max=max>fabs(result)?max:fabs(result);
//		mean+=result;
//		square+=result*result;
	}
//	for(i=0; i<1000; i++)
//	{
//		a[i]/=160000.0;
//		output<<setprecision(20)<<a[i]<<endl;
//	}
//	mean/=KK;
//	square/=KK;
//	square-=mean*mean;
//	cout<<mean<<endl<<square<<endl;
//	cout<<max<<endl;
	//cout<<clock()-begin<<endl;
}