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_constant_distribution
14
15// template<class _URNG> result_type operator()(_URNG& g);
16
17#include <random>
18#include <vector>
19#include <iterator>
20#include <numeric>
21#include <cassert>
22
23template <class T>
24inline
25T
26sqr(T x)
27{
28    return x*x;
29}
30
31int main()
32{
33    {
34        typedef std::piecewise_constant_distribution<> D;
35        typedef std::mt19937_64 G;
36        G g;
37        double b[] = {10, 14, 16, 17};
38        double p[] = {25, 62.5, 12.5};
39        const size_t Np = sizeof(p) / sizeof(p[0]);
40        D d(b, b+Np+1, p);
41        const int N = 1000000;
42        std::vector<D::result_type> u;
43        for (int i = 0; i < N; ++i)
44        {
45            D::result_type v = d(g);
46            assert(d.min() <= v && v < d.max());
47            u.push_back(v);
48        }
49        std::vector<double> prob(std::begin(p), std::end(p));
50        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
51        for (int i = 0; i < prob.size(); ++i)
52            prob[i] /= s;
53        std::sort(u.begin(), u.end());
54        for (int i = 0; i < Np; ++i)
55        {
56            typedef std::vector<D::result_type>::iterator I;
57            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
58            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
59            const size_t Ni = ub - lb;
60            if (prob[i] == 0)
61                assert(Ni == 0);
62            else
63            {
64                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
65                double mean = std::accumulate(lb, ub, 0.0) / Ni;
66                double var = 0;
67                double skew = 0;
68                double kurtosis = 0;
69                for (I j = lb; j != ub; ++j)
70                {
71                    double d = (*j - mean);
72                    double d2 = sqr(d);
73                    var += d2;
74                    skew += d * d2;
75                    kurtosis += d2 * d2;
76                }
77                var /= Ni;
78                double dev = std::sqrt(var);
79                skew /= Ni * dev * var;
80                kurtosis /= Ni * var * var;
81                kurtosis -= 3;
82                double x_mean = (b[i+1] + b[i]) / 2;
83                double x_var = sqr(b[i+1] - b[i]) / 12;
84                double x_skew = 0;
85                double x_kurtosis = -6./5;
86                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
87                assert(std::abs((var - x_var) / x_var) < 0.01);
88                assert(std::abs(skew - x_skew) < 0.01);
89                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
90            }
91        }
92    }
93    {
94        typedef std::piecewise_constant_distribution<> D;
95        typedef std::mt19937_64 G;
96        G g;
97        double b[] = {10, 14, 16, 17};
98        double p[] = {0, 62.5, 12.5};
99        const size_t Np = sizeof(p) / sizeof(p[0]);
100        D d(b, b+Np+1, p);
101        const int N = 1000000;
102        std::vector<D::result_type> u;
103        for (int i = 0; i < N; ++i)
104        {
105            D::result_type v = d(g);
106            assert(d.min() <= v && v < d.max());
107            u.push_back(v);
108        }
109        std::vector<double> prob(std::begin(p), std::end(p));
110        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
111        for (int i = 0; i < prob.size(); ++i)
112            prob[i] /= s;
113        std::sort(u.begin(), u.end());
114        for (int i = 0; i < Np; ++i)
115        {
116            typedef std::vector<D::result_type>::iterator I;
117            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
118            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
119            const size_t Ni = ub - lb;
120            if (prob[i] == 0)
121                assert(Ni == 0);
122            else
123            {
124                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
125                double mean = std::accumulate(lb, ub, 0.0) / Ni;
126                double var = 0;
127                double skew = 0;
128                double kurtosis = 0;
129                for (I j = lb; j != ub; ++j)
130                {
131                    double d = (*j - mean);
132                    double d2 = sqr(d);
133                    var += d2;
134                    skew += d * d2;
135                    kurtosis += d2 * d2;
136                }
137                var /= Ni;
138                double dev = std::sqrt(var);
139                skew /= Ni * dev * var;
140                kurtosis /= Ni * var * var;
141                kurtosis -= 3;
142                double x_mean = (b[i+1] + b[i]) / 2;
143                double x_var = sqr(b[i+1] - b[i]) / 12;
144                double x_skew = 0;
145                double x_kurtosis = -6./5;
146                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
147                assert(std::abs((var - x_var) / x_var) < 0.01);
148                assert(std::abs(skew - x_skew) < 0.01);
149                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
150            }
151        }
152    }
153    {
154        typedef std::piecewise_constant_distribution<> D;
155        typedef std::mt19937_64 G;
156        G g;
157        double b[] = {10, 14, 16, 17};
158        double p[] = {25, 0, 12.5};
159        const size_t Np = sizeof(p) / sizeof(p[0]);
160        D d(b, b+Np+1, p);
161        const int N = 1000000;
162        std::vector<D::result_type> u;
163        for (int i = 0; i < N; ++i)
164        {
165            D::result_type v = d(g);
166            assert(d.min() <= v && v < d.max());
167            u.push_back(v);
168        }
169        std::vector<double> prob(std::begin(p), std::end(p));
170        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
171        for (int i = 0; i < prob.size(); ++i)
172            prob[i] /= s;
173        std::sort(u.begin(), u.end());
174        for (int i = 0; i < Np; ++i)
175        {
176            typedef std::vector<D::result_type>::iterator I;
177            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
178            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
179            const size_t Ni = ub - lb;
180            if (prob[i] == 0)
181                assert(Ni == 0);
182            else
183            {
184                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
185                double mean = std::accumulate(lb, ub, 0.0) / Ni;
186                double var = 0;
187                double skew = 0;
188                double kurtosis = 0;
189                for (I j = lb; j != ub; ++j)
190                {
191                    double d = (*j - mean);
192                    double d2 = sqr(d);
193                    var += d2;
194                    skew += d * d2;
195                    kurtosis += d2 * d2;
196                }
197                var /= Ni;
198                double dev = std::sqrt(var);
199                skew /= Ni * dev * var;
200                kurtosis /= Ni * var * var;
201                kurtosis -= 3;
202                double x_mean = (b[i+1] + b[i]) / 2;
203                double x_var = sqr(b[i+1] - b[i]) / 12;
204                double x_skew = 0;
205                double x_kurtosis = -6./5;
206                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
207                assert(std::abs((var - x_var) / x_var) < 0.01);
208                assert(std::abs(skew - x_skew) < 0.01);
209                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
210            }
211        }
212    }
213    {
214        typedef std::piecewise_constant_distribution<> D;
215        typedef std::mt19937_64 G;
216        G g;
217        double b[] = {10, 14, 16, 17};
218        double p[] = {25, 62.5, 0};
219        const size_t Np = sizeof(p) / sizeof(p[0]);
220        D d(b, b+Np+1, p);
221        const int N = 1000000;
222        std::vector<D::result_type> u;
223        for (int i = 0; i < N; ++i)
224        {
225            D::result_type v = d(g);
226            assert(d.min() <= v && v < d.max());
227            u.push_back(v);
228        }
229        std::vector<double> prob(std::begin(p), std::end(p));
230        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
231        for (int i = 0; i < prob.size(); ++i)
232            prob[i] /= s;
233        std::sort(u.begin(), u.end());
234        for (int i = 0; i < Np; ++i)
235        {
236            typedef std::vector<D::result_type>::iterator I;
237            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
238            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
239            const size_t Ni = ub - lb;
240            if (prob[i] == 0)
241                assert(Ni == 0);
242            else
243            {
244                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
245                double mean = std::accumulate(lb, ub, 0.0) / Ni;
246                double var = 0;
247                double skew = 0;
248                double kurtosis = 0;
249                for (I j = lb; j != ub; ++j)
250                {
251                    double d = (*j - mean);
252                    double d2 = sqr(d);
253                    var += d2;
254                    skew += d * d2;
255                    kurtosis += d2 * d2;
256                }
257                var /= Ni;
258                double dev = std::sqrt(var);
259                skew /= Ni * dev * var;
260                kurtosis /= Ni * var * var;
261                kurtosis -= 3;
262                double x_mean = (b[i+1] + b[i]) / 2;
263                double x_var = sqr(b[i+1] - b[i]) / 12;
264                double x_skew = 0;
265                double x_kurtosis = -6./5;
266                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
267                assert(std::abs((var - x_var) / x_var) < 0.01);
268                assert(std::abs(skew - x_skew) < 0.01);
269                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
270            }
271        }
272    }
273    {
274        typedef std::piecewise_constant_distribution<> D;
275        typedef std::mt19937_64 G;
276        G g;
277        double b[] = {10, 14, 16, 17};
278        double p[] = {25, 0, 0};
279        const size_t Np = sizeof(p) / sizeof(p[0]);
280        D d(b, b+Np+1, p);
281        const int N = 100000;
282        std::vector<D::result_type> u;
283        for (int i = 0; i < N; ++i)
284        {
285            D::result_type v = d(g);
286            assert(d.min() <= v && v < d.max());
287            u.push_back(v);
288        }
289        std::vector<double> prob(std::begin(p), std::end(p));
290        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
291        for (int i = 0; i < prob.size(); ++i)
292            prob[i] /= s;
293        std::sort(u.begin(), u.end());
294        for (int i = 0; i < Np; ++i)
295        {
296            typedef std::vector<D::result_type>::iterator I;
297            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
298            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
299            const size_t Ni = ub - lb;
300            if (prob[i] == 0)
301                assert(Ni == 0);
302            else
303            {
304                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
305                double mean = std::accumulate(lb, ub, 0.0) / Ni;
306                double var = 0;
307                double skew = 0;
308                double kurtosis = 0;
309                for (I j = lb; j != ub; ++j)
310                {
311                    double d = (*j - mean);
312                    double d2 = sqr(d);
313                    var += d2;
314                    skew += d * d2;
315                    kurtosis += d2 * d2;
316                }
317                var /= Ni;
318                double dev = std::sqrt(var);
319                skew /= Ni * dev * var;
320                kurtosis /= Ni * var * var;
321                kurtosis -= 3;
322                double x_mean = (b[i+1] + b[i]) / 2;
323                double x_var = sqr(b[i+1] - b[i]) / 12;
324                double x_skew = 0;
325                double x_kurtosis = -6./5;
326                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
327                assert(std::abs((var - x_var) / x_var) < 0.01);
328                assert(std::abs(skew - x_skew) < 0.01);
329                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
330            }
331        }
332    }
333    {
334        typedef std::piecewise_constant_distribution<> D;
335        typedef std::mt19937_64 G;
336        G g;
337        double b[] = {10, 14, 16, 17};
338        double p[] = {0, 25, 0};
339        const size_t Np = sizeof(p) / sizeof(p[0]);
340        D d(b, b+Np+1, p);
341        const int N = 100000;
342        std::vector<D::result_type> u;
343        for (int i = 0; i < N; ++i)
344        {
345            D::result_type v = d(g);
346            assert(d.min() <= v && v < d.max());
347            u.push_back(v);
348        }
349        std::vector<double> prob(std::begin(p), std::end(p));
350        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
351        for (int i = 0; i < prob.size(); ++i)
352            prob[i] /= s;
353        std::sort(u.begin(), u.end());
354        for (int i = 0; i < Np; ++i)
355        {
356            typedef std::vector<D::result_type>::iterator I;
357            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
358            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
359            const size_t Ni = ub - lb;
360            if (prob[i] == 0)
361                assert(Ni == 0);
362            else
363            {
364                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
365                double mean = std::accumulate(lb, ub, 0.0) / Ni;
366                double var = 0;
367                double skew = 0;
368                double kurtosis = 0;
369                for (I j = lb; j != ub; ++j)
370                {
371                    double d = (*j - mean);
372                    double d2 = sqr(d);
373                    var += d2;
374                    skew += d * d2;
375                    kurtosis += d2 * d2;
376                }
377                var /= Ni;
378                double dev = std::sqrt(var);
379                skew /= Ni * dev * var;
380                kurtosis /= Ni * var * var;
381                kurtosis -= 3;
382                double x_mean = (b[i+1] + b[i]) / 2;
383                double x_var = sqr(b[i+1] - b[i]) / 12;
384                double x_skew = 0;
385                double x_kurtosis = -6./5;
386                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
387                assert(std::abs((var - x_var) / x_var) < 0.01);
388                assert(std::abs(skew - x_skew) < 0.01);
389                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
390            }
391        }
392    }
393    {
394        typedef std::piecewise_constant_distribution<> D;
395        typedef std::mt19937_64 G;
396        G g;
397        double b[] = {10, 14, 16, 17};
398        double p[] = {0, 0, 1};
399        const size_t Np = sizeof(p) / sizeof(p[0]);
400        D d(b, b+Np+1, p);
401        const int N = 100000;
402        std::vector<D::result_type> u;
403        for (int i = 0; i < N; ++i)
404        {
405            D::result_type v = d(g);
406            assert(d.min() <= v && v < d.max());
407            u.push_back(v);
408        }
409        std::vector<double> prob(std::begin(p), std::end(p));
410        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
411        for (int i = 0; i < prob.size(); ++i)
412            prob[i] /= s;
413        std::sort(u.begin(), u.end());
414        for (int i = 0; i < Np; ++i)
415        {
416            typedef std::vector<D::result_type>::iterator I;
417            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
418            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
419            const size_t Ni = ub - lb;
420            if (prob[i] == 0)
421                assert(Ni == 0);
422            else
423            {
424                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
425                double mean = std::accumulate(lb, ub, 0.0) / Ni;
426                double var = 0;
427                double skew = 0;
428                double kurtosis = 0;
429                for (I j = lb; j != ub; ++j)
430                {
431                    double d = (*j - mean);
432                    double d2 = sqr(d);
433                    var += d2;
434                    skew += d * d2;
435                    kurtosis += d2 * d2;
436                }
437                var /= Ni;
438                double dev = std::sqrt(var);
439                skew /= Ni * dev * var;
440                kurtosis /= Ni * var * var;
441                kurtosis -= 3;
442                double x_mean = (b[i+1] + b[i]) / 2;
443                double x_var = sqr(b[i+1] - b[i]) / 12;
444                double x_skew = 0;
445                double x_kurtosis = -6./5;
446                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
447                assert(std::abs((var - x_var) / x_var) < 0.01);
448                assert(std::abs(skew - x_skew) < 0.01);
449                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
450            }
451        }
452    }
453    {
454        typedef std::piecewise_constant_distribution<> D;
455        typedef std::mt19937_64 G;
456        G g;
457        double b[] = {10, 14, 16};
458        double p[] = {75, 25};
459        const size_t Np = sizeof(p) / sizeof(p[0]);
460        D d(b, b+Np+1, p);
461        const int N = 100000;
462        std::vector<D::result_type> u;
463        for (int i = 0; i < N; ++i)
464        {
465            D::result_type v = d(g);
466            assert(d.min() <= v && v < d.max());
467            u.push_back(v);
468        }
469        std::vector<double> prob(std::begin(p), std::end(p));
470        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
471        for (int i = 0; i < prob.size(); ++i)
472            prob[i] /= s;
473        std::sort(u.begin(), u.end());
474        for (int i = 0; i < Np; ++i)
475        {
476            typedef std::vector<D::result_type>::iterator I;
477            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
478            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
479            const size_t Ni = ub - lb;
480            if (prob[i] == 0)
481                assert(Ni == 0);
482            else
483            {
484                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
485                double mean = std::accumulate(lb, ub, 0.0) / Ni;
486                double var = 0;
487                double skew = 0;
488                double kurtosis = 0;
489                for (I j = lb; j != ub; ++j)
490                {
491                    double d = (*j - mean);
492                    double d2 = sqr(d);
493                    var += d2;
494                    skew += d * d2;
495                    kurtosis += d2 * d2;
496                }
497                var /= Ni;
498                double dev = std::sqrt(var);
499                skew /= Ni * dev * var;
500                kurtosis /= Ni * var * var;
501                kurtosis -= 3;
502                double x_mean = (b[i+1] + b[i]) / 2;
503                double x_var = sqr(b[i+1] - b[i]) / 12;
504                double x_skew = 0;
505                double x_kurtosis = -6./5;
506                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
507                assert(std::abs((var - x_var) / x_var) < 0.01);
508                assert(std::abs(skew - x_skew) < 0.01);
509                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
510            }
511        }
512    }
513    {
514        typedef std::piecewise_constant_distribution<> D;
515        typedef std::mt19937_64 G;
516        G g;
517        double b[] = {10, 14, 16};
518        double p[] = {0, 25};
519        const size_t Np = sizeof(p) / sizeof(p[0]);
520        D d(b, b+Np+1, p);
521        const int N = 100000;
522        std::vector<D::result_type> u;
523        for (int i = 0; i < N; ++i)
524        {
525            D::result_type v = d(g);
526            assert(d.min() <= v && v < d.max());
527            u.push_back(v);
528        }
529        std::vector<double> prob(std::begin(p), std::end(p));
530        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
531        for (int i = 0; i < prob.size(); ++i)
532            prob[i] /= s;
533        std::sort(u.begin(), u.end());
534        for (int i = 0; i < Np; ++i)
535        {
536            typedef std::vector<D::result_type>::iterator I;
537            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
538            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
539            const size_t Ni = ub - lb;
540            if (prob[i] == 0)
541                assert(Ni == 0);
542            else
543            {
544                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
545                double mean = std::accumulate(lb, ub, 0.0) / Ni;
546                double var = 0;
547                double skew = 0;
548                double kurtosis = 0;
549                for (I j = lb; j != ub; ++j)
550                {
551                    double d = (*j - mean);
552                    double d2 = sqr(d);
553                    var += d2;
554                    skew += d * d2;
555                    kurtosis += d2 * d2;
556                }
557                var /= Ni;
558                double dev = std::sqrt(var);
559                skew /= Ni * dev * var;
560                kurtosis /= Ni * var * var;
561                kurtosis -= 3;
562                double x_mean = (b[i+1] + b[i]) / 2;
563                double x_var = sqr(b[i+1] - b[i]) / 12;
564                double x_skew = 0;
565                double x_kurtosis = -6./5;
566                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
567                assert(std::abs((var - x_var) / x_var) < 0.01);
568                assert(std::abs(skew - x_skew) < 0.01);
569                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
570            }
571        }
572    }
573    {
574        typedef std::piecewise_constant_distribution<> D;
575        typedef std::mt19937_64 G;
576        G g;
577        double b[] = {10, 14, 16};
578        double p[] = {1, 0};
579        const size_t Np = sizeof(p) / sizeof(p[0]);
580        D d(b, b+Np+1, p);
581        const int N = 100000;
582        std::vector<D::result_type> u;
583        for (int i = 0; i < N; ++i)
584        {
585            D::result_type v = d(g);
586            assert(d.min() <= v && v < d.max());
587            u.push_back(v);
588        }
589        std::vector<double> prob(std::begin(p), std::end(p));
590        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
591        for (int i = 0; i < prob.size(); ++i)
592            prob[i] /= s;
593        std::sort(u.begin(), u.end());
594        for (int i = 0; i < Np; ++i)
595        {
596            typedef std::vector<D::result_type>::iterator I;
597            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
598            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
599            const size_t Ni = ub - lb;
600            if (prob[i] == 0)
601                assert(Ni == 0);
602            else
603            {
604                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
605                double mean = std::accumulate(lb, ub, 0.0) / Ni;
606                double var = 0;
607                double skew = 0;
608                double kurtosis = 0;
609                for (I j = lb; j != ub; ++j)
610                {
611                    double d = (*j - mean);
612                    double d2 = sqr(d);
613                    var += d2;
614                    skew += d * d2;
615                    kurtosis += d2 * d2;
616                }
617                var /= Ni;
618                double dev = std::sqrt(var);
619                skew /= Ni * dev * var;
620                kurtosis /= Ni * var * var;
621                kurtosis -= 3;
622                double x_mean = (b[i+1] + b[i]) / 2;
623                double x_var = sqr(b[i+1] - b[i]) / 12;
624                double x_skew = 0;
625                double x_kurtosis = -6./5;
626                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
627                assert(std::abs((var - x_var) / x_var) < 0.01);
628                assert(std::abs(skew - x_skew) < 0.01);
629                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
630            }
631        }
632    }
633    {
634        typedef std::piecewise_constant_distribution<> D;
635        typedef std::mt19937_64 G;
636        G g;
637        double b[] = {10, 14};
638        double p[] = {1};
639        const size_t Np = sizeof(p) / sizeof(p[0]);
640        D d(b, b+Np+1, p);
641        const int N = 100000;
642        std::vector<D::result_type> u;
643        for (int i = 0; i < N; ++i)
644        {
645            D::result_type v = d(g);
646            assert(d.min() <= v && v < d.max());
647            u.push_back(v);
648        }
649        std::vector<double> prob(std::begin(p), std::end(p));
650        double s = std::accumulate(prob.begin(), prob.end(), 0.0);
651        for (int i = 0; i < prob.size(); ++i)
652            prob[i] /= s;
653        std::sort(u.begin(), u.end());
654        for (int i = 0; i < Np; ++i)
655        {
656            typedef std::vector<D::result_type>::iterator I;
657            I lb = std::lower_bound(u.begin(), u.end(), b[i]);
658            I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
659            const size_t Ni = ub - lb;
660            if (prob[i] == 0)
661                assert(Ni == 0);
662            else
663            {
664                assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
665                double mean = std::accumulate(lb, ub, 0.0) / Ni;
666                double var = 0;
667                double skew = 0;
668                double kurtosis = 0;
669                for (I j = lb; j != ub; ++j)
670                {
671                    double d = (*j - mean);
672                    double d2 = sqr(d);
673                    var += d2;
674                    skew += d * d2;
675                    kurtosis += d2 * d2;
676                }
677                var /= Ni;
678                double dev = std::sqrt(var);
679                skew /= Ni * dev * var;
680                kurtosis /= Ni * var * var;
681                kurtosis -= 3;
682                double x_mean = (b[i+1] + b[i]) / 2;
683                double x_var = sqr(b[i+1] - b[i]) / 12;
684                double x_skew = 0;
685                double x_kurtosis = -6./5;
686                assert(std::abs((mean - x_mean) / x_mean) < 0.01);
687                assert(std::abs((var - x_var) / x_var) < 0.01);
688                assert(std::abs(skew - x_skew) < 0.01);
689                assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
690            }
691        }
692    }
693}
694