1 // Copyright 2017 The Abseil Authors.
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // https://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 // Generates gaussian_distribution.cc
16 //
17 // $ blaze run :gaussian_distribution_gentables > gaussian_distribution.cc
18 //
19 #include "absl/random/gaussian_distribution.h"
20
21 #include <cmath>
22 #include <cstddef>
23 #include <iostream>
24 #include <limits>
25 #include <string>
26
27 #include "absl/base/macros.h"
28
29 namespace absl {
30 ABSL_NAMESPACE_BEGIN
31 namespace random_internal {
32 namespace {
33
34 template <typename T, size_t N>
FormatArrayContents(std::ostream * os,T (& data)[N])35 void FormatArrayContents(std::ostream* os, T (&data)[N]) {
36 if (!std::numeric_limits<T>::is_exact) {
37 // Note: T is either an integer or a float.
38 // float requires higher precision to ensure that values are
39 // reproduced exactly.
40 // Trivia: C99 has hexadecimal floating point literals, but C++11 does not.
41 // Using them would remove all concern of precision loss.
42 os->precision(std::numeric_limits<T>::max_digits10 + 2);
43 }
44 *os << " {";
45 std::string separator = "";
46 for (size_t i = 0; i < N; ++i) {
47 *os << separator << data[i];
48 if ((i + 1) % 3 != 0) {
49 separator = ", ";
50 } else {
51 separator = ",\n ";
52 }
53 }
54 *os << "}";
55 }
56
57 } // namespace
58
59 class TableGenerator : public gaussian_distribution_base {
60 public:
61 TableGenerator();
62 void Print(std::ostream* os);
63
64 using gaussian_distribution_base::kMask;
65 using gaussian_distribution_base::kR;
66 using gaussian_distribution_base::kV;
67
68 private:
69 Tables tables_;
70 };
71
72 // Ziggurat gaussian initialization. For an explanation of the algorithm, see
73 // the Marsaglia paper, "The Ziggurat Method for Generating Random Variables".
74 // http://www.jstatsoft.org/v05/i08/
75 //
76 // Further details are available in the Doornik paper
77 // https://www.doornik.com/research/ziggurat.pdf
78 //
TableGenerator()79 TableGenerator::TableGenerator() {
80 // The constants here should match the values in gaussian_distribution.h
81 static constexpr int kC = kMask + 1;
82
83 static_assert((ABSL_ARRAYSIZE(tables_.x) == kC + 1),
84 "xArray must be length kMask + 2");
85
86 static_assert((ABSL_ARRAYSIZE(tables_.x) == ABSL_ARRAYSIZE(tables_.f)),
87 "fx and x arrays must be identical length");
88
89 auto f = [](double x) { return std::exp(-0.5 * x * x); };
90 auto f_inv = [](double x) { return std::sqrt(-2.0 * std::log(x)); };
91
92 tables_.x[0] = kV / f(kR);
93 tables_.f[0] = f(tables_.x[0]);
94
95 tables_.x[1] = kR;
96 tables_.f[1] = f(tables_.x[1]);
97
98 tables_.x[kC] = 0.0;
99 tables_.f[kC] = f(tables_.x[kC]); // 1.0
100
101 for (int i = 2; i < kC; i++) {
102 double v = (kV / tables_.x[i - 1]) + tables_.f[i - 1];
103 tables_.x[i] = f_inv(v);
104 tables_.f[i] = v;
105 }
106 }
107
Print(std::ostream * os)108 void TableGenerator::Print(std::ostream* os) {
109 *os << "// BEGIN GENERATED CODE; DO NOT EDIT\n"
110 "// clang-format off\n"
111 "\n"
112 "#include \"absl/random/gaussian_distribution.h\"\n"
113 "\n"
114 "namespace absl {\n"
115 "ABSL_NAMESPACE_BEGIN\n"
116 "namespace random_internal {\n"
117 "\n"
118 "const gaussian_distribution_base::Tables\n"
119 " gaussian_distribution_base::zg_ = {\n";
120 FormatArrayContents(os, tables_.x);
121 *os << ",\n";
122 FormatArrayContents(os, tables_.f);
123 *os << "};\n"
124 "\n"
125 "} // namespace random_internal\n"
126 "ABSL_NAMESPACE_END\n"
127 "} // namespace absl\n"
128 "\n"
129 "// clang-format on\n"
130 "// END GENERATED CODE";
131 *os << std::endl;
132 }
133
134 } // namespace random_internal
135 ABSL_NAMESPACE_END
136 } // namespace absl
137
main(int,char **)138 int main(int, char**) {
139 std::cerr << "\nCopy the output to gaussian_distribution.cc" << std::endl;
140 absl::random_internal::TableGenerator generator;
141 generator.Print(&std::cout);
142 return 0;
143 }
144