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