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