xref: /llvm-project/libcxx/test/std/numerics/rand/rand.dist/rand.dist.uni/rand.dist.uni.real/eval.pass.cpp (revision 09e3a360581dc36d0820d3fb6da9bd7cfed87b5d)
1eb1c5037SNikolas Klauser //===----------------------------------------------------------------------===//
2eb1c5037SNikolas Klauser //
3eb1c5037SNikolas Klauser // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4eb1c5037SNikolas Klauser // See https://llvm.org/LICENSE.txt for license information.
5eb1c5037SNikolas Klauser // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6eb1c5037SNikolas Klauser //
7eb1c5037SNikolas Klauser //===----------------------------------------------------------------------===//
8eb1c5037SNikolas Klauser //
9eb1c5037SNikolas Klauser // REQUIRES: long_tests
10eb1c5037SNikolas Klauser 
11eb1c5037SNikolas Klauser // <random>
12eb1c5037SNikolas Klauser 
13eb1c5037SNikolas Klauser // template<class RealType = double>
14eb1c5037SNikolas Klauser // class uniform_real_distribution
15eb1c5037SNikolas Klauser 
16eb1c5037SNikolas Klauser // template<class _URNG> result_type operator()(_URNG& g);
17eb1c5037SNikolas Klauser 
18eb1c5037SNikolas Klauser #include <random>
19eb1c5037SNikolas Klauser #include <cassert>
20*09e3a360SLouis Dionne #include <cmath>
21eb1c5037SNikolas Klauser #include <cstddef>
22*09e3a360SLouis Dionne #include <numeric>
23*09e3a360SLouis Dionne #include <vector>
24eb1c5037SNikolas Klauser 
25eb1c5037SNikolas Klauser #include "test_macros.h"
26eb1c5037SNikolas Klauser 
27eb1c5037SNikolas Klauser template <class T>
28eb1c5037SNikolas Klauser inline
29eb1c5037SNikolas Klauser T
30eb1c5037SNikolas Klauser sqr(T x)
31eb1c5037SNikolas Klauser {
32eb1c5037SNikolas Klauser     return x * x;
33eb1c5037SNikolas Klauser }
34eb1c5037SNikolas Klauser 
35eb1c5037SNikolas Klauser int main(int, char**)
36eb1c5037SNikolas Klauser {
37eb1c5037SNikolas Klauser     {
38eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
39eb1c5037SNikolas Klauser         typedef std::minstd_rand0 G;
40eb1c5037SNikolas Klauser         G g;
41eb1c5037SNikolas Klauser         D d;
42eb1c5037SNikolas Klauser         const int N = 100000;
43eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
44eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
45eb1c5037SNikolas Klauser         {
46eb1c5037SNikolas Klauser             D::result_type v = d(g);
47eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
48eb1c5037SNikolas Klauser             u.push_back(v);
49eb1c5037SNikolas Klauser         }
50eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
51eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
52eb1c5037SNikolas Klauser         D::result_type var = 0;
53eb1c5037SNikolas Klauser         D::result_type skew = 0;
54eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
55eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
56eb1c5037SNikolas Klauser         {
57eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
58eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
59eb1c5037SNikolas Klauser             var += d2;
60eb1c5037SNikolas Klauser             skew += dbl * d2;
61eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
62eb1c5037SNikolas Klauser         }
63eb1c5037SNikolas Klauser         var /= u.size();
64eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
65eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
66eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
67eb1c5037SNikolas Klauser         kurtosis -= 3;
68eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
69eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
70eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
71eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
72eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
73eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
74eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
75eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
76eb1c5037SNikolas Klauser     }
77eb1c5037SNikolas Klauser     {
78eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
79eb1c5037SNikolas Klauser         typedef std::minstd_rand G;
80eb1c5037SNikolas Klauser         G g;
81eb1c5037SNikolas Klauser         D d;
82eb1c5037SNikolas Klauser         const int N = 100000;
83eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
84eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
85eb1c5037SNikolas Klauser         {
86eb1c5037SNikolas Klauser             D::result_type v = d(g);
87eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
88eb1c5037SNikolas Klauser             u.push_back(v);
89eb1c5037SNikolas Klauser         }
90eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
91eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
92eb1c5037SNikolas Klauser         D::result_type var = 0;
93eb1c5037SNikolas Klauser         D::result_type skew = 0;
94eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
95eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
96eb1c5037SNikolas Klauser         {
97eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
98eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
99eb1c5037SNikolas Klauser             var += d2;
100eb1c5037SNikolas Klauser             skew += dbl * d2;
101eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
102eb1c5037SNikolas Klauser         }
103eb1c5037SNikolas Klauser         var /= u.size();
104eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
105eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
106eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
107eb1c5037SNikolas Klauser         kurtosis -= 3;
108eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
109eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
110eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
111eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
112eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
113eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
114eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
115eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
116eb1c5037SNikolas Klauser     }
117eb1c5037SNikolas Klauser     {
118eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
119eb1c5037SNikolas Klauser         typedef std::mt19937 G;
120eb1c5037SNikolas Klauser         G g;
121eb1c5037SNikolas Klauser         D d;
122eb1c5037SNikolas Klauser         const int N = 100000;
123eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
124eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
125eb1c5037SNikolas Klauser         {
126eb1c5037SNikolas Klauser             D::result_type v = d(g);
127eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
128eb1c5037SNikolas Klauser             u.push_back(v);
129eb1c5037SNikolas Klauser         }
130eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
131eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
132eb1c5037SNikolas Klauser         D::result_type var = 0;
133eb1c5037SNikolas Klauser         D::result_type skew = 0;
134eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
135eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
136eb1c5037SNikolas Klauser         {
137eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
138eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
139eb1c5037SNikolas Klauser             var += d2;
140eb1c5037SNikolas Klauser             skew += dbl * d2;
141eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
142eb1c5037SNikolas Klauser         }
143eb1c5037SNikolas Klauser         var /= u.size();
144eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
145eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
146eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
147eb1c5037SNikolas Klauser         kurtosis -= 3;
148eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
149eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
150eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
151eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
152eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
153eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
154eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
155eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
156eb1c5037SNikolas Klauser     }
157eb1c5037SNikolas Klauser     {
158eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
159eb1c5037SNikolas Klauser         typedef std::mt19937_64 G;
160eb1c5037SNikolas Klauser         G g;
161eb1c5037SNikolas Klauser         D d;
162eb1c5037SNikolas Klauser         const int N = 100000;
163eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
164eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
165eb1c5037SNikolas Klauser         {
166eb1c5037SNikolas Klauser             D::result_type v = d(g);
167eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
168eb1c5037SNikolas Klauser             u.push_back(v);
169eb1c5037SNikolas Klauser         }
170eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
171eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
172eb1c5037SNikolas Klauser         D::result_type var = 0;
173eb1c5037SNikolas Klauser         D::result_type skew = 0;
174eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
175eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
176eb1c5037SNikolas Klauser         {
177eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
178eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
179eb1c5037SNikolas Klauser             var += d2;
180eb1c5037SNikolas Klauser             skew += dbl * d2;
181eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
182eb1c5037SNikolas Klauser         }
183eb1c5037SNikolas Klauser         var /= u.size();
184eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
185eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
186eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
187eb1c5037SNikolas Klauser         kurtosis -= 3;
188eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
189eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
190eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
191eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
192eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
193eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
194eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
195eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
196eb1c5037SNikolas Klauser     }
197eb1c5037SNikolas Klauser     {
198eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
199eb1c5037SNikolas Klauser         typedef std::ranlux24_base G;
200eb1c5037SNikolas Klauser         G g;
201eb1c5037SNikolas Klauser         D d;
202eb1c5037SNikolas Klauser         const int N = 100000;
203eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
204eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
205eb1c5037SNikolas Klauser         {
206eb1c5037SNikolas Klauser             D::result_type v = d(g);
207eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
208eb1c5037SNikolas Klauser             u.push_back(v);
209eb1c5037SNikolas Klauser         }
210eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
211eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
212eb1c5037SNikolas Klauser         D::result_type var = 0;
213eb1c5037SNikolas Klauser         D::result_type skew = 0;
214eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
215eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
216eb1c5037SNikolas Klauser         {
217eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
218eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
219eb1c5037SNikolas Klauser             var += d2;
220eb1c5037SNikolas Klauser             skew += dbl * d2;
221eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
222eb1c5037SNikolas Klauser         }
223eb1c5037SNikolas Klauser         var /= u.size();
224eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
225eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
226eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
227eb1c5037SNikolas Klauser         kurtosis -= 3;
228eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
229eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
230eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
231eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
232eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
233eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
234eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.02);
235eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
236eb1c5037SNikolas Klauser     }
237eb1c5037SNikolas Klauser     {
238eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
239eb1c5037SNikolas Klauser         typedef std::ranlux48_base G;
240eb1c5037SNikolas Klauser         G g;
241eb1c5037SNikolas Klauser         D d;
242eb1c5037SNikolas Klauser         const int N = 100000;
243eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
244eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
245eb1c5037SNikolas Klauser         {
246eb1c5037SNikolas Klauser             D::result_type v = d(g);
247eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
248eb1c5037SNikolas Klauser             u.push_back(v);
249eb1c5037SNikolas Klauser         }
250eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
251eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
252eb1c5037SNikolas Klauser         D::result_type var = 0;
253eb1c5037SNikolas Klauser         D::result_type skew = 0;
254eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
255eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
256eb1c5037SNikolas Klauser         {
257eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
258eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
259eb1c5037SNikolas Klauser             var += d2;
260eb1c5037SNikolas Klauser             skew += dbl * d2;
261eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
262eb1c5037SNikolas Klauser         }
263eb1c5037SNikolas Klauser         var /= u.size();
264eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
265eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
266eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
267eb1c5037SNikolas Klauser         kurtosis -= 3;
268eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
269eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
270eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
271eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
272eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
273eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
274eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
275eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
276eb1c5037SNikolas Klauser     }
277eb1c5037SNikolas Klauser     {
278eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
279eb1c5037SNikolas Klauser         typedef std::ranlux24 G;
280eb1c5037SNikolas Klauser         G g;
281eb1c5037SNikolas Klauser         D d;
282eb1c5037SNikolas Klauser         const int N = 100000;
283eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
284eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
285eb1c5037SNikolas Klauser         {
286eb1c5037SNikolas Klauser             D::result_type v = d(g);
287eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
288eb1c5037SNikolas Klauser             u.push_back(v);
289eb1c5037SNikolas Klauser         }
290eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
291eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
292eb1c5037SNikolas Klauser         D::result_type var = 0;
293eb1c5037SNikolas Klauser         D::result_type skew = 0;
294eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
295eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
296eb1c5037SNikolas Klauser         {
297eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
298eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
299eb1c5037SNikolas Klauser             var += d2;
300eb1c5037SNikolas Klauser             skew += dbl * d2;
301eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
302eb1c5037SNikolas Klauser         }
303eb1c5037SNikolas Klauser         var /= u.size();
304eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
305eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
306eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
307eb1c5037SNikolas Klauser         kurtosis -= 3;
308eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
309eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
310eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
311eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
312eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
313eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
314eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
315eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
316eb1c5037SNikolas Klauser     }
317eb1c5037SNikolas Klauser     {
318eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
319eb1c5037SNikolas Klauser         typedef std::ranlux48 G;
320eb1c5037SNikolas Klauser         G g;
321eb1c5037SNikolas Klauser         D d;
322eb1c5037SNikolas Klauser         const int N = 100000;
323eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
324eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
325eb1c5037SNikolas Klauser         {
326eb1c5037SNikolas Klauser             D::result_type v = d(g);
327eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
328eb1c5037SNikolas Klauser             u.push_back(v);
329eb1c5037SNikolas Klauser         }
330eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
331eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
332eb1c5037SNikolas Klauser         D::result_type var = 0;
333eb1c5037SNikolas Klauser         D::result_type skew = 0;
334eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
335eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
336eb1c5037SNikolas Klauser         {
337eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
338eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
339eb1c5037SNikolas Klauser             var += d2;
340eb1c5037SNikolas Klauser             skew += dbl * d2;
341eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
342eb1c5037SNikolas Klauser         }
343eb1c5037SNikolas Klauser         var /= u.size();
344eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
345eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
346eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
347eb1c5037SNikolas Klauser         kurtosis -= 3;
348eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
349eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
350eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
351eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
352eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
353eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
354eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
355eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
356eb1c5037SNikolas Klauser     }
357eb1c5037SNikolas Klauser     {
358eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
359eb1c5037SNikolas Klauser         typedef std::knuth_b G;
360eb1c5037SNikolas Klauser         G g;
361eb1c5037SNikolas Klauser         D d;
362eb1c5037SNikolas Klauser         const int N = 100000;
363eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
364eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
365eb1c5037SNikolas Klauser         {
366eb1c5037SNikolas Klauser             D::result_type v = d(g);
367eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
368eb1c5037SNikolas Klauser             u.push_back(v);
369eb1c5037SNikolas Klauser         }
370eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
371eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
372eb1c5037SNikolas Klauser         D::result_type var = 0;
373eb1c5037SNikolas Klauser         D::result_type skew = 0;
374eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
375eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
376eb1c5037SNikolas Klauser         {
377eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
378eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
379eb1c5037SNikolas Klauser             var += d2;
380eb1c5037SNikolas Klauser             skew += dbl * d2;
381eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
382eb1c5037SNikolas Klauser         }
383eb1c5037SNikolas Klauser         var /= u.size();
384eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
385eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
386eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
387eb1c5037SNikolas Klauser         kurtosis -= 3;
388eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
389eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
390eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
391eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
392eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
393eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
394eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
395eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
396eb1c5037SNikolas Klauser     }
397eb1c5037SNikolas Klauser     {
398eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
399eb1c5037SNikolas Klauser         typedef std::minstd_rand G;
400eb1c5037SNikolas Klauser         G g;
401eb1c5037SNikolas Klauser         D d(-1, 1);
402eb1c5037SNikolas Klauser         const int N = 100000;
403eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
404eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
405eb1c5037SNikolas Klauser         {
406eb1c5037SNikolas Klauser             D::result_type v = d(g);
407eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
408eb1c5037SNikolas Klauser             u.push_back(v);
409eb1c5037SNikolas Klauser         }
410eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
411eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
412eb1c5037SNikolas Klauser         D::result_type var = 0;
413eb1c5037SNikolas Klauser         D::result_type skew = 0;
414eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
415eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
416eb1c5037SNikolas Klauser         {
417eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
418eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
419eb1c5037SNikolas Klauser             var += d2;
420eb1c5037SNikolas Klauser             skew += dbl * d2;
421eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
422eb1c5037SNikolas Klauser         }
423eb1c5037SNikolas Klauser         var /= u.size();
424eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
425eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
426eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
427eb1c5037SNikolas Klauser         kurtosis -= 3;
428eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
429eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
430eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
431eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
432eb1c5037SNikolas Klauser         assert(std::abs(mean - x_mean) < 0.01);
433eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
434eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
435eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
436eb1c5037SNikolas Klauser     }
437eb1c5037SNikolas Klauser     {
438eb1c5037SNikolas Klauser         typedef std::uniform_real_distribution<> D;
439eb1c5037SNikolas Klauser         typedef std::minstd_rand G;
440eb1c5037SNikolas Klauser         G g;
441eb1c5037SNikolas Klauser         D d(5.5, 25);
442eb1c5037SNikolas Klauser         const int N = 100000;
443eb1c5037SNikolas Klauser         std::vector<D::result_type> u;
444eb1c5037SNikolas Klauser         for (int i = 0; i < N; ++i)
445eb1c5037SNikolas Klauser         {
446eb1c5037SNikolas Klauser             D::result_type v = d(g);
447eb1c5037SNikolas Klauser             assert(d.a() <= v && v < d.b());
448eb1c5037SNikolas Klauser             u.push_back(v);
449eb1c5037SNikolas Klauser         }
450eb1c5037SNikolas Klauser         D::result_type mean = std::accumulate(u.begin(), u.end(),
451eb1c5037SNikolas Klauser                                               D::result_type(0)) / u.size();
452eb1c5037SNikolas Klauser         D::result_type var = 0;
453eb1c5037SNikolas Klauser         D::result_type skew = 0;
454eb1c5037SNikolas Klauser         D::result_type kurtosis = 0;
455eb1c5037SNikolas Klauser         for (std::size_t i = 0; i < u.size(); ++i)
456eb1c5037SNikolas Klauser         {
457eb1c5037SNikolas Klauser             D::result_type dbl = (u[i] - mean);
458eb1c5037SNikolas Klauser             D::result_type d2 = sqr(dbl);
459eb1c5037SNikolas Klauser             var += d2;
460eb1c5037SNikolas Klauser             skew += dbl * d2;
461eb1c5037SNikolas Klauser             kurtosis += d2 * d2;
462eb1c5037SNikolas Klauser         }
463eb1c5037SNikolas Klauser         var /= u.size();
464eb1c5037SNikolas Klauser         D::result_type dev = std::sqrt(var);
465eb1c5037SNikolas Klauser         skew /= u.size() * dev * var;
466eb1c5037SNikolas Klauser         kurtosis /= u.size() * var * var;
467eb1c5037SNikolas Klauser         kurtosis -= 3;
468eb1c5037SNikolas Klauser         D::result_type x_mean = (d.a() + d.b()) / 2;
469eb1c5037SNikolas Klauser         D::result_type x_var = sqr(d.b() - d.a()) / 12;
470eb1c5037SNikolas Klauser         D::result_type x_skew = 0;
471eb1c5037SNikolas Klauser         D::result_type x_kurtosis = -6./5;
472eb1c5037SNikolas Klauser         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
473eb1c5037SNikolas Klauser         assert(std::abs((var - x_var) / x_var) < 0.01);
474eb1c5037SNikolas Klauser         assert(std::abs(skew - x_skew) < 0.01);
475eb1c5037SNikolas Klauser         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
476eb1c5037SNikolas Klauser     }
477eb1c5037SNikolas Klauser 
478eb1c5037SNikolas Klauser   return 0;
479eb1c5037SNikolas Klauser }
480