#!/usr/local/bin/tcc -run
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
#define NUMBER_OF_POINTS 1000
static double pi = 3.14159265358979323846;
/*
#######################################################################
# #
# Box-Muller Algorithm to generate gaussian random numbers N(0, 1) #
# #
#######################################################################
*/
/* --- 一様乱数を生成する関数のプロタイプ宣言 --- */
void iid(double x[]);
main()
{
int i;
int n;
double r1[NUMBER_OF_POINTS];
double r2[NUMBER_OF_POINTS];
double x1[NUMBER_OF_POINTS];
double x2[NUMBER_OF_POINTS];
FILE *fp;
time_t st;
time_t et;
/* --- 時刻を計測する --- */
time(&st);
printf("\n%s\n",ctime(&st));
/* --- 一様乱数を生成する --- */
iid(x1);
printf("\n If you want to go, press 0.\n");
scanf("%d", &n);
iid(x2);
/* --- 正規分布乱数を生成する --- */
for ( n = 0 ; n < NUMBER_OF_POINTS ; n++ ) {
r1[n] = sqrt(-2.0*log(x1[n])) * cos(2.0*pi*x2[n]);
r2[n] = sqrt(-2.0*log(x2[n])) * sin(2.0*pi*x1[n]);
}
/* --- 計算結果を保存する --- */
if ( ( fp = fopen("gauss_1.dat", "w" ) ) == 0 ) {
fprintf(stderr, "\n\tSorry! I can't open gauss_1.dat.\n");
}
for ( n = 0 ; n < NUMBER_OF_POINTS ; n++ )
fprintf(fp, "%f\n", r1[n] );
fclose(fp);
if ( ( fp = fopen("gauss_2.dat", "w" ) ) == 0 ) {
fprintf(stderr, "\n\tSorry! I can't open gauss_2.dat.\n");
}
for ( n = 0 ; n < NUMBER_OF_POINTS ; n++ )
fprintf(fp, "%f\n", r2[n] );
fclose(fp);
time(&et);
printf("\n%s\n",ctime(&st));
printf("%s\n",ctime(&et));
printf("\n\t Gaussian random numbers are stored in gauss_1.dat and gauss_2.dat.\n");
return 0;
}
/* --- 一様乱数を生成する関数 --- */
void iid(double x[])
{
int j;
int n;
long idum;
long k;
static long idum2 = 123456789;
static long iy = 0;
static long iv[NTAB];
double temp;
idum = - time(NULL);
for ( n = 0 ; n < NUMBER_OF_POINTS ; n++ ) {
if ( idum <= 0 ) {
if ( -idum < 1 )
idum = 1;
else
idum = -idum;
idum2 = idum;
for ( j = NTAB + 7 ; j >= 0 ; j-- ) {
k = idum / IQ1;
idum = IA1 * ( idum - k*IQ1) - k*IR1;
if ( idum < 0 )
idum += IM1;
if ( j < NTAB )
iv[j] = idum;
}
iy = iv[0];
}
k = idum / IQ1;
idum = IA1 * (idum - k*IQ1) - k*IR1;
if ( idum < 0 )
idum += IM1;
k = idum2 / IQ2;
idum2 = IA2 * ( idum2 - k*IQ2 ) - k*IR2;
if ( idum2 < 0 )
idum2 += IM2;
j = iy / NDIV;
iy = iv[j] - idum2;
iv[j] = idum;
if ( iy < 1 )
iy += IMM1;
if ( ( temp = AM * iy ) > RNMX )
x[n] = RNMX;
else
x[n] = temp;
}
}