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