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 <random>
19#include <algorithm>
20#include <vector>
21#include <iterator>
22#include <numeric>
23#include <cassert>
24#include <limits>
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 int 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 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
80 if (k != kp)
81 {
82 a = 0;
83 for (int j = 0; j < k; ++j)
84 a += areas[j];
85 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
86 bk = b[k];
87 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
88 kp = k;
89 }
90 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
91 }
92}
93
94void
95test2()
96{
97 typedef std::piecewise_linear_distribution<> D;
98 typedef std::mt19937_64 G;
99 G g;
100 double b[] = {10, 14, 16, 17};
101 double p[] = {0, 0, 1, 0};
102 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
103 D d(b, b+Np+1, p);
104 const int N = 1000000;
105 std::vector<D::result_type> u;
106 for (std::size_t i = 0; i < N; ++i)
107 {
108 D::result_type v = d(g);
109 assert(d.min() <= v && v < d.max());
110 u.push_back(x: v);
111 }
112 std::sort(first: u.begin(), last: u.end());
113 int kp = -1;
114 double a = std::numeric_limits<double>::quiet_NaN();
115 double m = std::numeric_limits<double>::quiet_NaN();
116 double bk = std::numeric_limits<double>::quiet_NaN();
117 double c = std::numeric_limits<double>::quiet_NaN();
118 std::vector<double> areas(Np);
119 double S = 0;
120 for (std::size_t i = 0; i < areas.size(); ++i)
121 {
122 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
123 S += areas[i];
124 }
125 for (std::size_t i = 0; i < areas.size(); ++i)
126 areas[i] /= S;
127 for (std::size_t i = 0; i < Np+1; ++i)
128 p[i] /= S;
129 for (std::size_t i = 0; i < N; ++i)
130 {
131 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
132 if (k != kp)
133 {
134 a = 0;
135 for (int j = 0; j < k; ++j)
136 a += areas[j];
137 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
138 bk = b[k];
139 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
140 kp = k;
141 }
142 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
143 }
144}
145
146void
147test3()
148{
149 typedef std::piecewise_linear_distribution<> D;
150 typedef std::mt19937_64 G;
151 G g;
152 double b[] = {10, 14, 16, 17};
153 double p[] = {1, 0, 0, 0};
154 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
155 D d(b, b+Np+1, p);
156 const std::size_t N = 1000000;
157 std::vector<D::result_type> u;
158 for (std::size_t i = 0; i < N; ++i)
159 {
160 D::result_type v = d(g);
161 assert(d.min() <= v && v < d.max());
162 u.push_back(x: v);
163 }
164 std::sort(first: u.begin(), last: u.end());
165 int kp = -1;
166 double a = std::numeric_limits<double>::quiet_NaN();
167 double m = std::numeric_limits<double>::quiet_NaN();
168 double bk = std::numeric_limits<double>::quiet_NaN();
169 double c = std::numeric_limits<double>::quiet_NaN();
170 std::vector<double> areas(Np);
171 double S = 0;
172 for (std::size_t i = 0; i < areas.size(); ++i)
173 {
174 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
175 S += areas[i];
176 }
177 for (std::size_t i = 0; i < areas.size(); ++i)
178 areas[i] /= S;
179 for (std::size_t i = 0; i < Np+1; ++i)
180 p[i] /= S;
181 for (std::size_t i = 0; i < N; ++i)
182 {
183 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
184 if (k != kp)
185 {
186 a = 0;
187 for (int j = 0; j < k; ++j)
188 a += areas[j];
189 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
190 bk = b[k];
191 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
192 kp = k;
193 }
194 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
195 }
196}
197
198void
199test4()
200{
201 typedef std::piecewise_linear_distribution<> D;
202 typedef std::mt19937_64 G;
203 G g;
204 double b[] = {10, 14, 16};
205 double p[] = {0, 1, 0};
206 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
207 D d(b, b+Np+1, p);
208 const int N = 1000000;
209 std::vector<D::result_type> u;
210 for (std::size_t i = 0; i < N; ++i)
211 {
212 D::result_type v = d(g);
213 assert(d.min() <= v && v < d.max());
214 u.push_back(x: v);
215 }
216 std::sort(first: u.begin(), last: u.end());
217 int kp = -1;
218 double a = std::numeric_limits<double>::quiet_NaN();
219 double m = std::numeric_limits<double>::quiet_NaN();
220 double bk = std::numeric_limits<double>::quiet_NaN();
221 double c = std::numeric_limits<double>::quiet_NaN();
222 std::vector<double> areas(Np);
223 double S = 0;
224 for (std::size_t i = 0; i < areas.size(); ++i)
225 {
226 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
227 S += areas[i];
228 }
229 for (std::size_t i = 0; i < areas.size(); ++i)
230 areas[i] /= S;
231 for (std::size_t i = 0; i < Np+1; ++i)
232 p[i] /= S;
233 for (std::size_t i = 0; i < N; ++i)
234 {
235 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
236 if (k != kp)
237 {
238 a = 0;
239 for (int j = 0; j < k; ++j)
240 a += areas[j];
241 assert(k < static_cast<int>(Np));
242 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
243 bk = b[k];
244 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
245 kp = k;
246 }
247 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
248 }
249}
250
251void
252test5()
253{
254 typedef std::piecewise_linear_distribution<> D;
255 typedef std::mt19937_64 G;
256 G g;
257 double b[] = {10, 14};
258 double p[] = {1, 1};
259 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
260 D d(b, b+Np+1, p);
261 const int N = 1000000;
262 std::vector<D::result_type> u;
263 for (std::size_t i = 0; i < N; ++i)
264 {
265 D::result_type v = d(g);
266 assert(d.min() <= v && v < d.max());
267 u.push_back(x: v);
268 }
269 std::sort(first: u.begin(), last: u.end());
270 int kp = -1;
271 double a = std::numeric_limits<double>::quiet_NaN();
272 double m = std::numeric_limits<double>::quiet_NaN();
273 double bk = std::numeric_limits<double>::quiet_NaN();
274 double c = std::numeric_limits<double>::quiet_NaN();
275 std::vector<double> areas(Np);
276 double S = 0;
277 for (std::size_t i = 0; i < areas.size(); ++i)
278 {
279 assert(i < Np);
280 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
281 S += areas[i];
282 }
283 for (std::size_t i = 0; i < areas.size(); ++i)
284 areas[i] /= S;
285 for (std::size_t i = 0; i < Np+1; ++i)
286 p[i] /= S;
287 for (std::size_t i = 0; i < N; ++i)
288 {
289 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
290 if (k != kp)
291 {
292 a = 0;
293 for (int j = 0; j < k; ++j)
294 a += areas[j];
295 assert(k < static_cast<int>(Np));
296 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
297 bk = b[k];
298 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
299 kp = k;
300 }
301 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
302 }
303}
304
305void
306test6()
307{
308 typedef std::piecewise_linear_distribution<> D;
309 typedef std::mt19937_64 G;
310 G g;
311 double b[] = {10, 14, 16, 17};
312 double p[] = {25, 62.5, 12.5, 0};
313 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;
314 D d(b, b+Np+1, p);
315 const int N = 1000000;
316 std::vector<D::result_type> u;
317 for (std::size_t i = 0; i < N; ++i)
318 {
319 D::result_type v = d(g);
320 assert(d.min() <= v && v < d.max());
321 u.push_back(x: v);
322 }
323 std::sort(first: u.begin(), last: u.end());
324 int kp = -1;
325 double a = std::numeric_limits<double>::quiet_NaN();
326 double m = std::numeric_limits<double>::quiet_NaN();
327 double bk = std::numeric_limits<double>::quiet_NaN();
328 double c = std::numeric_limits<double>::quiet_NaN();
329 std::vector<double> areas(Np);
330 double S = 0;
331 for (std::size_t i = 0; i < areas.size(); ++i)
332 {
333 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
334 S += areas[i];
335 }
336 for (std::size_t i = 0; i < areas.size(); ++i)
337 areas[i] /= S;
338 for (std::size_t i = 0; i < Np+1; ++i)
339 p[i] /= S;
340 for (std::size_t i = 0; i < N; ++i)
341 {
342 int k = std::lower_bound(first: b, last: b+Np+1, val: u[i]) - b - 1;
343 if (k != kp)
344 {
345 a = 0;
346 for (int j = 0; j < k; ++j)
347 a += areas[j];
348 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
349 bk = b[k];
350 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
351 kp = k;
352 }
353 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
354 }
355}
356
357int main(int, char**)
358{
359 test1();
360 test2();
361 test3();
362 test4();
363 test5();
364 test6();
365
366 return 0;
367}
368

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