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
|
// Copyright 2017 The Abseil Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// https://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// Generates gaussian_distribution.cc
//
// $ blaze run :gaussian_distribution_gentables > gaussian_distribution.cc
//
#include "absl/random/gaussian_distribution.h"
#include <cmath>
#include <cstddef>
#include <iostream>
#include <limits>
#include <string>
#include "absl/base/macros.h"
namespace absl {
namespace random_internal {
namespace {
template <typename T, size_t N>
void FormatArrayContents(std::ostream* os, T (&data)[N]) {
if (!std::numeric_limits<T>::is_exact) {
// Note: T is either an integer or a float.
// float requires higher precision to ensure that values are
// reproduced exactly.
// Trivia: C99 has hexadecimal floating point literals, but C++11 does not.
// Using them would remove all concern of precision loss.
os->precision(std::numeric_limits<T>::max_digits10 + 2);
}
*os << " {";
std::string separator = "";
for (size_t i = 0; i < N; ++i) {
*os << separator << data[i];
if ((i + 1) % 3 != 0) {
separator = ", ";
} else {
separator = ",\n ";
}
}
*os << "}";
}
} // namespace
class TableGenerator : public gaussian_distribution_base {
public:
TableGenerator();
void Print(std::ostream* os);
using gaussian_distribution_base::kMask;
using gaussian_distribution_base::kR;
using gaussian_distribution_base::kV;
private:
Tables tables_;
};
// Ziggurat gaussian initialization. For an explanation of the algorithm, see
// the Marsaglia paper, "The Ziggurat Method for Generating Random Variables".
// http://www.jstatsoft.org/v05/i08/
//
// Further details are available in the Doornik paper
// https://www.doornik.com/research/ziggurat.pdf
//
TableGenerator::TableGenerator() {
// The constants here should match the values in gaussian_distribution.h
static constexpr int kC = kMask + 1;
static_assert((ABSL_ARRAYSIZE(tables_.x) == kC + 1),
"xArray must be length kMask + 2");
static_assert((ABSL_ARRAYSIZE(tables_.x) == ABSL_ARRAYSIZE(tables_.f)),
"fx and x arrays must be identical length");
auto f = [](double x) { return std::exp(-0.5 * x * x); };
auto f_inv = [](double x) { return std::sqrt(-2.0 * std::log(x)); };
tables_.x[0] = kV / f(kR);
tables_.f[0] = f(tables_.x[0]);
tables_.x[1] = kR;
tables_.f[1] = f(tables_.x[1]);
tables_.x[kC] = 0.0;
tables_.f[kC] = f(tables_.x[kC]); // 1.0
for (int i = 2; i < kC; i++) {
double v = (kV / tables_.x[i - 1]) + tables_.f[i - 1];
tables_.x[i] = f_inv(v);
tables_.f[i] = v;
}
}
void TableGenerator::Print(std::ostream* os) {
*os << "// BEGIN GENERATED CODE; DO NOT EDIT\n"
"// clang-format off\n"
"\n"
"#include \"absl/random/gaussian_distribution.h\"\n"
"\n"
"namespace absl {\n"
"namespace random_internal {\n"
"\n"
"const gaussian_distribution_base::Tables\n"
" gaussian_distribution_base::zg_ = {\n";
FormatArrayContents(os, tables_.x);
*os << ",\n";
FormatArrayContents(os, tables_.f);
*os << "};\n"
"\n"
"} // namespace random_internal\n"
"} // namespace absl\n"
"\n"
"// clang-format on\n"
"// END GENERATED CODE";
*os << std::endl;
}
} // namespace random_internal
} // namespace absl
int main(int, char**) {
std::cerr << "\nCopy the output to gaussian_distribution.cc" << std::endl;
absl::random_internal::TableGenerator generator;
generator.Print(&std::cout);
return 0;
}
|
