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 <vector> 23 #include <iterator> 24 #include <numeric> 25 #include <cassert> 26 27 template <class T> 28 inline 29 T 30 sqr(T x) 31 { 32 return x*x; 33 } 34 35 double 36 f(double x, double a, double m, double b, double c) 37 { 38 return a + m*(sqr(x) - sqr(b))/2 + c*(x-b); 39 } 40 41 int main() 42 { 43 { 44 typedef std::piecewise_linear_distribution<> D; 45 typedef D::param_type P; 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 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 (int 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 int kp = -1; 62 double a; 63 double m; 64 double bk; 65 double c; 66 std::vector<double> areas(Np); 67 double S = 0; 68 for (int 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 (int i = 0; i < areas.size(); ++i) 74 areas[i] /= S; 75 for (int i = 0; i < Np+1; ++i) 76 p[i] /= S; 77 for (int i = 0; i < N; ++i) 78 { 79 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 80 if (k != kp) 81 { 82 a = 0; 83 for (int j = 0; j < k; ++j) 84 a += areas[j]; 85 m = (p[k+1] - p[k]) / (b[k+1] - b[k]); 86 bk = b[k]; 87 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]); 88 kp = k; 89 } 90 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001); 91 } 92 } 93 { 94 typedef std::piecewise_linear_distribution<> D; 95 typedef D::param_type P; 96 typedef std::mt19937_64 G; 97 G g; 98 double b[] = {10, 14, 16, 17}; 99 double p[] = {0, 0, 1, 0}; 100 const size_t Np = sizeof(p) / sizeof(p[0]) - 1; 101 D d(b, b+Np+1, p); 102 const int N = 1000000; 103 std::vector<D::result_type> u; 104 for (int i = 0; i < N; ++i) 105 { 106 D::result_type v = d(g); 107 assert(d.min() <= v && v < d.max()); 108 u.push_back(v); 109 } 110 std::sort(u.begin(), u.end()); 111 int kp = -1; 112 double a; 113 double m; 114 double bk; 115 double c; 116 std::vector<double> areas(Np); 117 double S = 0; 118 for (int i = 0; i < areas.size(); ++i) 119 { 120 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2; 121 S += areas[i]; 122 } 123 for (int i = 0; i < areas.size(); ++i) 124 areas[i] /= S; 125 for (int i = 0; i < Np+1; ++i) 126 p[i] /= S; 127 for (int i = 0; i < N; ++i) 128 { 129 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 130 if (k != kp) 131 { 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) < .001); 141 } 142 } 143 { 144 typedef std::piecewise_linear_distribution<> D; 145 typedef D::param_type P; 146 typedef std::mt19937_64 G; 147 G g; 148 double b[] = {10, 14, 16, 17}; 149 double p[] = {1, 0, 0, 0}; 150 const size_t Np = sizeof(p) / sizeof(p[0]) - 1; 151 D d(b, b+Np+1, p); 152 const int N = 1000000; 153 std::vector<D::result_type> u; 154 for (int i = 0; i < N; ++i) 155 { 156 D::result_type v = d(g); 157 assert(d.min() <= v && v < d.max()); 158 u.push_back(v); 159 } 160 std::sort(u.begin(), u.end()); 161 int kp = -1; 162 double a; 163 double m; 164 double bk; 165 double c; 166 std::vector<double> areas(Np); 167 double S = 0; 168 for (int i = 0; i < areas.size(); ++i) 169 { 170 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2; 171 S += areas[i]; 172 } 173 for (int i = 0; i < areas.size(); ++i) 174 areas[i] /= S; 175 for (int i = 0; i < Np+1; ++i) 176 p[i] /= S; 177 for (int i = 0; i < N; ++i) 178 { 179 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 180 if (k != kp) 181 { 182 a = 0; 183 for (int j = 0; j < k; ++j) 184 a += areas[j]; 185 m = (p[k+1] - p[k]) / (b[k+1] - b[k]); 186 bk = b[k]; 187 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]); 188 kp = k; 189 } 190 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001); 191 } 192 } 193 { 194 typedef std::piecewise_linear_distribution<> D; 195 typedef D::param_type P; 196 typedef std::mt19937_64 G; 197 G g; 198 double b[] = {10, 14, 16}; 199 double p[] = {0, 1, 0}; 200 const size_t Np = sizeof(p) / sizeof(p[0]) - 1; 201 D d(b, b+Np+1, p); 202 const int N = 1000000; 203 std::vector<D::result_type> u; 204 for (int i = 0; i < N; ++i) 205 { 206 D::result_type v = d(g); 207 assert(d.min() <= v && v < d.max()); 208 u.push_back(v); 209 } 210 std::sort(u.begin(), u.end()); 211 int kp = -1; 212 double a; 213 double m; 214 double bk; 215 double c; 216 std::vector<double> areas(Np); 217 double S = 0; 218 for (int i = 0; i < areas.size(); ++i) 219 { 220 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2; 221 S += areas[i]; 222 } 223 for (int i = 0; i < areas.size(); ++i) 224 areas[i] /= S; 225 for (int i = 0; i < Np+1; ++i) 226 p[i] /= S; 227 for (int i = 0; i < N; ++i) 228 { 229 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 230 if (k != kp) 231 { 232 a = 0; 233 for (int j = 0; j < k; ++j) 234 a += areas[j]; 235 m = (p[k+1] - p[k]) / (b[k+1] - b[k]); 236 bk = b[k]; 237 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]); 238 kp = k; 239 } 240 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001); 241 } 242 } 243 { 244 typedef std::piecewise_linear_distribution<> D; 245 typedef D::param_type P; 246 typedef std::mt19937_64 G; 247 G g; 248 double b[] = {10, 14}; 249 double p[] = {1, 1}; 250 const size_t Np = sizeof(p) / sizeof(p[0]) - 1; 251 D d(b, b+Np+1, p); 252 const int N = 1000000; 253 std::vector<D::result_type> u; 254 for (int i = 0; i < N; ++i) 255 { 256 D::result_type v = d(g); 257 assert(d.min() <= v && v < d.max()); 258 u.push_back(v); 259 } 260 std::sort(u.begin(), u.end()); 261 int kp = -1; 262 double a; 263 double m; 264 double bk; 265 double c; 266 std::vector<double> areas(Np); 267 double S = 0; 268 for (int i = 0; i < areas.size(); ++i) 269 { 270 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2; 271 S += areas[i]; 272 } 273 for (int i = 0; i < areas.size(); ++i) 274 areas[i] /= S; 275 for (int i = 0; i < Np+1; ++i) 276 p[i] /= S; 277 for (int i = 0; i < N; ++i) 278 { 279 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 280 if (k != kp) 281 { 282 a = 0; 283 for (int j = 0; j < k; ++j) 284 a += areas[j]; 285 m = (p[k+1] - p[k]) / (b[k+1] - b[k]); 286 bk = b[k]; 287 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]); 288 kp = k; 289 } 290 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001); 291 } 292 } 293 { 294 typedef std::piecewise_linear_distribution<> D; 295 typedef D::param_type P; 296 typedef std::mt19937_64 G; 297 G g; 298 double b[] = {10, 14, 16, 17}; 299 double p[] = {25, 62.5, 12.5, 0}; 300 const size_t Np = sizeof(p) / sizeof(p[0]) - 1; 301 D d(b, b+Np+1, p); 302 const int N = 1000000; 303 std::vector<D::result_type> u; 304 for (int i = 0; i < N; ++i) 305 { 306 D::result_type v = d(g); 307 assert(d.min() <= v && v < d.max()); 308 u.push_back(v); 309 } 310 std::sort(u.begin(), u.end()); 311 int kp = -1; 312 double a; 313 double m; 314 double bk; 315 double c; 316 std::vector<double> areas(Np); 317 double S = 0; 318 for (int i = 0; i < areas.size(); ++i) 319 { 320 areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2; 321 S += areas[i]; 322 } 323 for (int i = 0; i < areas.size(); ++i) 324 areas[i] /= S; 325 for (int i = 0; i < Np+1; ++i) 326 p[i] /= S; 327 for (int i = 0; i < N; ++i) 328 { 329 int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1; 330 if (k != kp) 331 { 332 a = 0; 333 for (int j = 0; j < k; ++j) 334 a += areas[j]; 335 m = (p[k+1] - p[k]) / (b[k+1] - b[k]); 336 bk = b[k]; 337 c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]); 338 kp = k; 339 } 340 assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001); 341 } 342 } 343 } 344
