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