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 <vector>
23 #include <iterator>
24 #include <numeric>
25 #include <cassert>
26
27 template <class T>
28 inline
29 T
sqr(T x)30 sqr(T x)
31 {
32 return x*x;
33 }
34
35 double
f(double x,double a,double m,double b,double c)36 f(double x, double a, double m, double b, double c)
37 {
38 return a + m*(sqr(x) - sqr(b))/2 + c*(x-b);
39 }
40
main()41 int main()
42 {
43 {
44 typedef std::piecewise_linear_distribution<> D;
45 typedef D::param_type P;
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 (int 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;
63 double m;
64 double bk;
65 double c;
66 std::vector<double> areas(Np);
67 double S = 0;
68 for (int 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 (int i = 0; i < areas.size(); ++i)
74 areas[i] /= S;
75 for (int i = 0; i < Np+1; ++i)
76 p[i] /= S;
77 for (int 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 typedef std::piecewise_linear_distribution<> D;
95 typedef D::param_type P;
96 typedef std::mt19937_64 G;
97 G g;
98 double b[] = {10, 14, 16, 17};
99 double p[] = {0, 0, 1, 0};
100 const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
101 D d(b, b+Np+1, p);
102 const int N = 1000000;
103 std::vector<D::result_type> u;
104 for (int i = 0; i < N; ++i)
105 {
106 D::result_type v = d(g);
107 assert(d.min() <= v && v < d.max());
108 u.push_back(v);
109 }
110 std::sort(u.begin(), u.end());
111 int kp = -1;
112 double a;
113 double m;
114 double bk;
115 double c;
116 std::vector<double> areas(Np);
117 double S = 0;
118 for (int i = 0; i < areas.size(); ++i)
119 {
120 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
121 S += areas[i];
122 }
123 for (int i = 0; i < areas.size(); ++i)
124 areas[i] /= S;
125 for (int i = 0; i < Np+1; ++i)
126 p[i] /= S;
127 for (int i = 0; i < N; ++i)
128 {
129 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
130 if (k != kp)
131 {
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) < .001);
141 }
142 }
143 {
144 typedef std::piecewise_linear_distribution<> D;
145 typedef D::param_type P;
146 typedef std::mt19937_64 G;
147 G g;
148 double b[] = {10, 14, 16, 17};
149 double p[] = {1, 0, 0, 0};
150 const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
151 D d(b, b+Np+1, p);
152 const int N = 1000000;
153 std::vector<D::result_type> u;
154 for (int i = 0; i < N; ++i)
155 {
156 D::result_type v = d(g);
157 assert(d.min() <= v && v < d.max());
158 u.push_back(v);
159 }
160 std::sort(u.begin(), u.end());
161 int kp = -1;
162 double a;
163 double m;
164 double bk;
165 double c;
166 std::vector<double> areas(Np);
167 double S = 0;
168 for (int i = 0; i < areas.size(); ++i)
169 {
170 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
171 S += areas[i];
172 }
173 for (int i = 0; i < areas.size(); ++i)
174 areas[i] /= S;
175 for (int i = 0; i < Np+1; ++i)
176 p[i] /= S;
177 for (int i = 0; i < N; ++i)
178 {
179 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
180 if (k != kp)
181 {
182 a = 0;
183 for (int j = 0; j < k; ++j)
184 a += areas[j];
185 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
186 bk = b[k];
187 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
188 kp = k;
189 }
190 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
191 }
192 }
193 {
194 typedef std::piecewise_linear_distribution<> D;
195 typedef D::param_type P;
196 typedef std::mt19937_64 G;
197 G g;
198 double b[] = {10, 14, 16};
199 double p[] = {0, 1, 0};
200 const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
201 D d(b, b+Np+1, p);
202 const int N = 1000000;
203 std::vector<D::result_type> u;
204 for (int i = 0; i < N; ++i)
205 {
206 D::result_type v = d(g);
207 assert(d.min() <= v && v < d.max());
208 u.push_back(v);
209 }
210 std::sort(u.begin(), u.end());
211 int kp = -1;
212 double a;
213 double m;
214 double bk;
215 double c;
216 std::vector<double> areas(Np);
217 double S = 0;
218 for (int i = 0; i < areas.size(); ++i)
219 {
220 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
221 S += areas[i];
222 }
223 for (int i = 0; i < areas.size(); ++i)
224 areas[i] /= S;
225 for (int i = 0; i < Np+1; ++i)
226 p[i] /= S;
227 for (int i = 0; i < N; ++i)
228 {
229 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
230 if (k != kp)
231 {
232 a = 0;
233 for (int j = 0; j < k; ++j)
234 a += areas[j];
235 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
236 bk = b[k];
237 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
238 kp = k;
239 }
240 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
241 }
242 }
243 {
244 typedef std::piecewise_linear_distribution<> D;
245 typedef D::param_type P;
246 typedef std::mt19937_64 G;
247 G g;
248 double b[] = {10, 14};
249 double p[] = {1, 1};
250 const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
251 D d(b, b+Np+1, p);
252 const int N = 1000000;
253 std::vector<D::result_type> u;
254 for (int i = 0; i < N; ++i)
255 {
256 D::result_type v = d(g);
257 assert(d.min() <= v && v < d.max());
258 u.push_back(v);
259 }
260 std::sort(u.begin(), u.end());
261 int kp = -1;
262 double a;
263 double m;
264 double bk;
265 double c;
266 std::vector<double> areas(Np);
267 double S = 0;
268 for (int i = 0; i < areas.size(); ++i)
269 {
270 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
271 S += areas[i];
272 }
273 for (int i = 0; i < areas.size(); ++i)
274 areas[i] /= S;
275 for (int i = 0; i < Np+1; ++i)
276 p[i] /= S;
277 for (int i = 0; i < N; ++i)
278 {
279 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
280 if (k != kp)
281 {
282 a = 0;
283 for (int j = 0; j < k; ++j)
284 a += areas[j];
285 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
286 bk = b[k];
287 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
288 kp = k;
289 }
290 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
291 }
292 }
293 {
294 typedef std::piecewise_linear_distribution<> D;
295 typedef D::param_type P;
296 typedef std::mt19937_64 G;
297 G g;
298 double b[] = {10, 14, 16, 17};
299 double p[] = {25, 62.5, 12.5, 0};
300 const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
301 D d(b, b+Np+1, p);
302 const int N = 1000000;
303 std::vector<D::result_type> u;
304 for (int i = 0; i < N; ++i)
305 {
306 D::result_type v = d(g);
307 assert(d.min() <= v && v < d.max());
308 u.push_back(v);
309 }
310 std::sort(u.begin(), u.end());
311 int kp = -1;
312 double a;
313 double m;
314 double bk;
315 double c;
316 std::vector<double> areas(Np);
317 double S = 0;
318 for (int i = 0; i < areas.size(); ++i)
319 {
320 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
321 S += areas[i];
322 }
323 for (int i = 0; i < areas.size(); ++i)
324 areas[i] /= S;
325 for (int i = 0; i < Np+1; ++i)
326 p[i] /= S;
327 for (int i = 0; i < N; ++i)
328 {
329 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
330 if (k != kp)
331 {
332 a = 0;
333 for (int j = 0; j < k; ++j)
334 a += areas[j];
335 m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
336 bk = b[k];
337 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
338 kp = k;
339 }
340 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
341 }
342 }
343 }
344