#include <stdio.h>
#include <stdlib.h>
#include <math.h>
static int ix = 1, iy = 1, iz = 1; /* 1--30000の任意の整数 */
/*
* 改良版乱数(Wichmann--Hillの乱数発生法)
*/
void init_rnd(int x, int y, int z) /* x, y, z = 1..30000 */
{
ix = x;
iy = y;
iz = z;
}
double Rand(void) /* 0 <= rnd() < 1 */
{
double r;
ix = 171 * (ix % 177) - 2 * (ix / 177);
iy = 172 * (iy % 176) - 35 * (iy / 176);
iz = 170 * (iz % 178) - 63 * (iz / 178);
if (ix < 0) ix += 30269;
if (iy < 0) iy += 30307;
if (iz < 0) iz += 30323;
r = ix / 30269.0 + iy / 30307.0 + iz / 30323.0;
while (r >= 1.0)
r = r - 1.0;
return r;
}
/*
* ボックス=ミューラー法
*/
double NormalRandom(double mu, double sigma)
{
double t, u, r;
static int sw = 0;
static double rs;
if (sw == 0) {
t
= sqrt(-2.0 * log(1.0 - Rand
())); u = 2.0 * M_PI * Rand();
sw = 1;
} else {
r = rs;
sw = 0;
}
return sigma * r + mu;
}
int main(void)
{
int i, n, ix, *histo;
double mu, sigma, x, s1 = 0.0, s2 = 0.0;
double *array;
init_rnd
((time(NULL
) % 30000) + 1, ((time(NULL
) + 23456) % 30000) + 1, ((time(NULL
) + 45678) % 30000) + 1); /* 1..30000 の任意の整数で初期化. */// printf("個数 n = ");
// scanf("%d", &n);
n = 10000;
// printf("平均 μ = ");
// scanf("%lf", &mu);
mu = 1.5;
// printf("分散 σ = ");
// scanf("%lf", &sigma);
sigma = 0.5;
if ((histo
= (int *)malloc(sizeof(int) * n
)) == NULL
) if ((array
= (double *)malloc(sizeof(double) * n
)) == NULL
)
/*
* 乱数はarray[]に保存する
* 度数分布を同時に求める
*/
for (i = 0; i < n; i++) {
x = array[i] = NormalRandom(mu, sigma);
ix
= (int)floor(2.0 * x
) + 10; if (ix >= 0 && ix < 20)
histo[ix]++;
s1 += x;
s2 += x * x;
}
/*
* 平均と標準偏差の計算
*/
for (i = 0; i < 20; i++)
printf("%4.1f -- %4.1f: %5.1f%%\n", 0.5 * (i
- 10), 0.5 * (i
- 9), 100.0 * histo
[i
] / n
); s1 /= n;
s2
= sqrt((s2
- n
* s1
* s1
) / (n
- 1)); printf("平均 %g 標準偏差 %g\n", s1
, s2
);
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCnN0YXRpYyBpbnQgaXggPSAxLCBpeSA9IDEsIGl6ID0gMTsgIC8qIDEtLTMwMDAw44Gu5Lu75oSP44Gu5pW05pWwICovCgovKgogKiDmlLnoia/niYjkubHmlbAoV2ljaG1hbm4tLUhpbGzjga7kubHmlbDnmbrnlJ/ms5UpCiAqLwp2b2lkIGluaXRfcm5kKGludCB4LCBpbnQgeSwgaW50IHopICAvKiB4LCB5LCB6ID0gMS4uMzAwMDAgKi8KewogIGl4ID0geDsKICBpeSA9IHk7CiAgaXogPSB6Owp9Cgpkb3VibGUgUmFuZCh2b2lkKSAgLyogMCA8PSBybmQoKSA8IDEgKi8KewogIGRvdWJsZSByOwoKICBpeCA9IDE3MSAqIChpeCAlIDE3NykgLSAgMiAqIChpeCAvIDE3Nyk7CiAgaXkgPSAxNzIgKiAoaXkgJSAxNzYpIC0gMzUgKiAoaXkgLyAxNzYpOwogIGl6ID0gMTcwICogKGl6ICUgMTc4KSAtIDYzICogKGl6IC8gMTc4KTsKICBpZiAoaXggPCAwKSBpeCArPSAzMDI2OTsKICBpZiAoaXkgPCAwKSBpeSArPSAzMDMwNzsKICBpZiAoaXogPCAwKSBpeiArPSAzMDMyMzsKICByID0gaXggLyAzMDI2OS4wICsgaXkgLyAzMDMwNy4wICsgaXogLyAzMDMyMy4wOwogIHdoaWxlIChyID49IDEuMCkKICAgIHIgPSByIC0gMS4wOwogIAogIHJldHVybiByOwp9CgovKgogKiDjg5zjg4Pjgq/jgrnvvJ3jg5/jg6Xjg7zjg6njg7zms5UKICovCmRvdWJsZSBOb3JtYWxSYW5kb20oZG91YmxlIG11LCBkb3VibGUgc2lnbWEpCnsKICBkb3VibGUgdCwgdSwgcjsKICBzdGF0aWMgaW50IHN3ID0gMDsKICBzdGF0aWMgZG91YmxlIHJzOwogIAogIGlmIChzdyA9PSAwKSB7CiAgICB0ID0gc3FydCgtMi4wICogbG9nKDEuMCAtIFJhbmQoKSkpOwogICAgdSA9IDIuMCAqIE1fUEkgKiBSYW5kKCk7CiAgICByID0gdCAqIGNvcyh1KTsKICAgIHJzID0gdCAqIHNpbih1KTsKICAgIHN3ID0gMTsKICB9IGVsc2UgewogICAgciA9IHJzOwogICAgc3cgPSAwOwogIH0KICAKICByZXR1cm4gc2lnbWEgKiByICsgbXU7Cn0KCmludCBtYWluKHZvaWQpCnsKICBpbnQgaSwgbiwgaXgsICpoaXN0bzsKICBkb3VibGUgbXUsIHNpZ21hLCB4LCBzMSA9IDAuMCwgczIgPSAwLjA7CiAgZG91YmxlICphcnJheTsKICAKICBpbml0X3JuZCgodGltZShOVUxMKSAlIDMwMDAwKSArIDEsICgodGltZShOVUxMKSArIDIzNDU2KSAlIDMwMDAwKSArIDEsICgodGltZShOVUxMKSArIDQ1Njc4KSAlIDMwMDAwKSArIDEpOyAgLyogMS4uMzAwMDAg44Gu5Lu75oSP44Gu5pW05pWw44Gn5Yid5pyf5YyWLiAqLwovLyAgcHJpbnRmKCLlgIvmlbAgbiA9ICIpOwovLyAgc2NhbmYoIiVkIiwgJm4pOwogIG4gPSAxMDAwMDsKLy8gIHByaW50Zigi5bmz5Z2HIM68ID0gIik7Ci8vICBzY2FuZigiJWxmIiwgJm11KTsKICBtdSA9IDEuNTsKLy8gIHByaW50Zigi5YiG5pWjIM+DID0gIik7Ci8vICBzY2FuZigiJWxmIiwgJnNpZ21hKTsKICBzaWdtYSA9IDAuNTsKICAKICBpZiAoKGhpc3RvID0gKGludCAqKW1hbGxvYyhzaXplb2YoaW50KSAqIG4pKSA9PSBOVUxMKQogICAgZXhpdCgxKTsKICBpZiAoKGFycmF5ID0gKGRvdWJsZSAqKW1hbGxvYyhzaXplb2YoZG91YmxlKSAqIG4pKSA9PSBOVUxMKQogICAgZXhpdCgxKTsKICAKICAvKgogICAqIOS5seaVsOOBr2FycmF5W13jgavkv53lrZjjgZnjgosKICAgKiDluqbmlbDliIbluIPjgpLlkIzmmYLjgavmsYLjgoHjgosKICAgKi8KICBmb3IgKGkgPSAwOyBpIDwgbjsgaSsrKSB7CiAgICB4ID0gYXJyYXlbaV0gPSBOb3JtYWxSYW5kb20obXUsIHNpZ21hKTsKICAgIGl4ID0gKGludClmbG9vcigyLjAgKiB4KSArIDEwOwogICAgaWYgKGl4ID49IDAgJiYgaXggPCAyMCkKICAgICAgaGlzdG9baXhdKys7CiAgICBzMSArPSB4OwogICAgczIgKz0geCAqIHg7CiAgfQoKICAvKgogICAqIOW5s+Wdh+OBqOaomea6luWBj+W3ruOBruioiOeulwogICAqLwogIGZvciAoaSA9IDA7IGkgPCAyMDsgaSsrKQogICAgcHJpbnRmKCIlNC4xZiAtLSAlNC4xZjogJTUuMWYlJVxuIiwgMC41ICogKGkgLSAxMCksIDAuNSAqIChpIC0gOSksIDEwMC4wICogaGlzdG9baV0gLyBuKTsKICBzMSAvPSBuOwogIHMyID0gc3FydCgoczIgLSBuICogczEgKiBzMSkgLyAobiAtIDEpKTsKICBwcmludGYoIuW5s+WdhyAlZyAg5qiZ5rqW5YGP5beuICVnXG4iLCBzMSwgczIpOwogIAogIGZyZWUoYXJyYXkpOwogIGZyZWUoaGlzdG8pOwogIAogIHJldHVybiAwOwp9Cg==