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
28 template <class T>
29 inline
30 T
sqr(T x)31 sqr(T x)
32 {
33 return x*x;
34 }
35
36 double
f(double x,double a,double m,double b,double c)37 f(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
42 void
test1()43 test1()
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 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 (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(v);
59 }
60 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
74 areas[i] /= S;
75 for (size_t i = 0; i < Np+1; ++i)
76 p[i] /= S;
77 for (size_t i = 0; i < N; ++i)
78 {
79 int k = std::lower_bound(b, b+Np+1, 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
94 void
test2()95 test2()
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 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 (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(v);
111 }
112 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
126 areas[i] /= S;
127 for (size_t i = 0; i < Np+1; ++i)
128 p[i] /= S;
129 for (size_t i = 0; i < N; ++i)
130 {
131 int k = std::lower_bound(b, b+Np+1, 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
146 void
test3()147 test3()
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 size_t Np = sizeof(p) / sizeof(p[0]) - 1;
155 D d(b, b+Np+1, p);
156 const size_t N = 1000000;
157 std::vector<D::result_type> u;
158 for (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(v);
163 }
164 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
178 areas[i] /= S;
179 for (size_t i = 0; i < Np+1; ++i)
180 p[i] /= S;
181 for (size_t i = 0; i < N; ++i)
182 {
183 int k = std::lower_bound(b, b+Np+1, 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
198 void
test4()199 test4()
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 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 (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(v);
215 }
216 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
230 areas[i] /= S;
231 for (size_t i = 0; i < Np+1; ++i)
232 p[i] /= S;
233 for (size_t i = 0; i < N; ++i)
234 {
235 int k = std::lower_bound(b, b+Np+1, 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
251 void
test5()252 test5()
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 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 (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(v);
268 }
269 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
284 areas[i] /= S;
285 for (size_t i = 0; i < Np+1; ++i)
286 p[i] /= S;
287 for (size_t i = 0; i < N; ++i)
288 {
289 int k = std::lower_bound(b, b+Np+1, 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
305 void
test6()306 test6()
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 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 (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(v);
322 }
323 std::sort(u.begin(), 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 (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 (size_t i = 0; i < areas.size(); ++i)
337 areas[i] /= S;
338 for (size_t i = 0; i < Np+1; ++i)
339 p[i] /= S;
340 for (size_t i = 0; i < N; ++i)
341 {
342 int k = std::lower_bound(b, b+Np+1, 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
main(int,char **)357 int main(int, char**)
358 {
359 test1();
360 test2();
361 test3();
362 test4();
363 test5();
364 test6();
365
366 return 0;
367 }
368