方法一:
使用由Box和Muller提供的,在Knuth的网上讨论过的方法:
#include <stdio.h>
#include <math.h>
#define PI 3.1415926
double gaussrand()
{
static double U,V;
static int phase = 0;
double e;
if(phase == 0)
{
U = (rand()+1.0)/(RAND_MAX+2.0);
V = rand()/(RAND_MAX+1.0);
z = sqrt(-2 * log(U)) * sin(2* PI *V);
}else
{
z = sqrt(-2 * log(U)) * cos(2* PI *V);
}
phase = 1 - phase;
return z;
}
方法二:
使用最初由Marsaglia提供的方法:
#include <stdio.h>
#include <math.h>
double gaussrand()
{
static double V1,V2,S = 0;
static int phase = 0;
double X;
if(phase == 0)
{
do{
U1 = (rand()+0.0)/(RAND_MAX+0.0);
U2 = (rand()+0.0)/(RAND_MAX+0.0);
V1 = 2* U1 - 1;
V2 = 2* U2 - 1;
S = V1*V1 + V2*V2;
}while( S >= 1 || S == 0);
X = sqrt( -2 * log(S) / S ) * V1 ;
}else
{
X = sqrt( -2 * log(S) / S ) * V2 ;
}
phase = 1 - phase;
return X;
}