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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
/*
* 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 <float.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>
#include "random.h"
#include "std.h"
/*
* NOTE: I choose a coordinate system in units of d so that the distance
* between the straight lines is 1.
*/
#define MIN_COORD -10.0
#define MAX_COORD 10.0
#define TOTAL 10000
// Length of the needle with respect to the distance between the lines, this
// is totally arbitrary.
#define L 0.7
double
random_angle(Random rnd[static 1])
{
/* Sampling of all of the directions in the unit circle by a simple,
* yet inefficient, accept-reject algorithm. Using trigonometric
* functions is allowed because they are entire functions i.e. they
* are equal to their Taylor expansion around the origin, which does
* not require pi to be calculated.
*/
double x, y;
do {
x = random_rannyu_range(rnd, -1.0, 1.0);
y = random_rannyu_range(rnd, -1.0, 1.0);
} while (pow(x, 2) + pow(y, 2) >= 1);
/* FIXME(mario): Don't use this, just return x, y */
return atan2(y, x);
}
bool
needle_hits_line(double x, double y, double alpha)
{
/* The needle 'hits' a line if the top y value is one unit bigger than
* the bottom y value. See the picture below. */
double yy = sin(alpha) * (double)(L) / 2.0;
double ytop = y + yy;
double ybottom = y - yy;
// clang-format off
/* y = 1 ---------------------------------------------
*
* / y = 0.49 THE NEEDLE 'HITS'
* y = 0 -------/-------------------------------------
* / y = -0.49
*
* y = -1 ---------------------------------------------
*
*
* y = 1 ---------------------------------------------
* / y = 0.7
* / y = 0.35 THE NEEDLE DOES NOT 'HIT'
* y = 0 ---------------------------------------------
*
*
* y = -1 ---------------------------------------------
*/
// clang-format on
return (int)fabs(floor(ytop) - floor(ybottom)) > DBL_EPSILON;
}
static int nblocks[] = { 5, 10, 20, 25, 50, 100, 200 };
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 (size_t i = 0; i < LEN(nblocks); i++) {
int nblock = nblocks[i];
int block_steps = TOTAL / nblock;
double sum = 0.0, sum2 = 0.0;
for (int block = 0; block < nblock; block++) {
int hits = 0;
double pi = 0.0;
for (int step = 0; step < block_steps; step++) {
double x = random_rannyu_range(
&rnd, MIN_COORD, MAX_COORD);
double y = random_rannyu_range(
&rnd, MIN_COORD, MAX_COORD);
if (needle_hits_line(
x, y, random_angle(&rnd))) {
hits++;
}
}
pi = 2 * L * (double)(block_steps) / (double)(hits);
sum += pi;
sum2 += pi * pi;
}
sum /= (double)(nblock);
printf("%d %lf %lf\n", nblock, sum,
sqrt(sum2 / (double)(nblock)-sum * sum));
}
random_save_seed(&rnd, "seed.out");
return 0;
}
|
