RandomThread.cpp
5.15 KB
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
/**
* \file RandomThread.cpp
* \brief Example of parallelization with %RandomLib using threads
*
* Compile/link with, e.g.,\n
* g++ -I../include -O2 -funroll-loops -fopenmp
* -o RandomThread RandomThread.cpp ../src/Random.cpp\n
* ./RandomThread
*
* See \ref parallel, for a description of this example.
*
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#include <iostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <string>
#include <sstream>
#include <numeric>
#include <RandomLib/Random.hpp>
#if HAVE_OPENMP
#include <omp.h>
#endif
using RandomLib::Random;
double dotrials(Random& r, int d, long long n) {
// Require d > 0, n > 0
double w = 0;
for ( ; n--; ) { // Iterate n times
double h = 0;
for (int i = 1; i < d; ++i) { // Iterate d-1 times
double x = 2 * r.FixedS(); // x is in (-1, 1)
h += x * x; // cumulative radius^2
if (h >= 1) break; // Point can't be in sphere; bail out,
}
// If h < 1 then inside a (d-1) dimensional unit sphere at radius
// sqrt(h), so extent of last dimension is +/- sqrt(1 -h)
w += h < 1 ? sqrt(1 - h) : 0;
}
return w;
}
double result(int d, long long n, double w) {
// Volume of (d-1) dimensional box = 2^(d-1).
// Multiply by another 2 to account for +/- extent in last dimension.
return double(1U << d) * w / double(n);
}
int usage(const std::string& name, int retval) {
( retval == 0 ? std::cout : std::cerr )
<< "Usage: \n" << name
<< " [-d dim] [-n nsamp] [-k ntask] [-l stride] [-s seed] [-h]\n"
<< "Estimate volume of n-dimensional unit sphere\n";
return retval;
}
int main(int argc, char* argv[]) {
int d = 4; // Number of dimensions
long long n = 10000000; // Number of trials 10^7
int k = 100; // Number of tasks
int l = 4; // The leapfrogging stride
std::string seedstr;
bool seedgiven = false;
for (int m = 1; m < argc; ++m) {
std::string arg(argv[m]);
if (arg == "-d") {
if (++m == argc) return usage(argv[0], true);
std::istringstream str(argv[m]);
char c;
if (!(str >> d) || (str >> c) || d <= 0) {
std::cerr << "Number of dimensions " << argv[m]
<< " is not a positive number\n";
return 1;
}
} else if (arg == "-n") {
if (++m == argc) return usage(argv[0], true);
std::istringstream str(argv[m]);
char c;
double fn;
if (!(str >> fn) || (str >> c) || fn < 1 ||
fn > double(std::numeric_limits<long long>::max())) {
std::cerr << "Total count " << argv[m] << " is not a positive number\n";
return 1;
}
n = (long long)(fn);
} else if (arg == "-k") {
if (++m == argc) return usage(argv[0], true);
std::istringstream str(argv[m]);
char c;
if (!(str >> k) || (str >> c) || k <= 0) {
std::cerr << "Number of tasks " << argv[m]
<< " is not a positive number\n";
return 1;
}
} else if (arg == "-l") {
if (++m == argc) return usage(argv[0], true);
std::istringstream str(argv[m]);
char c;
if (!(str >> l) || (str >> c) || l <= 0) {
std::cerr << "Leapfrog stride " << argv[m]
<< " is not a positive number\n";
return 1;
}
} else if (arg == "-s") {
seedgiven = true;
if (++m == argc) return usage(argv[0], 1);
seedstr = std::string(argv[m]);
} else
return usage(argv[0], arg != "-h");
}
std::vector<unsigned long> master_seed =
seedgiven ? Random::StringToVector(seedstr) : Random::SeedVector();
std::cout << "Estimate volume of a " << d
<< "-dimensional sphere;\nsamples = -n " << n
<< "; tasks = -k " << k
<< "; leapfrog stride = -l " << l << ";\nusing " << Random::Name()
<< "\nwith master seed = -s "
<< Random::VectorToString(master_seed) << ".\n";
std::cout << "Estimated volume = "
<< std::fixed << std::setprecision(8) << std::flush;
// Reserve room for a task id at the end of the seed
master_seed.push_back(0);
// Fill weight vector with NaNs to verify that all tasks have completed
std::vector<double> w(k, std::numeric_limits<double>::quiet_NaN());
#if HAVE_OPENMP
#pragma omp parallel for
#endif
for (int i = 0; i < k; ++i) { // the main loop over tasks
Random r; // task specific Random
{
std::vector<unsigned long> seed(master_seed); // task specific seed
seed.back() = i / l; // include task id in seed
r.Reseed(seed);
// Turn on leapfrogging with an offset that depends on the task id
r.SetStride(l, i % l);
}
// Do the work; last argument splits n exactly into k pieces
w[i] = dotrials(r, d, (n * (i + 1))/k - (n * i)/k);
}
// Sum up the weights from the individual tasks
double weight = accumulate(w.begin(), w.end(), 0.0);
// Compute the result
std::cout << result(d, n, weight) << "\n";
return 0;
}