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