1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
/*
* Copyright ©️ 2024 Mario Forzanini <mf@marioforzanini.com>
*
* This file is part of my exercises for LSN.
*
* My exercises for LSN are free software: you can redistribute them
* and/or modify them under the terms of the GNU General Public License
* as published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* My exercises for LSN are distributed in the hope that they will be
* useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
* See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with my exercises for LSN. If not, see
* <https://www.gnu.org/licenses/>.
*
*/
#include <errno.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>
#include "random.h"
#include "std.h"
// Total number of "throws"
#define TOTAL 10000
// Number of blocks
#define NBLOCKS 100
#define BLOCK_STEPS (TOTAL / NBLOCKS)
/*! Array in which results are stored for each block of the blocking
* method. */
static double blocks[NBLOCKS] = { 0 };
/*! Average to calculate the standard deviation of the mean. */
static double prog_avg[NBLOCKS] = { 0 };
/*! Squared average to calculate the standard deviation of the mean */
static double prog_avg2[NBLOCKS] = { 0 };
/*! Successive values of the variance as the number of blocks considered
* enlarges. */
static double variances[NBLOCKS] = { 0 };
int
main(void)
{
FILE *primes, *input;
Random rnd = { 0 };
if (!(primes = fopen("Primes", "r"))) {
STD_ERROR("Error opening Primes: %s\n", strerror(errno));
return 1;
} else if (!(input = fopen("seed.in", "r"))) {
STD_ERROR("Error opening seed.in: %s\n", strerror(errno));
return 1;
}
if (!random_init(&rnd, primes, input)) {
return 1;
}
for (int block = 0; block < NBLOCKS; block++) {
double sum = 0.0;
for (int i = 0; i < TOTAL / NBLOCKS; i++) {
sum += random_rannyu(&rnd);
}
blocks[block] = sum / (double)(BLOCK_STEPS);
}
for (int block = 0; block < NBLOCKS; block++) {
for (int i = 0; i < block + 1; i++) {
double Ai = blocks[i];
prog_avg[block] += Ai;
prog_avg2[block] += Ai * Ai;
}
prog_avg[block] /= block + 1;
prog_avg2[block] /= block + 1;
}
for (int block = 1; block < NBLOCKS; block++) {
variances[block] = sqrt(
(prog_avg2[block] - prog_avg[block] * prog_avg[block])
/ (double)block);
}
for (int block = 0; block < NBLOCKS; block++) {
printf("%lf %lf\n", prog_avg[block] - 0.5, variances[block]);
}
if (!random_save_seed(&rnd, "seed.out")) {
return 1;
}
return 0;
}
|
