### 隨機數生成器

Park和MIller提供的“最小標准”隨機數生成器

```#define a 16807
#define m 2147483647    // 2^31 -1
#define q (m / a)
#define r (m % a)

static long int seed = 1;

/*
X <-- (aX + c) mod m;  (c = 0)
return [1, 2147483647]隨機數
*/
long int PMrand()
{
long int hi = seed / q;
long int lo = seed % q;
long int test = a * lo - r * hi;
if(test > 0)    seed = test;
else     seed = test + m;

return seed;
}

/*
X <-- (aX + c) mod m;  (c = 0)
return [0, 1]隨機數
*/
double PMrandDouble()
{
long int hi = seed / q;
long int lo = seed % q;
long int test = a * lo - r * hi;
if(test > 0)    seed = test;
else     seed = test + m;

return (double)seed / m;
}
```

```#include <stdlib.h>
#include <math.h>
#define NSUM 25

double GaussRand()
{
double x = 0;
for(int i = 0; i < NSUM; i ++){
x += (double)rand() / RAND_MAX;
}
x -= NSUM / 2.0;
x /= sqrt(NSUM / 12.0);

return x;
}```

```#include <stdlib.h>
#include <math.h>
#define PI 3.141592654

double GaussRand()
{
static double U, V;
static int phase = 0;
double Z;
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;
}```

```#include <stdlib.h>
#include <math.h>

double GaussRand()
{
static doubel V1, V2, S;
static int phase = 0;
double X;
if(phase == 0){
do{
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;

V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
} while(S >= 1 || S == 0);

X = V1 * sqrt(-2 * log(S) / S);
}
else{
X = V2 * sqrt(-2 * log(S) / S);
}
phase = 1 - phase;

return X;
}

```