1//===----------------------------------------------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8//
9// REQUIRES: long_tests
10
11// <random>
12
13// template<class RealType = double>
14// class piecewise_linear_distribution
15
16// template<class _URNG> result_type operator()(_URNG& g);
17
18#include <algorithm>
19#include <cassert>
20#include <cmath>
21#include <cstddef>
22#include <limits>
23#include <random>
24#include <vector>
25
26#include "test_macros.h"
27
28template <class T>
29inline
30T
31sqr(T x)
32{
33 return x*x;
34}
35
36double
37f(double x, double a, double m, double b, double c)
38{
39 return a + m*(sqr(x) - sqr(b))/2 + c*(x-b);
40}
41
42void
43test1()
44{
45 typedef std::piecewise_linear_distribution<> D;
46 typedef std::mt19937_64 G;
47 G g;
48 double b[] = {10, 14, 16, 17};
49 double p[] = {0, 1, 1, 0};
50 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
51 D d(b, b+Np+1, p);
52 const int N = 1000000;
53 std::vector<D::result_type> u;
54 for (std::size_t i = 0; i < N; ++i)
55 {
56 D::result_type v = d(g);
57 assert(d.min() <= v && v < d.max());
58 u.push_back(x: v);
59 }
60 std::sort(first: u.begin(), last: u.end());
61 std::ptrdiff_t kp = -1;
62 double a = std::numeric_limits<double>::quiet_NaN();
63 double m = std::numeric_limits<double>::quiet_NaN();
64 double bk = std::numeric_limits<double>::quiet_NaN();
65 double c = std::numeric_limits<double>::quiet_NaN();
66 std::vector<double> areas(Np);
67 double S = 0;
68 for (std::size_t i = 0; i < areas.size(); ++i)
69 {
70 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
71 S += areas[i];
72 }
73 for (std::size_t i = 0; i < areas.size(); ++i)
74 areas[i] /= S;
75 for (std::size_t i = 0; i < Np+1; ++i)
76 p[i] /= S;
77 for (std::size_t i = 0; i < N; ++i)
78 {
79 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
80 if (k != kp) {
81 a = 0;
82 for (int j = 0; j < k; ++j)
83 a += areas[j];
84 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
85 bk = b[k];
86 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
87 kp = k;
88 }
89 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
90 }
91}
92
93void
94test2()
95{
96 typedef std::piecewise_linear_distribution<> D;
97 typedef std::mt19937_64 G;
98 G g;
99 double b[] = {10, 14, 16, 17};
100 double p[] = {0, 0, 1, 0};
101 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
102 D d(b, b+Np+1, p);
103 const int N = 1000000;
104 std::vector<D::result_type> u;
105 for (std::size_t i = 0; i < N; ++i)
106 {
107 D::result_type v = d(g);
108 assert(d.min() <= v && v < d.max());
109 u.push_back(x: v);
110 }
111 std::sort(first: u.begin(), last: u.end());
112 std::ptrdiff_t kp = -1;
113 double a = std::numeric_limits<double>::quiet_NaN();
114 double m = std::numeric_limits<double>::quiet_NaN();
115 double bk = std::numeric_limits<double>::quiet_NaN();
116 double c = std::numeric_limits<double>::quiet_NaN();
117 std::vector<double> areas(Np);
118 double S = 0;
119 for (std::size_t i = 0; i < areas.size(); ++i)
120 {
121 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
122 S += areas[i];
123 }
124 for (std::size_t i = 0; i < areas.size(); ++i)
125 areas[i] /= S;
126 for (std::size_t i = 0; i < Np+1; ++i)
127 p[i] /= S;
128 for (std::size_t i = 0; i < N; ++i)
129 {
130 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
131 if (k != kp) {
132 a = 0;
133 for (int j = 0; j < k; ++j)
134 a += areas[j];
135 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
136 bk = b[k];
137 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
138 kp = k;
139 }
140 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
141 }
142}
143
144void
145test3()
146{
147 typedef std::piecewise_linear_distribution<> D;
148 typedef std::mt19937_64 G;
149 G g;
150 double b[] = {10, 14, 16, 17};
151 double p[] = {1, 0, 0, 0};
152 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
153 D d(b, b+Np+1, p);
154 const std::size_t N = 1000000;
155 std::vector<D::result_type> u;
156 for (std::size_t i = 0; i < N; ++i)
157 {
158 D::result_type v = d(g);
159 assert(d.min() <= v && v < d.max());
160 u.push_back(x: v);
161 }
162 std::sort(first: u.begin(), last: u.end());
163 std::ptrdiff_t kp = -1;
164 double a = std::numeric_limits<double>::quiet_NaN();
165 double m = std::numeric_limits<double>::quiet_NaN();
166 double bk = std::numeric_limits<double>::quiet_NaN();
167 double c = std::numeric_limits<double>::quiet_NaN();
168 std::vector<double> areas(Np);
169 double S = 0;
170 for (std::size_t i = 0; i < areas.size(); ++i)
171 {
172 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
173 S += areas[i];
174 }
175 for (std::size_t i = 0; i < areas.size(); ++i)
176 areas[i] /= S;
177 for (std::size_t i = 0; i < Np+1; ++i)
178 p[i] /= S;
179 for (std::size_t i = 0; i < N; ++i)
180 {
181 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
182 if (k != kp) {
183 a = 0;
184 for (int j = 0; j < k; ++j)
185 a += areas[j];
186 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
187 bk = b[k];
188 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
189 kp = k;
190 }
191 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
192 }
193}
194
195void
196test4()
197{
198 typedef std::piecewise_linear_distribution<> D;
199 typedef std::mt19937_64 G;
200 G g;
201 double b[] = {10, 14, 16};
202 double p[] = {0, 1, 0};
203 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
204 D d(b, b+Np+1, p);
205 const int N = 1000000;
206 std::vector<D::result_type> u;
207 for (std::size_t i = 0; i < N; ++i)
208 {
209 D::result_type v = d(g);
210 assert(d.min() <= v && v < d.max());
211 u.push_back(x: v);
212 }
213 std::sort(first: u.begin(), last: u.end());
214 std::ptrdiff_t kp = -1;
215 double a = std::numeric_limits<double>::quiet_NaN();
216 double m = std::numeric_limits<double>::quiet_NaN();
217 double bk = std::numeric_limits<double>::quiet_NaN();
218 double c = std::numeric_limits<double>::quiet_NaN();
219 std::vector<double> areas(Np);
220 double S = 0;
221 for (std::size_t i = 0; i < areas.size(); ++i)
222 {
223 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
224 S += areas[i];
225 }
226 for (std::size_t i = 0; i < areas.size(); ++i)
227 areas[i] /= S;
228 for (std::size_t i = 0; i < Np+1; ++i)
229 p[i] /= S;
230 for (std::size_t i = 0; i < N; ++i)
231 {
232 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
233 if (k != kp) {
234 a = 0;
235 for (int j = 0; j < k; ++j)
236 a += areas[j];
237 assert(k < static_cast<int>(Np));
238 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
239 bk = b[k];
240 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
241 kp = k;
242 }
243 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
244 }
245}
246
247void
248test5()
249{
250 typedef std::piecewise_linear_distribution<> D;
251 typedef std::mt19937_64 G;
252 G g;
253 double b[] = {10, 14};
254 double p[] = {1, 1};
255 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
256 D d(b, b+Np+1, p);
257 const int N = 1000000;
258 std::vector<D::result_type> u;
259 for (std::size_t i = 0; i < N; ++i)
260 {
261 D::result_type v = d(g);
262 assert(d.min() <= v && v < d.max());
263 u.push_back(x: v);
264 }
265 std::sort(first: u.begin(), last: u.end());
266 std::ptrdiff_t kp = -1;
267 double a = std::numeric_limits<double>::quiet_NaN();
268 double m = std::numeric_limits<double>::quiet_NaN();
269 double bk = std::numeric_limits<double>::quiet_NaN();
270 double c = std::numeric_limits<double>::quiet_NaN();
271 std::vector<double> areas(Np);
272 double S = 0;
273 for (std::size_t i = 0; i < areas.size(); ++i)
274 {
275 assert(i < Np);
276 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
277 S += areas[i];
278 }
279 for (std::size_t i = 0; i < areas.size(); ++i)
280 areas[i] /= S;
281 for (std::size_t i = 0; i < Np+1; ++i)
282 p[i] /= S;
283 for (std::size_t i = 0; i < N; ++i)
284 {
285 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
286 if (k != kp) {
287 a = 0;
288 for (int j = 0; j < k; ++j)
289 a += areas[j];
290 assert(k < static_cast<int>(Np));
291 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
292 bk = b[k];
293 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
294 kp = k;
295 }
296 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
297 }
298}
299
300void
301test6()
302{
303 typedef std::piecewise_linear_distribution<> D;
304 typedef std::mt19937_64 G;
305 G g;
306 double b[] = {10, 14, 16, 17};
307 double p[] = {25, 62.5, 12.5, 0};
308 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
309 D d(b, b+Np+1, p);
310 const int N = 1000000;
311 std::vector<D::result_type> u;
312 for (std::size_t i = 0; i < N; ++i)
313 {
314 D::result_type v = d(g);
315 assert(d.min() <= v && v < d.max());
316 u.push_back(x: v);
317 }
318 std::sort(first: u.begin(), last: u.end());
319 std::ptrdiff_t kp = -1;
320 double a = std::numeric_limits<double>::quiet_NaN();
321 double m = std::numeric_limits<double>::quiet_NaN();
322 double bk = std::numeric_limits<double>::quiet_NaN();
323 double c = std::numeric_limits<double>::quiet_NaN();
324 std::vector<double> areas(Np);
325 double S = 0;
326 for (std::size_t i = 0; i < areas.size(); ++i)
327 {
328 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
329 S += areas[i];
330 }
331 for (std::size_t i = 0; i < areas.size(); ++i)
332 areas[i] /= S;
333 for (std::size_t i = 0; i < Np+1; ++i)
334 p[i] /= S;
335 for (std::size_t i = 0; i < N; ++i)
336 {
337 std::ptrdiff_t k = std::lower_bound(first: b, last: b + Np + 1, val: u[i]) - b - 1;
338 if (k != kp) {
339 a = 0;
340 for (int j = 0; j < k; ++j)
341 a += areas[j];
342 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);
343 bk = b[k];
344 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);
345 kp = k;
346 }
347 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);
348 }
349}
350
351int main(int, char**)
352{
353 test1();
354 test2();
355 test3();
356 test4();
357 test5();
358 test6();
359
360 return 0;
361}
362

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.samp/rand.dist.samp.plinear/eval.pass.cpp