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
29template <class T>
30inline
31T
32sqr(T x)
33{
34    return x*x;
35}
36
37double
38f(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
43void
44test1()
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
95void
96test2()
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
147void
148test3()
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
199void
200test4()
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
252void
253test5()
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
306void
307test6()
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
358int main()
359{
360    test1();
361    test2();
362    test3();
363    test4();
364    test5();
365    test6();
366}
367