362 lines · cpp
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8//9// REQUIRES: long_tests10 11// <random>12 13// template<class RealType = double>14// class piecewise_linear_distribution15 16// template<class _URNG> result_type operator()(_URNG& g);17 18#include <algorithm>19#include <cassert>20#include <cmath>21#include <cstddef>22#include <limits>23#include <random>24#include <vector>25 26#include "test_macros.h"27 28template <class T>29inline30T31sqr(T x)32{33 return x*x;34}35 36double37f(double x, double a, double m, double b, double c)38{39 return a + m*(sqr(x) - sqr(b))/2 + c*(x-b);40}41 42void43test1()44{45 typedef std::piecewise_linear_distribution<> D;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 std::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 (std::size_t 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 std::ptrdiff_t kp = -1;62 double a = std::numeric_limits<double>::quiet_NaN();63 double m = std::numeric_limits<double>::quiet_NaN();64 double bk = std::numeric_limits<double>::quiet_NaN();65 double c = std::numeric_limits<double>::quiet_NaN();66 std::vector<double> areas(Np);67 double S = 0;68 for (std::size_t 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 (std::size_t i = 0; i < areas.size(); ++i)74 areas[i] /= S;75 for (std::size_t i = 0; i < Np+1; ++i)76 p[i] /= S;77 for (std::size_t i = 0; i < N; ++i)78 {79 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;80 if (k != kp) {81 a = 0;82 for (int j = 0; j < k; ++j)83 a += areas[j];84 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);85 bk = b[k];86 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);87 kp = k;88 }89 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);90 }91}92 93void94test2()95{96 typedef std::piecewise_linear_distribution<> D;97 typedef std::mt19937_64 G;98 G g;99 double b[] = {10, 14, 16, 17};100 double p[] = {0, 0, 1, 0};101 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;102 D d(b, b+Np+1, p);103 const int N = 1000000;104 std::vector<D::result_type> u;105 for (std::size_t i = 0; i < N; ++i)106 {107 D::result_type v = d(g);108 assert(d.min() <= v && v < d.max());109 u.push_back(v);110 }111 std::sort(u.begin(), u.end());112 std::ptrdiff_t kp = -1;113 double a = std::numeric_limits<double>::quiet_NaN();114 double m = std::numeric_limits<double>::quiet_NaN();115 double bk = std::numeric_limits<double>::quiet_NaN();116 double c = std::numeric_limits<double>::quiet_NaN();117 std::vector<double> areas(Np);118 double S = 0;119 for (std::size_t i = 0; i < areas.size(); ++i)120 {121 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;122 S += areas[i];123 }124 for (std::size_t i = 0; i < areas.size(); ++i)125 areas[i] /= S;126 for (std::size_t i = 0; i < Np+1; ++i)127 p[i] /= S;128 for (std::size_t i = 0; i < N; ++i)129 {130 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;131 if (k != kp) {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) < .0013);141 }142}143 144void145test3()146{147 typedef std::piecewise_linear_distribution<> D;148 typedef std::mt19937_64 G;149 G g;150 double b[] = {10, 14, 16, 17};151 double p[] = {1, 0, 0, 0};152 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;153 D d(b, b+Np+1, p);154 const std::size_t N = 1000000;155 std::vector<D::result_type> u;156 for (std::size_t i = 0; i < N; ++i)157 {158 D::result_type v = d(g);159 assert(d.min() <= v && v < d.max());160 u.push_back(v);161 }162 std::sort(u.begin(), u.end());163 std::ptrdiff_t kp = -1;164 double a = std::numeric_limits<double>::quiet_NaN();165 double m = std::numeric_limits<double>::quiet_NaN();166 double bk = std::numeric_limits<double>::quiet_NaN();167 double c = std::numeric_limits<double>::quiet_NaN();168 std::vector<double> areas(Np);169 double S = 0;170 for (std::size_t i = 0; i < areas.size(); ++i)171 {172 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;173 S += areas[i];174 }175 for (std::size_t i = 0; i < areas.size(); ++i)176 areas[i] /= S;177 for (std::size_t i = 0; i < Np+1; ++i)178 p[i] /= S;179 for (std::size_t i = 0; i < N; ++i)180 {181 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;182 if (k != kp) {183 a = 0;184 for (int j = 0; j < k; ++j)185 a += areas[j];186 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);187 bk = b[k];188 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);189 kp = k;190 }191 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);192 }193}194 195void196test4()197{198 typedef std::piecewise_linear_distribution<> D;199 typedef std::mt19937_64 G;200 G g;201 double b[] = {10, 14, 16};202 double p[] = {0, 1, 0};203 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;204 D d(b, b+Np+1, p);205 const int N = 1000000;206 std::vector<D::result_type> u;207 for (std::size_t i = 0; i < N; ++i)208 {209 D::result_type v = d(g);210 assert(d.min() <= v && v < d.max());211 u.push_back(v);212 }213 std::sort(u.begin(), u.end());214 std::ptrdiff_t kp = -1;215 double a = std::numeric_limits<double>::quiet_NaN();216 double m = std::numeric_limits<double>::quiet_NaN();217 double bk = std::numeric_limits<double>::quiet_NaN();218 double c = std::numeric_limits<double>::quiet_NaN();219 std::vector<double> areas(Np);220 double S = 0;221 for (std::size_t i = 0; i < areas.size(); ++i)222 {223 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;224 S += areas[i];225 }226 for (std::size_t i = 0; i < areas.size(); ++i)227 areas[i] /= S;228 for (std::size_t i = 0; i < Np+1; ++i)229 p[i] /= S;230 for (std::size_t i = 0; i < N; ++i)231 {232 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;233 if (k != kp) {234 a = 0;235 for (int j = 0; j < k; ++j)236 a += areas[j];237 assert(k < static_cast<int>(Np));238 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);239 bk = b[k];240 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);241 kp = k;242 }243 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);244 }245}246 247void248test5()249{250 typedef std::piecewise_linear_distribution<> D;251 typedef std::mt19937_64 G;252 G g;253 double b[] = {10, 14};254 double p[] = {1, 1};255 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;256 D d(b, b+Np+1, p);257 const int N = 1000000;258 std::vector<D::result_type> u;259 for (std::size_t i = 0; i < N; ++i)260 {261 D::result_type v = d(g);262 assert(d.min() <= v && v < d.max());263 u.push_back(v);264 }265 std::sort(u.begin(), u.end());266 std::ptrdiff_t kp = -1;267 double a = std::numeric_limits<double>::quiet_NaN();268 double m = std::numeric_limits<double>::quiet_NaN();269 double bk = std::numeric_limits<double>::quiet_NaN();270 double c = std::numeric_limits<double>::quiet_NaN();271 std::vector<double> areas(Np);272 double S = 0;273 for (std::size_t i = 0; i < areas.size(); ++i)274 {275 assert(i < Np);276 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;277 S += areas[i];278 }279 for (std::size_t i = 0; i < areas.size(); ++i)280 areas[i] /= S;281 for (std::size_t i = 0; i < Np+1; ++i)282 p[i] /= S;283 for (std::size_t i = 0; i < N; ++i)284 {285 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;286 if (k != kp) {287 a = 0;288 for (int j = 0; j < k; ++j)289 a += areas[j];290 assert(k < static_cast<int>(Np));291 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);292 bk = b[k];293 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);294 kp = k;295 }296 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);297 }298}299 300void301test6()302{303 typedef std::piecewise_linear_distribution<> D;304 typedef std::mt19937_64 G;305 G g;306 double b[] = {10, 14, 16, 17};307 double p[] = {25, 62.5, 12.5, 0};308 const std::size_t Np = sizeof(p) / sizeof(p[0]) - 1;309 D d(b, b+Np+1, p);310 const int N = 1000000;311 std::vector<D::result_type> u;312 for (std::size_t i = 0; i < N; ++i)313 {314 D::result_type v = d(g);315 assert(d.min() <= v && v < d.max());316 u.push_back(v);317 }318 std::sort(u.begin(), u.end());319 std::ptrdiff_t kp = -1;320 double a = std::numeric_limits<double>::quiet_NaN();321 double m = std::numeric_limits<double>::quiet_NaN();322 double bk = std::numeric_limits<double>::quiet_NaN();323 double c = std::numeric_limits<double>::quiet_NaN();324 std::vector<double> areas(Np);325 double S = 0;326 for (std::size_t i = 0; i < areas.size(); ++i)327 {328 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;329 S += areas[i];330 }331 for (std::size_t i = 0; i < areas.size(); ++i)332 areas[i] /= S;333 for (std::size_t i = 0; i < Np+1; ++i)334 p[i] /= S;335 for (std::size_t i = 0; i < N; ++i)336 {337 std::ptrdiff_t k = std::lower_bound(b, b + Np + 1, u[i]) - b - 1;338 if (k != kp) {339 a = 0;340 for (int j = 0; j < k; ++j)341 a += areas[j];342 m = (p[k + 1] - p[k]) / (b[k + 1] - b[k]);343 bk = b[k];344 c = (b[k + 1] * p[k] - b[k] * p[k + 1]) / (b[k + 1] - b[k]);345 kp = k;346 }347 assert(std::abs(f(u[i], a, m, bk, c) - double(i) / N) < .0013);348 }349}350 351int main(int, char**)352{353 test1();354 test2();355 test3();356 test4();357 test5();358 test6();359 360 return 0;361}362