RandomExact.cpp
6.57 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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
/**
* \file RandomExact.cpp
* \brief Using %RandomLib to generate exact random results
*
* Compile/link with, e.g.,\n
* g++ -I../include -O2 -funroll-loops
* -o RandomExact RandomExact.cpp ../src/Random.cpp\n
* ./RandomExact
*
* See \ref otherdist, for more information.
*
* Copyright (c) Charles Karney (2006-2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#include <iostream>
#include <utility>
#include <RandomLib/Random.hpp>
#include <RandomLib/RandomNumber.hpp>
#include <RandomLib/ExactNormal.hpp>
#include <RandomLib/ExactExponential.hpp>
#include <RandomLib/UniformInteger.hpp>
#include <RandomLib/DiscreteNormal.hpp>
#include <RandomLib/DiscreteNormalAlt.hpp>
#include <RandomLib/InversePiProb.hpp>
#include <RandomLib/InverseEProb.hpp>
int main(int, char**) {
// Create r with a random seed
RandomLib::Random r; r.Reseed();
std::cout << "Using " << r.Name() << "\n"
<< "with seed " << r.SeedString() << "\n\n";
{
std::cout
<< "Sampling exactly from the normal distribution. First number is\n"
<< "in binary with ... indicating an infinite sequence of random\n"
<< "bits. Second number gives the corresponding interval. Third\n"
<< "number is the result of filling in the missing bits and rounding\n"
<< "exactly to the nearest representable double.\n";
const int bits = 1;
RandomLib::ExactNormal<bits> ndist;
long long num = 20000000ll;
long long bitcount = 0;
int numprint = 16;
for (long long i = 0; i < num; ++i) {
long long k = r.Count();
RandomLib::RandomNumber<bits> x = ndist(r); // Sample
bitcount += r.Count() - k;
if (i < numprint) {
std::pair<double, double> z = x.Range();
std::cout << x << " = " // Print in binary with ellipsis
<< "(" << z.first << "," << z.second << ")"; // Print range
double v = x.Value<double>(r); // Round exactly to nearest double
std::cout << " = " << v << "\n";
} else if (i == numprint)
std::cout << std::flush;
}
std::cout
<< "Number of bits needed to obtain the binary representation averaged\n"
<< "over " << num << " samples = " << bitcount/double(num) << "\n\n";
}
{
std::cout
<< "Sampling exactly from exp(-x). First number is in binary with\n"
<< "... indicating an infinite sequence of random bits. Second\n"
<< "number gives the corresponding interval. Third number is the\n"
<< "result of filling in the missing bits and rounding exactly to\n"
<< "the nearest representable double.\n";
const int bits = 1;
RandomLib::ExactExponential<bits> edist;
long long num = 50000000ll;
long long bitcount = 0;
int numprint = 16;
for (long long i = 0; i < num; ++i) {
long long k = r.Count();
RandomLib::RandomNumber<bits> x = edist(r); // Sample
bitcount += r.Count() - k;
if (i < numprint) {
std::pair<double, double> z = x.Range();
std::cout << x << " = " // Print in binary with ellipsis
<< "(" << z.first << "," << z.second << ")"; // Print range
double v = x.Value<double>(r); // Round exactly to nearest double
std::cout << " = " << v << "\n";
} else if (i == numprint)
std::cout << std::flush;
}
std::cout
<< "Number of bits needed to obtain the binary representation averaged\n"
<< "over " << num << " samples = " << bitcount/double(num) << "\n\n";
}
{
std::cout
<< "Sampling exactly from the discrete normal distribution with\n"
<< "sigma = 7 and mu = 1/2.\n";
RandomLib::DiscreteNormal<int> gdist(7,1,1,2);
long long num = 50000000ll;
long long count = r.Count();
int numprint = 16;
for (long long i = 0; i < num; ++i) {
int k = gdist(r); // Sample
if (i < numprint)
std::cout << k << " ";
else if (i == numprint)
std::cout << std::endl;
}
count = r.Count() - count;
std::cout
<< "Number of random variates needed averaged\n"
<< "over " << num << " samples = " << count/double(num) << "\n\n";
}
{
std::cout
<< "Sampling exactly from the discrete normal distribution with\n"
<< "sigma = 1024 and mu = 1/7. First result printed is a uniform\n"
<< "range (with covers a power of two). The second number is the\n"
<< "result of sampling additional bits within that range to obtain\n"
<< "a definite result.\n";
RandomLib::DiscreteNormalAlt<int,1> gdist(1024,1,1,7);
long long num = 20000000ll;
long long count = r.Count();
long long entropy = 0;
int numprint = 16;
for (long long i = 0; i < num; ++i) {
RandomLib::UniformInteger<int,1> u = gdist(r);
entropy += u.Entropy();
if (i < numprint)
std::cout << u << " = ";
int k = u(r);
if (i < numprint)
std::cout << k << "\n";
else if (i == numprint)
std::cout << std::flush;
}
count = r.Count() - count;
std::cout
<< "Number of random bits needed for full result (for range) averaged\n"
<< "over " << num << " samples = " << count/double(num) << " ("
<< (count - entropy)/double(num) << ")\n\n";
}
{
std::cout
<< "Random bits with 1 occurring with probability 1/pi exactly:\n";
long long num = 100000000ll;
int numprint = 72;
RandomLib::InversePiProb pp;
long long nbits = 0;
long long k = r.Count();
for (long long i = 0; i < num; ++i) {
bool b = pp(r);
nbits += int(b);
if (i < numprint) std::cout << int(b);
else if (i == numprint) std::cout << "..." << std::flush;
}
std::cout << "\n";
std::cout << "Frequency of 1 averaged over " << num << " samples = 1/"
<< double(num)/nbits << "\n"
<< "bits/sample = " << (r.Count() - k)/double(num) << "\n\n";
}
{
std::cout
<< "Random bits with 1 occurring with probability 1/e exactly:\n";
long long num = 200000000ll;
int numprint = 72;
RandomLib::InverseEProb ep;
long long nbits = 0;
long long k = r.Count();
for (long long i = 0; i < num; ++i) {
bool b = ep(r);
nbits += int(b);
if (i < numprint) std::cout << int(b);
else if (i == numprint) std::cout << "..." << std::flush;
}
std::cout << "\n";
std::cout << "Frequency of 1 averaged over " << num << " samples = 1/"
<< double(num)/nbits << "\n"
<< "bits/sample = " << (r.Count() - k)/double(num) << "\n";
}
return 0;
}