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-exception 6 // 7 //===----------------------------------------------------------------------===// 8 // 9 // REQUIRES: long_tests 10 11 // <random> 12 13 // template<class RealType = double> 14 // class piecewise_linear_distribution 15 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 28 template <class T> 29 inline 30 T 31 sqr(T x) 32 { 33 return x*x; 34 } 35 36 double 37 f(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 42 void 43 test1() 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 93 void 94 test2() 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 144 void 145 test3() 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 195 void 196 test4() 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 247 void 248 test5() 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 300 void 301 test6() 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 351 int main(int, char**) 352 { 353 test1(); 354 test2(); 355 test3(); 356 test4(); 357 test5(); 358 test6(); 359 360 return 0; 361 } 362