1// Random number extensions -*- C++ -*- 2 3// Copyright (C) 2012-2014 Free Software Foundation, Inc. 4// 5// This file is part of the GNU ISO C++ Library. This library is free 6// software; you can redistribute it and/or modify it under the 7// terms of the GNU General Public License as published by the 8// Free Software Foundation; either version 3, or (at your option) 9// any later version. 10 11// This library is distributed in the hope that it will be useful, 12// but WITHOUT ANY WARRANTY; without even the implied warranty of 13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14// GNU General Public License for more details. 15 16// Under Section 7 of GPL version 3, you are granted additional 17// permissions described in the GCC Runtime Library Exception, version 18// 3.1, as published by the Free Software Foundation. 19 20// You should have received a copy of the GNU General Public License and 21// a copy of the GCC Runtime Library Exception along with this program; 22// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 23// <http://www.gnu.org/licenses/>. 24 25/** @file ext/random 26 * This file is a GNU extension to the Standard C++ Library. 27 */ 28 29#ifndef _EXT_RANDOM 30#define _EXT_RANDOM 1 31 32#pragma GCC system_header 33 34#if __cplusplus < 201103L 35# include <bits/c++0x_warning.h> 36#else 37 38#include <random> 39#include <array> 40#include <ext/cmath> 41#ifdef __SSE2__ 42# include <x86intrin.h> 43#endif 44 45#ifdef _GLIBCXX_USE_C99_STDINT_TR1 46 47namespace __gnu_cxx _GLIBCXX_VISIBILITY(default) 48{ 49_GLIBCXX_BEGIN_NAMESPACE_VERSION 50 51#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 52 53 /* Mersenne twister implementation optimized for vector operations. 54 * 55 * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/ 56 */ 57 template<typename _UIntType, size_t __m, 58 size_t __pos1, size_t __sl1, size_t __sl2, 59 size_t __sr1, size_t __sr2, 60 uint32_t __msk1, uint32_t __msk2, 61 uint32_t __msk3, uint32_t __msk4, 62 uint32_t __parity1, uint32_t __parity2, 63 uint32_t __parity3, uint32_t __parity4> 64 class simd_fast_mersenne_twister_engine 65 { 66 static_assert(std::is_unsigned<_UIntType>::value, "template argument " 67 "substituting _UIntType not an unsigned integral type"); 68 static_assert(__sr1 < 32, "first right shift too large"); 69 static_assert(__sr2 < 16, "second right shift too large"); 70 static_assert(__sl1 < 32, "first left shift too large"); 71 static_assert(__sl2 < 16, "second left shift too large"); 72 73 public: 74 typedef _UIntType result_type; 75 76 private: 77 static constexpr size_t m_w = sizeof(result_type) * 8; 78 static constexpr size_t _M_nstate = __m / 128 + 1; 79 static constexpr size_t _M_nstate32 = _M_nstate * 4; 80 81 static_assert(std::is_unsigned<_UIntType>::value, "template argument " 82 "substituting _UIntType not an unsigned integral type"); 83 static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size"); 84 static_assert(16 % sizeof(_UIntType) == 0, 85 "UIntType size must divide 16"); 86 87 public: 88 static constexpr size_t state_size = _M_nstate * (16 89 / sizeof(result_type)); 90 static constexpr result_type default_seed = 5489u; 91 92 // constructors and member function 93 explicit 94 simd_fast_mersenne_twister_engine(result_type __sd = default_seed) 95 { seed(__sd); } 96 97 template<typename _Sseq, typename = typename 98 std::enable_if<!std::is_same<_Sseq, 99 simd_fast_mersenne_twister_engine>::value> 100 ::type> 101 explicit 102 simd_fast_mersenne_twister_engine(_Sseq& __q) 103 { seed(__q); } 104 105 void 106 seed(result_type __sd = default_seed); 107 108 template<typename _Sseq> 109 typename std::enable_if<std::is_class<_Sseq>::value>::type 110 seed(_Sseq& __q); 111 112 static constexpr result_type 113 min() 114 { return 0; }; 115 116 static constexpr result_type 117 max() 118 { return std::numeric_limits<result_type>::max(); } 119 120 void 121 discard(unsigned long long __z); 122 123 result_type 124 operator()() 125 { 126 if (__builtin_expect(_M_pos >= state_size, 0)) 127 _M_gen_rand(); 128 129 return _M_stateT[_M_pos++]; 130 } 131 132 template<typename _UIntType_2, size_t __m_2, 133 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2, 134 size_t __sr1_2, size_t __sr2_2, 135 uint32_t __msk1_2, uint32_t __msk2_2, 136 uint32_t __msk3_2, uint32_t __msk4_2, 137 uint32_t __parity1_2, uint32_t __parity2_2, 138 uint32_t __parity3_2, uint32_t __parity4_2> 139 friend bool 140 operator==(const simd_fast_mersenne_twister_engine<_UIntType_2, 141 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2, 142 __msk1_2, __msk2_2, __msk3_2, __msk4_2, 143 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs, 144 const simd_fast_mersenne_twister_engine<_UIntType_2, 145 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2, 146 __msk1_2, __msk2_2, __msk3_2, __msk4_2, 147 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs); 148 149 template<typename _UIntType_2, size_t __m_2, 150 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2, 151 size_t __sr1_2, size_t __sr2_2, 152 uint32_t __msk1_2, uint32_t __msk2_2, 153 uint32_t __msk3_2, uint32_t __msk4_2, 154 uint32_t __parity1_2, uint32_t __parity2_2, 155 uint32_t __parity3_2, uint32_t __parity4_2, 156 typename _CharT, typename _Traits> 157 friend std::basic_ostream<_CharT, _Traits>& 158 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 159 const __gnu_cxx::simd_fast_mersenne_twister_engine 160 <_UIntType_2, 161 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2, 162 __msk1_2, __msk2_2, __msk3_2, __msk4_2, 163 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x); 164 165 template<typename _UIntType_2, size_t __m_2, 166 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2, 167 size_t __sr1_2, size_t __sr2_2, 168 uint32_t __msk1_2, uint32_t __msk2_2, 169 uint32_t __msk3_2, uint32_t __msk4_2, 170 uint32_t __parity1_2, uint32_t __parity2_2, 171 uint32_t __parity3_2, uint32_t __parity4_2, 172 typename _CharT, typename _Traits> 173 friend std::basic_istream<_CharT, _Traits>& 174 operator>>(std::basic_istream<_CharT, _Traits>& __is, 175 __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2, 176 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2, 177 __msk1_2, __msk2_2, __msk3_2, __msk4_2, 178 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x); 179 180 private: 181 union 182 { 183#ifdef __SSE2__ 184 __m128i _M_state[_M_nstate]; 185#endif 186 uint32_t _M_state32[_M_nstate32]; 187 result_type _M_stateT[state_size]; 188 } __attribute__ ((__aligned__ (16))); 189 size_t _M_pos; 190 191 void _M_gen_rand(void); 192 void _M_period_certification(); 193 }; 194 195 196 template<typename _UIntType, size_t __m, 197 size_t __pos1, size_t __sl1, size_t __sl2, 198 size_t __sr1, size_t __sr2, 199 uint32_t __msk1, uint32_t __msk2, 200 uint32_t __msk3, uint32_t __msk4, 201 uint32_t __parity1, uint32_t __parity2, 202 uint32_t __parity3, uint32_t __parity4> 203 inline bool 204 operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 205 __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3, 206 __msk4, __parity1, __parity2, __parity3, __parity4>& __lhs, 207 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 208 __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3, 209 __msk4, __parity1, __parity2, __parity3, __parity4>& __rhs) 210 { return !(__lhs == __rhs); } 211 212 213 /* Definitions for the SIMD-oriented Fast Mersenne Twister as defined 214 * in the C implementation by Daito and Matsumoto, as both a 32-bit 215 * and 64-bit version. 216 */ 217 typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2, 218 15, 3, 13, 3, 219 0xfdff37ffU, 0xef7f3f7dU, 220 0xff777b7dU, 0x7ff7fb2fU, 221 0x00000001U, 0x00000000U, 222 0x00000000U, 0x5986f054U> 223 sfmt607; 224 225 typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2, 226 15, 3, 13, 3, 227 0xfdff37ffU, 0xef7f3f7dU, 228 0xff777b7dU, 0x7ff7fb2fU, 229 0x00000001U, 0x00000000U, 230 0x00000000U, 0x5986f054U> 231 sfmt607_64; 232 233 234 typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7, 235 14, 3, 5, 1, 236 0xf7fefffdU, 0x7fefcfffU, 237 0xaff3ef3fU, 0xb5ffff7fU, 238 0x00000001U, 0x00000000U, 239 0x00000000U, 0x20000000U> 240 sfmt1279; 241 242 typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7, 243 14, 3, 5, 1, 244 0xf7fefffdU, 0x7fefcfffU, 245 0xaff3ef3fU, 0xb5ffff7fU, 246 0x00000001U, 0x00000000U, 247 0x00000000U, 0x20000000U> 248 sfmt1279_64; 249 250 251 typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12, 252 19, 1, 5, 1, 253 0xbff7ffbfU, 0xfdfffffeU, 254 0xf7ffef7fU, 0xf2f7cbbfU, 255 0x00000001U, 0x00000000U, 256 0x00000000U, 0x41dfa600U> 257 sfmt2281; 258 259 typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12, 260 19, 1, 5, 1, 261 0xbff7ffbfU, 0xfdfffffeU, 262 0xf7ffef7fU, 0xf2f7cbbfU, 263 0x00000001U, 0x00000000U, 264 0x00000000U, 0x41dfa600U> 265 sfmt2281_64; 266 267 268 typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17, 269 20, 1, 7, 1, 270 0x9f7bffffU, 0x9fffff5fU, 271 0x3efffffbU, 0xfffff7bbU, 272 0xa8000001U, 0xaf5390a3U, 273 0xb740b3f8U, 0x6c11486dU> 274 sfmt4253; 275 276 typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17, 277 20, 1, 7, 1, 278 0x9f7bffffU, 0x9fffff5fU, 279 0x3efffffbU, 0xfffff7bbU, 280 0xa8000001U, 0xaf5390a3U, 281 0xb740b3f8U, 0x6c11486dU> 282 sfmt4253_64; 283 284 285 typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68, 286 14, 3, 7, 3, 287 0xeffff7fbU, 0xffffffefU, 288 0xdfdfbfffU, 0x7fffdbfdU, 289 0x00000001U, 0x00000000U, 290 0xe8148000U, 0xd0c7afa3U> 291 sfmt11213; 292 293 typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68, 294 14, 3, 7, 3, 295 0xeffff7fbU, 0xffffffefU, 296 0xdfdfbfffU, 0x7fffdbfdU, 297 0x00000001U, 0x00000000U, 298 0xe8148000U, 0xd0c7afa3U> 299 sfmt11213_64; 300 301 302 typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122, 303 18, 1, 11, 1, 304 0xdfffffefU, 0xddfecb7fU, 305 0xbffaffffU, 0xbffffff6U, 306 0x00000001U, 0x00000000U, 307 0x00000000U, 0x13c9e684U> 308 sfmt19937; 309 310 typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122, 311 18, 1, 11, 1, 312 0xdfffffefU, 0xddfecb7fU, 313 0xbffaffffU, 0xbffffff6U, 314 0x00000001U, 0x00000000U, 315 0x00000000U, 0x13c9e684U> 316 sfmt19937_64; 317 318 319 typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330, 320 5, 3, 9, 3, 321 0xeffffffbU, 0xdfbebfffU, 322 0xbfbf7befU, 0x9ffd7bffU, 323 0x00000001U, 0x00000000U, 324 0xa3ac4000U, 0xecc1327aU> 325 sfmt44497; 326 327 typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330, 328 5, 3, 9, 3, 329 0xeffffffbU, 0xdfbebfffU, 330 0xbfbf7befU, 0x9ffd7bffU, 331 0x00000001U, 0x00000000U, 332 0xa3ac4000U, 0xecc1327aU> 333 sfmt44497_64; 334 335 336 typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366, 337 6, 7, 19, 1, 338 0xfdbffbffU, 0xbff7ff3fU, 339 0xfd77efffU, 0xbf9ff3ffU, 340 0x00000001U, 0x00000000U, 341 0x00000000U, 0xe9528d85U> 342 sfmt86243; 343 344 typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366, 345 6, 7, 19, 1, 346 0xfdbffbffU, 0xbff7ff3fU, 347 0xfd77efffU, 0xbf9ff3ffU, 348 0x00000001U, 0x00000000U, 349 0x00000000U, 0xe9528d85U> 350 sfmt86243_64; 351 352 353 typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110, 354 19, 1, 21, 1, 355 0xffffbb5fU, 0xfb6ebf95U, 356 0xfffefffaU, 0xcff77fffU, 357 0x00000001U, 0x00000000U, 358 0xcb520000U, 0xc7e91c7dU> 359 sfmt132049; 360 361 typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110, 362 19, 1, 21, 1, 363 0xffffbb5fU, 0xfb6ebf95U, 364 0xfffefffaU, 0xcff77fffU, 365 0x00000001U, 0x00000000U, 366 0xcb520000U, 0xc7e91c7dU> 367 sfmt132049_64; 368 369 370 typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627, 371 11, 3, 10, 1, 372 0xbff7bff7U, 0xbfffffffU, 373 0xbffffa7fU, 0xffddfbfbU, 374 0xf8000001U, 0x89e80709U, 375 0x3bd2b64bU, 0x0c64b1e4U> 376 sfmt216091; 377 378 typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627, 379 11, 3, 10, 1, 380 0xbff7bff7U, 0xbfffffffU, 381 0xbffffa7fU, 0xffddfbfbU, 382 0xf8000001U, 0x89e80709U, 383 0x3bd2b64bU, 0x0c64b1e4U> 384 sfmt216091_64; 385 386#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 387 388 /** 389 * @brief A beta continuous distribution for random numbers. 390 * 391 * The formula for the beta probability density function is: 392 * @f[ 393 * p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)} 394 * x^{\alpha - 1} (1 - x)^{\beta - 1} 395 * @f] 396 */ 397 template<typename _RealType = double> 398 class beta_distribution 399 { 400 static_assert(std::is_floating_point<_RealType>::value, 401 "template argument not a floating point type"); 402 403 public: 404 /** The type of the range of the distribution. */ 405 typedef _RealType result_type; 406 /** Parameter type. */ 407 struct param_type 408 { 409 typedef beta_distribution<_RealType> distribution_type; 410 friend class beta_distribution<_RealType>; 411 412 explicit 413 param_type(_RealType __alpha_val = _RealType(1), 414 _RealType __beta_val = _RealType(1)) 415 : _M_alpha(__alpha_val), _M_beta(__beta_val) 416 { 417 _GLIBCXX_DEBUG_ASSERT(_M_alpha > _RealType(0)); 418 _GLIBCXX_DEBUG_ASSERT(_M_beta > _RealType(0)); 419 } 420 421 _RealType 422 alpha() const 423 { return _M_alpha; } 424 425 _RealType 426 beta() const 427 { return _M_beta; } 428 429 friend bool 430 operator==(const param_type& __p1, const param_type& __p2) 431 { return (__p1._M_alpha == __p2._M_alpha 432 && __p1._M_beta == __p2._M_beta); } 433 434 private: 435 void 436 _M_initialize(); 437 438 _RealType _M_alpha; 439 _RealType _M_beta; 440 }; 441 442 public: 443 /** 444 * @brief Constructs a beta distribution with parameters 445 * @f$\alpha@f$ and @f$\beta@f$. 446 */ 447 explicit 448 beta_distribution(_RealType __alpha_val = _RealType(1), 449 _RealType __beta_val = _RealType(1)) 450 : _M_param(__alpha_val, __beta_val) 451 { } 452 453 explicit 454 beta_distribution(const param_type& __p) 455 : _M_param(__p) 456 { } 457 458 /** 459 * @brief Resets the distribution state. 460 */ 461 void 462 reset() 463 { } 464 465 /** 466 * @brief Returns the @f$\alpha@f$ of the distribution. 467 */ 468 _RealType 469 alpha() const 470 { return _M_param.alpha(); } 471 472 /** 473 * @brief Returns the @f$\beta@f$ of the distribution. 474 */ 475 _RealType 476 beta() const 477 { return _M_param.beta(); } 478 479 /** 480 * @brief Returns the parameter set of the distribution. 481 */ 482 param_type 483 param() const 484 { return _M_param; } 485 486 /** 487 * @brief Sets the parameter set of the distribution. 488 * @param __param The new parameter set of the distribution. 489 */ 490 void 491 param(const param_type& __param) 492 { _M_param = __param; } 493 494 /** 495 * @brief Returns the greatest lower bound value of the distribution. 496 */ 497 result_type 498 min() const 499 { return result_type(0); } 500 501 /** 502 * @brief Returns the least upper bound value of the distribution. 503 */ 504 result_type 505 max() const 506 { return result_type(1); } 507 508 /** 509 * @brief Generating functions. 510 */ 511 template<typename _UniformRandomNumberGenerator> 512 result_type 513 operator()(_UniformRandomNumberGenerator& __urng) 514 { return this->operator()(__urng, _M_param); } 515 516 template<typename _UniformRandomNumberGenerator> 517 result_type 518 operator()(_UniformRandomNumberGenerator& __urng, 519 const param_type& __p); 520 521 template<typename _ForwardIterator, 522 typename _UniformRandomNumberGenerator> 523 void 524 __generate(_ForwardIterator __f, _ForwardIterator __t, 525 _UniformRandomNumberGenerator& __urng) 526 { this->__generate(__f, __t, __urng, _M_param); } 527 528 template<typename _ForwardIterator, 529 typename _UniformRandomNumberGenerator> 530 void 531 __generate(_ForwardIterator __f, _ForwardIterator __t, 532 _UniformRandomNumberGenerator& __urng, 533 const param_type& __p) 534 { this->__generate_impl(__f, __t, __urng, __p); } 535 536 template<typename _UniformRandomNumberGenerator> 537 void 538 __generate(result_type* __f, result_type* __t, 539 _UniformRandomNumberGenerator& __urng, 540 const param_type& __p) 541 { this->__generate_impl(__f, __t, __urng, __p); } 542 543 /** 544 * @brief Return true if two beta distributions have the same 545 * parameters and the sequences that would be generated 546 * are equal. 547 */ 548 friend bool 549 operator==(const beta_distribution& __d1, 550 const beta_distribution& __d2) 551 { return __d1._M_param == __d2._M_param; } 552 553 /** 554 * @brief Inserts a %beta_distribution random number distribution 555 * @p __x into the output stream @p __os. 556 * 557 * @param __os An output stream. 558 * @param __x A %beta_distribution random number distribution. 559 * 560 * @returns The output stream with the state of @p __x inserted or in 561 * an error state. 562 */ 563 template<typename _RealType1, typename _CharT, typename _Traits> 564 friend std::basic_ostream<_CharT, _Traits>& 565 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 566 const __gnu_cxx::beta_distribution<_RealType1>& __x); 567 568 /** 569 * @brief Extracts a %beta_distribution random number distribution 570 * @p __x from the input stream @p __is. 571 * 572 * @param __is An input stream. 573 * @param __x A %beta_distribution random number generator engine. 574 * 575 * @returns The input stream with @p __x extracted or in an error state. 576 */ 577 template<typename _RealType1, typename _CharT, typename _Traits> 578 friend std::basic_istream<_CharT, _Traits>& 579 operator>>(std::basic_istream<_CharT, _Traits>& __is, 580 __gnu_cxx::beta_distribution<_RealType1>& __x); 581 582 private: 583 template<typename _ForwardIterator, 584 typename _UniformRandomNumberGenerator> 585 void 586 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 587 _UniformRandomNumberGenerator& __urng, 588 const param_type& __p); 589 590 param_type _M_param; 591 }; 592 593 /** 594 * @brief Return true if two beta distributions are different. 595 */ 596 template<typename _RealType> 597 inline bool 598 operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1, 599 const __gnu_cxx::beta_distribution<_RealType>& __d2) 600 { return !(__d1 == __d2); } 601 602 603 /** 604 * @brief A multi-variate normal continuous distribution for random numbers. 605 * 606 * The formula for the normal probability density function is 607 * @f[ 608 * p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) = 609 * \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}} 610 * e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T} 611 * \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})} 612 * @f] 613 * 614 * where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are 615 * vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance 616 * matrix (which must be positive-definite). 617 */ 618 template<std::size_t _Dimen, typename _RealType = double> 619 class normal_mv_distribution 620 { 621 static_assert(std::is_floating_point<_RealType>::value, 622 "template argument not a floating point type"); 623 static_assert(_Dimen != 0, "dimension is zero"); 624 625 public: 626 /** The type of the range of the distribution. */ 627 typedef std::array<_RealType, _Dimen> result_type; 628 /** Parameter type. */ 629 class param_type 630 { 631 static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2; 632 633 public: 634 typedef normal_mv_distribution<_Dimen, _RealType> distribution_type; 635 friend class normal_mv_distribution<_Dimen, _RealType>; 636 637 param_type() 638 { 639 std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0)); 640 auto __it = _M_t.begin(); 641 for (size_t __i = 0; __i < _Dimen; ++__i) 642 { 643 std::fill_n(__it, __i, _RealType(0)); 644 __it += __i; 645 *__it++ = _RealType(1); 646 } 647 } 648 649 template<typename _ForwardIterator1, typename _ForwardIterator2> 650 param_type(_ForwardIterator1 __meanbegin, 651 _ForwardIterator1 __meanend, 652 _ForwardIterator2 __varcovbegin, 653 _ForwardIterator2 __varcovend) 654 { 655 __glibcxx_function_requires(_ForwardIteratorConcept< 656 _ForwardIterator1>) 657 __glibcxx_function_requires(_ForwardIteratorConcept< 658 _ForwardIterator2>) 659 _GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend) 660 <= _Dimen); 661 const auto __dist = std::distance(__varcovbegin, __varcovend); 662 _GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen 663 || __dist == _Dimen * (_Dimen + 1) / 2 664 || __dist == _Dimen); 665 666 if (__dist == _Dimen * _Dimen) 667 _M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend); 668 else if (__dist == _Dimen * (_Dimen + 1) / 2) 669 _M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend); 670 else 671 _M_init_diagonal(__meanbegin, __meanend, 672 __varcovbegin, __varcovend); 673 } 674 675 param_type(std::initializer_list<_RealType> __mean, 676 std::initializer_list<_RealType> __varcov) 677 { 678 _GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen); 679 _GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen 680 || __varcov.size() == _Dimen * (_Dimen + 1) / 2 681 || __varcov.size() == _Dimen); 682 683 if (__varcov.size() == _Dimen * _Dimen) 684 _M_init_full(__mean.begin(), __mean.end(), 685 __varcov.begin(), __varcov.end()); 686 else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2) 687 _M_init_lower(__mean.begin(), __mean.end(), 688 __varcov.begin(), __varcov.end()); 689 else 690 _M_init_diagonal(__mean.begin(), __mean.end(), 691 __varcov.begin(), __varcov.end()); 692 } 693 694 std::array<_RealType, _Dimen> 695 mean() const 696 { return _M_mean; } 697 698 std::array<_RealType, _M_t_size> 699 varcov() const 700 { return _M_t; } 701 702 friend bool 703 operator==(const param_type& __p1, const param_type& __p2) 704 { return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; } 705 706 private: 707 template <typename _InputIterator1, typename _InputIterator2> 708 void _M_init_full(_InputIterator1 __meanbegin, 709 _InputIterator1 __meanend, 710 _InputIterator2 __varcovbegin, 711 _InputIterator2 __varcovend); 712 template <typename _InputIterator1, typename _InputIterator2> 713 void _M_init_lower(_InputIterator1 __meanbegin, 714 _InputIterator1 __meanend, 715 _InputIterator2 __varcovbegin, 716 _InputIterator2 __varcovend); 717 template <typename _InputIterator1, typename _InputIterator2> 718 void _M_init_diagonal(_InputIterator1 __meanbegin, 719 _InputIterator1 __meanend, 720 _InputIterator2 __varbegin, 721 _InputIterator2 __varend); 722 723 std::array<_RealType, _Dimen> _M_mean; 724 std::array<_RealType, _M_t_size> _M_t; 725 }; 726 727 public: 728 normal_mv_distribution() 729 : _M_param(), _M_nd() 730 { } 731 732 template<typename _ForwardIterator1, typename _ForwardIterator2> 733 normal_mv_distribution(_ForwardIterator1 __meanbegin, 734 _ForwardIterator1 __meanend, 735 _ForwardIterator2 __varcovbegin, 736 _ForwardIterator2 __varcovend) 737 : _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend), 738 _M_nd() 739 { } 740 741 normal_mv_distribution(std::initializer_list<_RealType> __mean, 742 std::initializer_list<_RealType> __varcov) 743 : _M_param(__mean, __varcov), _M_nd() 744 { } 745 746 explicit 747 normal_mv_distribution(const param_type& __p) 748 : _M_param(__p), _M_nd() 749 { } 750 751 /** 752 * @brief Resets the distribution state. 753 */ 754 void 755 reset() 756 { _M_nd.reset(); } 757 758 /** 759 * @brief Returns the mean of the distribution. 760 */ 761 result_type 762 mean() const 763 { return _M_param.mean(); } 764 765 /** 766 * @brief Returns the compact form of the variance/covariance 767 * matrix of the distribution. 768 */ 769 std::array<_RealType, _Dimen * (_Dimen + 1) / 2> 770 varcov() const 771 { return _M_param.varcov(); } 772 773 /** 774 * @brief Returns the parameter set of the distribution. 775 */ 776 param_type 777 param() const 778 { return _M_param; } 779 780 /** 781 * @brief Sets the parameter set of the distribution. 782 * @param __param The new parameter set of the distribution. 783 */ 784 void 785 param(const param_type& __param) 786 { _M_param = __param; } 787 788 /** 789 * @brief Returns the greatest lower bound value of the distribution. 790 */ 791 result_type 792 min() const 793 { result_type __res; 794 __res.fill(std::numeric_limits<_RealType>::lowest()); 795 return __res; } 796 797 /** 798 * @brief Returns the least upper bound value of the distribution. 799 */ 800 result_type 801 max() const 802 { result_type __res; 803 __res.fill(std::numeric_limits<_RealType>::max()); 804 return __res; } 805 806 /** 807 * @brief Generating functions. 808 */ 809 template<typename _UniformRandomNumberGenerator> 810 result_type 811 operator()(_UniformRandomNumberGenerator& __urng) 812 { return this->operator()(__urng, _M_param); } 813 814 template<typename _UniformRandomNumberGenerator> 815 result_type 816 operator()(_UniformRandomNumberGenerator& __urng, 817 const param_type& __p); 818 819 template<typename _ForwardIterator, 820 typename _UniformRandomNumberGenerator> 821 void 822 __generate(_ForwardIterator __f, _ForwardIterator __t, 823 _UniformRandomNumberGenerator& __urng) 824 { return this->__generate_impl(__f, __t, __urng, _M_param); } 825 826 template<typename _ForwardIterator, 827 typename _UniformRandomNumberGenerator> 828 void 829 __generate(_ForwardIterator __f, _ForwardIterator __t, 830 _UniformRandomNumberGenerator& __urng, 831 const param_type& __p) 832 { return this->__generate_impl(__f, __t, __urng, __p); } 833 834 /** 835 * @brief Return true if two multi-variant normal distributions have 836 * the same parameters and the sequences that would 837 * be generated are equal. 838 */ 839 template<size_t _Dimen1, typename _RealType1> 840 friend bool 841 operator==(const 842 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>& 843 __d1, 844 const 845 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>& 846 __d2); 847 848 /** 849 * @brief Inserts a %normal_mv_distribution random number distribution 850 * @p __x into the output stream @p __os. 851 * 852 * @param __os An output stream. 853 * @param __x A %normal_mv_distribution random number distribution. 854 * 855 * @returns The output stream with the state of @p __x inserted or in 856 * an error state. 857 */ 858 template<size_t _Dimen1, typename _RealType1, 859 typename _CharT, typename _Traits> 860 friend std::basic_ostream<_CharT, _Traits>& 861 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 862 const 863 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>& 864 __x); 865 866 /** 867 * @brief Extracts a %normal_mv_distribution random number distribution 868 * @p __x from the input stream @p __is. 869 * 870 * @param __is An input stream. 871 * @param __x A %normal_mv_distribution random number generator engine. 872 * 873 * @returns The input stream with @p __x extracted or in an error 874 * state. 875 */ 876 template<size_t _Dimen1, typename _RealType1, 877 typename _CharT, typename _Traits> 878 friend std::basic_istream<_CharT, _Traits>& 879 operator>>(std::basic_istream<_CharT, _Traits>& __is, 880 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>& 881 __x); 882 883 private: 884 template<typename _ForwardIterator, 885 typename _UniformRandomNumberGenerator> 886 void 887 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 888 _UniformRandomNumberGenerator& __urng, 889 const param_type& __p); 890 891 param_type _M_param; 892 std::normal_distribution<_RealType> _M_nd; 893 }; 894 895 /** 896 * @brief Return true if two multi-variate normal distributions are 897 * different. 898 */ 899 template<size_t _Dimen, typename _RealType> 900 inline bool 901 operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 902 __d1, 903 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 904 __d2) 905 { return !(__d1 == __d2); } 906 907 908 /** 909 * @brief A Rice continuous distribution for random numbers. 910 * 911 * The formula for the Rice probability density function is 912 * @f[ 913 * p(x|\nu,\sigma) = \frac{x}{\sigma^2} 914 * \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right) 915 * I_0\left(\frac{x \nu}{\sigma^2}\right) 916 * @f] 917 * where @f$I_0(z)@f$ is the modified Bessel function of the first kind 918 * of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$. 919 * 920 * <table border=1 cellpadding=10 cellspacing=0> 921 * <caption align=top>Distribution Statistics</caption> 922 * <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr> 923 * <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2 924 * + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr> 925 * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr> 926 * </table> 927 * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2. 928 */ 929 template<typename _RealType = double> 930 class 931 rice_distribution 932 { 933 static_assert(std::is_floating_point<_RealType>::value, 934 "template argument not a floating point type"); 935 public: 936 /** The type of the range of the distribution. */ 937 typedef _RealType result_type; 938 /** Parameter type. */ 939 struct param_type 940 { 941 typedef rice_distribution<result_type> distribution_type; 942 943 param_type(result_type __nu_val = result_type(0), 944 result_type __sigma_val = result_type(1)) 945 : _M_nu(__nu_val), _M_sigma(__sigma_val) 946 { 947 _GLIBCXX_DEBUG_ASSERT(_M_nu >= result_type(0)); 948 _GLIBCXX_DEBUG_ASSERT(_M_sigma > result_type(0)); 949 } 950 951 result_type 952 nu() const 953 { return _M_nu; } 954 955 result_type 956 sigma() const 957 { return _M_sigma; } 958 959 friend bool 960 operator==(const param_type& __p1, const param_type& __p2) 961 { return __p1._M_nu == __p2._M_nu 962 && __p1._M_sigma == __p2._M_sigma; } 963 964 private: 965 void _M_initialize(); 966 967 result_type _M_nu; 968 result_type _M_sigma; 969 }; 970 971 /** 972 * @brief Constructors. 973 */ 974 explicit 975 rice_distribution(result_type __nu_val = result_type(0), 976 result_type __sigma_val = result_type(1)) 977 : _M_param(__nu_val, __sigma_val), 978 _M_ndx(__nu_val, __sigma_val), 979 _M_ndy(result_type(0), __sigma_val) 980 { } 981 982 explicit 983 rice_distribution(const param_type& __p) 984 : _M_param(__p), 985 _M_ndx(__p.nu(), __p.sigma()), 986 _M_ndy(result_type(0), __p.sigma()) 987 { } 988 989 /** 990 * @brief Resets the distribution state. 991 */ 992 void 993 reset() 994 { 995 _M_ndx.reset(); 996 _M_ndy.reset(); 997 } 998 999 /** 1000 * @brief Return the parameters of the distribution. 1001 */ 1002 result_type 1003 nu() const 1004 { return _M_param.nu(); } 1005 1006 result_type 1007 sigma() const 1008 { return _M_param.sigma(); } 1009 1010 /** 1011 * @brief Returns the parameter set of the distribution. 1012 */ 1013 param_type 1014 param() const 1015 { return _M_param; } 1016 1017 /** 1018 * @brief Sets the parameter set of the distribution. 1019 * @param __param The new parameter set of the distribution. 1020 */ 1021 void 1022 param(const param_type& __param) 1023 { _M_param = __param; } 1024 1025 /** 1026 * @brief Returns the greatest lower bound value of the distribution. 1027 */ 1028 result_type 1029 min() const 1030 { return result_type(0); } 1031 1032 /** 1033 * @brief Returns the least upper bound value of the distribution. 1034 */ 1035 result_type 1036 max() const 1037 { return std::numeric_limits<result_type>::max(); } 1038 1039 /** 1040 * @brief Generating functions. 1041 */ 1042 template<typename _UniformRandomNumberGenerator> 1043 result_type 1044 operator()(_UniformRandomNumberGenerator& __urng) 1045 { 1046 result_type __x = this->_M_ndx(__urng); 1047 result_type __y = this->_M_ndy(__urng); 1048#if _GLIBCXX_USE_C99_MATH_TR1 1049 return std::hypot(__x, __y); 1050#else 1051 return std::sqrt(__x * __x + __y * __y); 1052#endif 1053 } 1054 1055 template<typename _UniformRandomNumberGenerator> 1056 result_type 1057 operator()(_UniformRandomNumberGenerator& __urng, 1058 const param_type& __p) 1059 { 1060 typename std::normal_distribution<result_type>::param_type 1061 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma()); 1062 result_type __x = this->_M_ndx(__px, __urng); 1063 result_type __y = this->_M_ndy(__py, __urng); 1064#if _GLIBCXX_USE_C99_MATH_TR1 1065 return std::hypot(__x, __y); 1066#else 1067 return std::sqrt(__x * __x + __y * __y); 1068#endif 1069 } 1070 1071 template<typename _ForwardIterator, 1072 typename _UniformRandomNumberGenerator> 1073 void 1074 __generate(_ForwardIterator __f, _ForwardIterator __t, 1075 _UniformRandomNumberGenerator& __urng) 1076 { this->__generate(__f, __t, __urng, _M_param); } 1077 1078 template<typename _ForwardIterator, 1079 typename _UniformRandomNumberGenerator> 1080 void 1081 __generate(_ForwardIterator __f, _ForwardIterator __t, 1082 _UniformRandomNumberGenerator& __urng, 1083 const param_type& __p) 1084 { this->__generate_impl(__f, __t, __urng, __p); } 1085 1086 template<typename _UniformRandomNumberGenerator> 1087 void 1088 __generate(result_type* __f, result_type* __t, 1089 _UniformRandomNumberGenerator& __urng, 1090 const param_type& __p) 1091 { this->__generate_impl(__f, __t, __urng, __p); } 1092 1093 /** 1094 * @brief Return true if two Rice distributions have 1095 * the same parameters and the sequences that would 1096 * be generated are equal. 1097 */ 1098 friend bool 1099 operator==(const rice_distribution& __d1, 1100 const rice_distribution& __d2) 1101 { return (__d1._M_param == __d2._M_param 1102 && __d1._M_ndx == __d2._M_ndx 1103 && __d1._M_ndy == __d2._M_ndy); } 1104 1105 /** 1106 * @brief Inserts a %rice_distribution random number distribution 1107 * @p __x into the output stream @p __os. 1108 * 1109 * @param __os An output stream. 1110 * @param __x A %rice_distribution random number distribution. 1111 * 1112 * @returns The output stream with the state of @p __x inserted or in 1113 * an error state. 1114 */ 1115 template<typename _RealType1, typename _CharT, typename _Traits> 1116 friend std::basic_ostream<_CharT, _Traits>& 1117 operator<<(std::basic_ostream<_CharT, _Traits>&, 1118 const rice_distribution<_RealType1>&); 1119 1120 /** 1121 * @brief Extracts a %rice_distribution random number distribution 1122 * @p __x from the input stream @p __is. 1123 * 1124 * @param __is An input stream. 1125 * @param __x A %rice_distribution random number 1126 * generator engine. 1127 * 1128 * @returns The input stream with @p __x extracted or in an error state. 1129 */ 1130 template<typename _RealType1, typename _CharT, typename _Traits> 1131 friend std::basic_istream<_CharT, _Traits>& 1132 operator>>(std::basic_istream<_CharT, _Traits>&, 1133 rice_distribution<_RealType1>&); 1134 1135 private: 1136 template<typename _ForwardIterator, 1137 typename _UniformRandomNumberGenerator> 1138 void 1139 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 1140 _UniformRandomNumberGenerator& __urng, 1141 const param_type& __p); 1142 1143 param_type _M_param; 1144 1145 std::normal_distribution<result_type> _M_ndx; 1146 std::normal_distribution<result_type> _M_ndy; 1147 }; 1148 1149 /** 1150 * @brief Return true if two Rice distributions are not equal. 1151 */ 1152 template<typename _RealType1> 1153 inline bool 1154 operator!=(const rice_distribution<_RealType1>& __d1, 1155 const rice_distribution<_RealType1>& __d2) 1156 { return !(__d1 == __d2); } 1157 1158 1159 /** 1160 * @brief A Nakagami continuous distribution for random numbers. 1161 * 1162 * The formula for the Nakagami probability density function is 1163 * @f[ 1164 * p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu} 1165 * x^{2\mu-1}e^{-\mu x / \omega} 1166 * @f] 1167 * where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$ 1168 * and @f$\omega > 0@f$. 1169 */ 1170 template<typename _RealType = double> 1171 class 1172 nakagami_distribution 1173 { 1174 static_assert(std::is_floating_point<_RealType>::value, 1175 "template argument not a floating point type"); 1176 1177 public: 1178 /** The type of the range of the distribution. */ 1179 typedef _RealType result_type; 1180 /** Parameter type. */ 1181 struct param_type 1182 { 1183 typedef nakagami_distribution<result_type> distribution_type; 1184 1185 param_type(result_type __mu_val = result_type(1), 1186 result_type __omega_val = result_type(1)) 1187 : _M_mu(__mu_val), _M_omega(__omega_val) 1188 { 1189 _GLIBCXX_DEBUG_ASSERT(_M_mu >= result_type(0.5L)); 1190 _GLIBCXX_DEBUG_ASSERT(_M_omega > result_type(0)); 1191 } 1192 1193 result_type 1194 mu() const 1195 { return _M_mu; } 1196 1197 result_type 1198 omega() const 1199 { return _M_omega; } 1200 1201 friend bool 1202 operator==(const param_type& __p1, const param_type& __p2) 1203 { return __p1._M_mu == __p2._M_mu 1204 && __p1._M_omega == __p2._M_omega; } 1205 1206 private: 1207 void _M_initialize(); 1208 1209 result_type _M_mu; 1210 result_type _M_omega; 1211 }; 1212 1213 /** 1214 * @brief Constructors. 1215 */ 1216 explicit 1217 nakagami_distribution(result_type __mu_val = result_type(1), 1218 result_type __omega_val = result_type(1)) 1219 : _M_param(__mu_val, __omega_val), 1220 _M_gd(__mu_val, __omega_val / __mu_val) 1221 { } 1222 1223 explicit 1224 nakagami_distribution(const param_type& __p) 1225 : _M_param(__p), 1226 _M_gd(__p.mu(), __p.omega() / __p.mu()) 1227 { } 1228 1229 /** 1230 * @brief Resets the distribution state. 1231 */ 1232 void 1233 reset() 1234 { _M_gd.reset(); } 1235 1236 /** 1237 * @brief Return the parameters of the distribution. 1238 */ 1239 result_type 1240 mu() const 1241 { return _M_param.mu(); } 1242 1243 result_type 1244 omega() const 1245 { return _M_param.omega(); } 1246 1247 /** 1248 * @brief Returns the parameter set of the distribution. 1249 */ 1250 param_type 1251 param() const 1252 { return _M_param; } 1253 1254 /** 1255 * @brief Sets the parameter set of the distribution. 1256 * @param __param The new parameter set of the distribution. 1257 */ 1258 void 1259 param(const param_type& __param) 1260 { _M_param = __param; } 1261 1262 /** 1263 * @brief Returns the greatest lower bound value of the distribution. 1264 */ 1265 result_type 1266 min() const 1267 { return result_type(0); } 1268 1269 /** 1270 * @brief Returns the least upper bound value of the distribution. 1271 */ 1272 result_type 1273 max() const 1274 { return std::numeric_limits<result_type>::max(); } 1275 1276 /** 1277 * @brief Generating functions. 1278 */ 1279 template<typename _UniformRandomNumberGenerator> 1280 result_type 1281 operator()(_UniformRandomNumberGenerator& __urng) 1282 { return std::sqrt(this->_M_gd(__urng)); } 1283 1284 template<typename _UniformRandomNumberGenerator> 1285 result_type 1286 operator()(_UniformRandomNumberGenerator& __urng, 1287 const param_type& __p) 1288 { 1289 typename std::gamma_distribution<result_type>::param_type 1290 __pg(__p.mu(), __p.omega() / __p.mu()); 1291 return std::sqrt(this->_M_gd(__pg, __urng)); 1292 } 1293 1294 template<typename _ForwardIterator, 1295 typename _UniformRandomNumberGenerator> 1296 void 1297 __generate(_ForwardIterator __f, _ForwardIterator __t, 1298 _UniformRandomNumberGenerator& __urng) 1299 { this->__generate(__f, __t, __urng, _M_param); } 1300 1301 template<typename _ForwardIterator, 1302 typename _UniformRandomNumberGenerator> 1303 void 1304 __generate(_ForwardIterator __f, _ForwardIterator __t, 1305 _UniformRandomNumberGenerator& __urng, 1306 const param_type& __p) 1307 { this->__generate_impl(__f, __t, __urng, __p); } 1308 1309 template<typename _UniformRandomNumberGenerator> 1310 void 1311 __generate(result_type* __f, result_type* __t, 1312 _UniformRandomNumberGenerator& __urng, 1313 const param_type& __p) 1314 { this->__generate_impl(__f, __t, __urng, __p); } 1315 1316 /** 1317 * @brief Return true if two Nakagami distributions have 1318 * the same parameters and the sequences that would 1319 * be generated are equal. 1320 */ 1321 friend bool 1322 operator==(const nakagami_distribution& __d1, 1323 const nakagami_distribution& __d2) 1324 { return (__d1._M_param == __d2._M_param 1325 && __d1._M_gd == __d2._M_gd); } 1326 1327 /** 1328 * @brief Inserts a %nakagami_distribution random number distribution 1329 * @p __x into the output stream @p __os. 1330 * 1331 * @param __os An output stream. 1332 * @param __x A %nakagami_distribution random number distribution. 1333 * 1334 * @returns The output stream with the state of @p __x inserted or in 1335 * an error state. 1336 */ 1337 template<typename _RealType1, typename _CharT, typename _Traits> 1338 friend std::basic_ostream<_CharT, _Traits>& 1339 operator<<(std::basic_ostream<_CharT, _Traits>&, 1340 const nakagami_distribution<_RealType1>&); 1341 1342 /** 1343 * @brief Extracts a %nakagami_distribution random number distribution 1344 * @p __x from the input stream @p __is. 1345 * 1346 * @param __is An input stream. 1347 * @param __x A %nakagami_distribution random number 1348 * generator engine. 1349 * 1350 * @returns The input stream with @p __x extracted or in an error state. 1351 */ 1352 template<typename _RealType1, typename _CharT, typename _Traits> 1353 friend std::basic_istream<_CharT, _Traits>& 1354 operator>>(std::basic_istream<_CharT, _Traits>&, 1355 nakagami_distribution<_RealType1>&); 1356 1357 private: 1358 template<typename _ForwardIterator, 1359 typename _UniformRandomNumberGenerator> 1360 void 1361 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 1362 _UniformRandomNumberGenerator& __urng, 1363 const param_type& __p); 1364 1365 param_type _M_param; 1366 1367 std::gamma_distribution<result_type> _M_gd; 1368 }; 1369 1370 /** 1371 * @brief Return true if two Nakagami distributions are not equal. 1372 */ 1373 template<typename _RealType> 1374 inline bool 1375 operator!=(const nakagami_distribution<_RealType>& __d1, 1376 const nakagami_distribution<_RealType>& __d2) 1377 { return !(__d1 == __d2); } 1378 1379 1380 /** 1381 * @brief A Pareto continuous distribution for random numbers. 1382 * 1383 * The formula for the Pareto cumulative probability function is 1384 * @f[ 1385 * P(x|\alpha,\mu) = 1 - \left(\frac{\mu}{x}\right)^\alpha 1386 * @f] 1387 * The formula for the Pareto probability density function is 1388 * @f[ 1389 * p(x|\alpha,\mu) = \frac{\alpha + 1}{\mu} 1390 * \left(\frac{\mu}{x}\right)^{\alpha + 1} 1391 * @f] 1392 * where @f$x >= \mu@f$ and @f$\mu > 0@f$, @f$\alpha > 0@f$. 1393 * 1394 * <table border=1 cellpadding=10 cellspacing=0> 1395 * <caption align=top>Distribution Statistics</caption> 1396 * <tr><td>Mean</td><td>@f$\alpha \mu / (\alpha - 1)@f$ 1397 * for @f$\alpha > 1@f$</td></tr> 1398 * <tr><td>Variance</td><td>@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$ 1399 * for @f$\alpha > 2@f$</td></tr> 1400 * <tr><td>Range</td><td>@f$[\mu, \infty)@f$</td></tr> 1401 * </table> 1402 */ 1403 template<typename _RealType = double> 1404 class 1405 pareto_distribution 1406 { 1407 static_assert(std::is_floating_point<_RealType>::value, 1408 "template argument not a floating point type"); 1409 1410 public: 1411 /** The type of the range of the distribution. */ 1412 typedef _RealType result_type; 1413 /** Parameter type. */ 1414 struct param_type 1415 { 1416 typedef pareto_distribution<result_type> distribution_type; 1417 1418 param_type(result_type __alpha_val = result_type(1), 1419 result_type __mu_val = result_type(1)) 1420 : _M_alpha(__alpha_val), _M_mu(__mu_val) 1421 { 1422 _GLIBCXX_DEBUG_ASSERT(_M_alpha > result_type(0)); 1423 _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0)); 1424 } 1425 1426 result_type 1427 alpha() const 1428 { return _M_alpha; } 1429 1430 result_type 1431 mu() const 1432 { return _M_mu; } 1433 1434 friend bool 1435 operator==(const param_type& __p1, const param_type& __p2) 1436 { return __p1._M_alpha == __p2._M_alpha && __p1._M_mu == __p2._M_mu; } 1437 1438 private: 1439 void _M_initialize(); 1440 1441 result_type _M_alpha; 1442 result_type _M_mu; 1443 }; 1444 1445 /** 1446 * @brief Constructors. 1447 */ 1448 explicit 1449 pareto_distribution(result_type __alpha_val = result_type(1), 1450 result_type __mu_val = result_type(1)) 1451 : _M_param(__alpha_val, __mu_val), 1452 _M_ud() 1453 { } 1454 1455 explicit 1456 pareto_distribution(const param_type& __p) 1457 : _M_param(__p), 1458 _M_ud() 1459 { } 1460 1461 /** 1462 * @brief Resets the distribution state. 1463 */ 1464 void 1465 reset() 1466 { 1467 _M_ud.reset(); 1468 } 1469 1470 /** 1471 * @brief Return the parameters of the distribution. 1472 */ 1473 result_type 1474 alpha() const 1475 { return _M_param.alpha(); } 1476 1477 result_type 1478 mu() const 1479 { return _M_param.mu(); } 1480 1481 /** 1482 * @brief Returns the parameter set of the distribution. 1483 */ 1484 param_type 1485 param() const 1486 { return _M_param; } 1487 1488 /** 1489 * @brief Sets the parameter set of the distribution. 1490 * @param __param The new parameter set of the distribution. 1491 */ 1492 void 1493 param(const param_type& __param) 1494 { _M_param = __param; } 1495 1496 /** 1497 * @brief Returns the greatest lower bound value of the distribution. 1498 */ 1499 result_type 1500 min() const 1501 { return this->mu(); } 1502 1503 /** 1504 * @brief Returns the least upper bound value of the distribution. 1505 */ 1506 result_type 1507 max() const 1508 { return std::numeric_limits<result_type>::max(); } 1509 1510 /** 1511 * @brief Generating functions. 1512 */ 1513 template<typename _UniformRandomNumberGenerator> 1514 result_type 1515 operator()(_UniformRandomNumberGenerator& __urng) 1516 { 1517 return this->mu() * std::pow(this->_M_ud(__urng), 1518 -result_type(1) / this->alpha()); 1519 } 1520 1521 template<typename _UniformRandomNumberGenerator> 1522 result_type 1523 operator()(_UniformRandomNumberGenerator& __urng, 1524 const param_type& __p) 1525 { 1526 return __p.mu() * std::pow(this->_M_ud(__urng), 1527 -result_type(1) / __p.alpha()); 1528 } 1529 1530 template<typename _ForwardIterator, 1531 typename _UniformRandomNumberGenerator> 1532 void 1533 __generate(_ForwardIterator __f, _ForwardIterator __t, 1534 _UniformRandomNumberGenerator& __urng) 1535 { this->__generate(__f, __t, __urng, _M_param); } 1536 1537 template<typename _ForwardIterator, 1538 typename _UniformRandomNumberGenerator> 1539 void 1540 __generate(_ForwardIterator __f, _ForwardIterator __t, 1541 _UniformRandomNumberGenerator& __urng, 1542 const param_type& __p) 1543 { this->__generate_impl(__f, __t, __urng, __p); } 1544 1545 template<typename _UniformRandomNumberGenerator> 1546 void 1547 __generate(result_type* __f, result_type* __t, 1548 _UniformRandomNumberGenerator& __urng, 1549 const param_type& __p) 1550 { this->__generate_impl(__f, __t, __urng, __p); } 1551 1552 /** 1553 * @brief Return true if two Pareto distributions have 1554 * the same parameters and the sequences that would 1555 * be generated are equal. 1556 */ 1557 friend bool 1558 operator==(const pareto_distribution& __d1, 1559 const pareto_distribution& __d2) 1560 { return (__d1._M_param == __d2._M_param 1561 && __d1._M_ud == __d2._M_ud); } 1562 1563 /** 1564 * @brief Inserts a %pareto_distribution random number distribution 1565 * @p __x into the output stream @p __os. 1566 * 1567 * @param __os An output stream. 1568 * @param __x A %pareto_distribution random number distribution. 1569 * 1570 * @returns The output stream with the state of @p __x inserted or in 1571 * an error state. 1572 */ 1573 template<typename _RealType1, typename _CharT, typename _Traits> 1574 friend std::basic_ostream<_CharT, _Traits>& 1575 operator<<(std::basic_ostream<_CharT, _Traits>&, 1576 const pareto_distribution<_RealType1>&); 1577 1578 /** 1579 * @brief Extracts a %pareto_distribution random number distribution 1580 * @p __x from the input stream @p __is. 1581 * 1582 * @param __is An input stream. 1583 * @param __x A %pareto_distribution random number 1584 * generator engine. 1585 * 1586 * @returns The input stream with @p __x extracted or in an error state. 1587 */ 1588 template<typename _RealType1, typename _CharT, typename _Traits> 1589 friend std::basic_istream<_CharT, _Traits>& 1590 operator>>(std::basic_istream<_CharT, _Traits>&, 1591 pareto_distribution<_RealType1>&); 1592 1593 private: 1594 template<typename _ForwardIterator, 1595 typename _UniformRandomNumberGenerator> 1596 void 1597 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 1598 _UniformRandomNumberGenerator& __urng, 1599 const param_type& __p); 1600 1601 param_type _M_param; 1602 1603 std::uniform_real_distribution<result_type> _M_ud; 1604 }; 1605 1606 /** 1607 * @brief Return true if two Pareto distributions are not equal. 1608 */ 1609 template<typename _RealType> 1610 inline bool 1611 operator!=(const pareto_distribution<_RealType>& __d1, 1612 const pareto_distribution<_RealType>& __d2) 1613 { return !(__d1 == __d2); } 1614 1615 1616 /** 1617 * @brief A K continuous distribution for random numbers. 1618 * 1619 * The formula for the K probability density function is 1620 * @f[ 1621 * p(x|\lambda, \mu, \nu) = \frac{2}{x} 1622 * \left(\frac{\lambda\nu x}{\mu}\right)^{\frac{\lambda + \nu}{2}} 1623 * \frac{1}{\Gamma(\lambda)\Gamma(\nu)} 1624 * K_{\nu - \lambda}\left(2\sqrt{\frac{\lambda\nu x}{\mu}}\right) 1625 * @f] 1626 * where @f$I_0(z)@f$ is the modified Bessel function of the second kind 1627 * of order @f$\nu - \lambda@f$ and @f$\lambda > 0@f$, @f$\mu > 0@f$ 1628 * and @f$\nu > 0@f$. 1629 * 1630 * <table border=1 cellpadding=10 cellspacing=0> 1631 * <caption align=top>Distribution Statistics</caption> 1632 * <tr><td>Mean</td><td>@f$\mu@f$</td></tr> 1633 * <tr><td>Variance</td><td>@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$</td></tr> 1634 * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr> 1635 * </table> 1636 */ 1637 template<typename _RealType = double> 1638 class 1639 k_distribution 1640 { 1641 static_assert(std::is_floating_point<_RealType>::value, 1642 "template argument not a floating point type"); 1643 1644 public: 1645 /** The type of the range of the distribution. */ 1646 typedef _RealType result_type; 1647 /** Parameter type. */ 1648 struct param_type 1649 { 1650 typedef k_distribution<result_type> distribution_type; 1651 1652 param_type(result_type __lambda_val = result_type(1), 1653 result_type __mu_val = result_type(1), 1654 result_type __nu_val = result_type(1)) 1655 : _M_lambda(__lambda_val), _M_mu(__mu_val), _M_nu(__nu_val) 1656 { 1657 _GLIBCXX_DEBUG_ASSERT(_M_lambda > result_type(0)); 1658 _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0)); 1659 _GLIBCXX_DEBUG_ASSERT(_M_nu > result_type(0)); 1660 } 1661 1662 result_type 1663 lambda() const 1664 { return _M_lambda; } 1665 1666 result_type 1667 mu() const 1668 { return _M_mu; } 1669 1670 result_type 1671 nu() const 1672 { return _M_nu; } 1673 1674 friend bool 1675 operator==(const param_type& __p1, const param_type& __p2) 1676 { return __p1._M_lambda == __p2._M_lambda 1677 && __p1._M_mu == __p2._M_mu 1678 && __p1._M_nu == __p2._M_nu; } 1679 1680 private: 1681 void _M_initialize(); 1682 1683 result_type _M_lambda; 1684 result_type _M_mu; 1685 result_type _M_nu; 1686 }; 1687 1688 /** 1689 * @brief Constructors. 1690 */ 1691 explicit 1692 k_distribution(result_type __lambda_val = result_type(1), 1693 result_type __mu_val = result_type(1), 1694 result_type __nu_val = result_type(1)) 1695 : _M_param(__lambda_val, __mu_val, __nu_val), 1696 _M_gd1(__lambda_val, result_type(1) / __lambda_val), 1697 _M_gd2(__nu_val, __mu_val / __nu_val) 1698 { } 1699 1700 explicit 1701 k_distribution(const param_type& __p) 1702 : _M_param(__p), 1703 _M_gd1(__p.lambda(), result_type(1) / __p.lambda()), 1704 _M_gd2(__p.nu(), __p.mu() / __p.nu()) 1705 { } 1706 1707 /** 1708 * @brief Resets the distribution state. 1709 */ 1710 void 1711 reset() 1712 { 1713 _M_gd1.reset(); 1714 _M_gd2.reset(); 1715 } 1716 1717 /** 1718 * @brief Return the parameters of the distribution. 1719 */ 1720 result_type 1721 lambda() const 1722 { return _M_param.lambda(); } 1723 1724 result_type 1725 mu() const 1726 { return _M_param.mu(); } 1727 1728 result_type 1729 nu() const 1730 { return _M_param.nu(); } 1731 1732 /** 1733 * @brief Returns the parameter set of the distribution. 1734 */ 1735 param_type 1736 param() const 1737 { return _M_param; } 1738 1739 /** 1740 * @brief Sets the parameter set of the distribution. 1741 * @param __param The new parameter set of the distribution. 1742 */ 1743 void 1744 param(const param_type& __param) 1745 { _M_param = __param; } 1746 1747 /** 1748 * @brief Returns the greatest lower bound value of the distribution. 1749 */ 1750 result_type 1751 min() const 1752 { return result_type(0); } 1753 1754 /** 1755 * @brief Returns the least upper bound value of the distribution. 1756 */ 1757 result_type 1758 max() const 1759 { return std::numeric_limits<result_type>::max(); } 1760 1761 /** 1762 * @brief Generating functions. 1763 */ 1764 template<typename _UniformRandomNumberGenerator> 1765 result_type 1766 operator()(_UniformRandomNumberGenerator&); 1767 1768 template<typename _UniformRandomNumberGenerator> 1769 result_type 1770 operator()(_UniformRandomNumberGenerator&, const param_type&); 1771 1772 template<typename _ForwardIterator, 1773 typename _UniformRandomNumberGenerator> 1774 void 1775 __generate(_ForwardIterator __f, _ForwardIterator __t, 1776 _UniformRandomNumberGenerator& __urng) 1777 { this->__generate(__f, __t, __urng, _M_param); } 1778 1779 template<typename _ForwardIterator, 1780 typename _UniformRandomNumberGenerator> 1781 void 1782 __generate(_ForwardIterator __f, _ForwardIterator __t, 1783 _UniformRandomNumberGenerator& __urng, 1784 const param_type& __p) 1785 { this->__generate_impl(__f, __t, __urng, __p); } 1786 1787 template<typename _UniformRandomNumberGenerator> 1788 void 1789 __generate(result_type* __f, result_type* __t, 1790 _UniformRandomNumberGenerator& __urng, 1791 const param_type& __p) 1792 { this->__generate_impl(__f, __t, __urng, __p); } 1793 1794 /** 1795 * @brief Return true if two K distributions have 1796 * the same parameters and the sequences that would 1797 * be generated are equal. 1798 */ 1799 friend bool 1800 operator==(const k_distribution& __d1, 1801 const k_distribution& __d2) 1802 { return (__d1._M_param == __d2._M_param 1803 && __d1._M_gd1 == __d2._M_gd1 1804 && __d1._M_gd2 == __d2._M_gd2); } 1805 1806 /** 1807 * @brief Inserts a %k_distribution random number distribution 1808 * @p __x into the output stream @p __os. 1809 * 1810 * @param __os An output stream. 1811 * @param __x A %k_distribution random number distribution. 1812 * 1813 * @returns The output stream with the state of @p __x inserted or in 1814 * an error state. 1815 */ 1816 template<typename _RealType1, typename _CharT, typename _Traits> 1817 friend std::basic_ostream<_CharT, _Traits>& 1818 operator<<(std::basic_ostream<_CharT, _Traits>&, 1819 const k_distribution<_RealType1>&); 1820 1821 /** 1822 * @brief Extracts a %k_distribution random number distribution 1823 * @p __x from the input stream @p __is. 1824 * 1825 * @param __is An input stream. 1826 * @param __x A %k_distribution random number 1827 * generator engine. 1828 * 1829 * @returns The input stream with @p __x extracted or in an error state. 1830 */ 1831 template<typename _RealType1, typename _CharT, typename _Traits> 1832 friend std::basic_istream<_CharT, _Traits>& 1833 operator>>(std::basic_istream<_CharT, _Traits>&, 1834 k_distribution<_RealType1>&); 1835 1836 private: 1837 template<typename _ForwardIterator, 1838 typename _UniformRandomNumberGenerator> 1839 void 1840 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 1841 _UniformRandomNumberGenerator& __urng, 1842 const param_type& __p); 1843 1844 param_type _M_param; 1845 1846 std::gamma_distribution<result_type> _M_gd1; 1847 std::gamma_distribution<result_type> _M_gd2; 1848 }; 1849 1850 /** 1851 * @brief Return true if two K distributions are not equal. 1852 */ 1853 template<typename _RealType> 1854 inline bool 1855 operator!=(const k_distribution<_RealType>& __d1, 1856 const k_distribution<_RealType>& __d2) 1857 { return !(__d1 == __d2); } 1858 1859 1860 /** 1861 * @brief An arcsine continuous distribution for random numbers. 1862 * 1863 * The formula for the arcsine probability density function is 1864 * @f[ 1865 * p(x|a,b) = \frac{1}{\pi \sqrt{(x - a)(b - x)}} 1866 * @f] 1867 * where @f$x >= a@f$ and @f$x <= b@f$. 1868 * 1869 * <table border=1 cellpadding=10 cellspacing=0> 1870 * <caption align=top>Distribution Statistics</caption> 1871 * <tr><td>Mean</td><td>@f$ (a + b) / 2 @f$</td></tr> 1872 * <tr><td>Variance</td><td>@f$ (b - a)^2 / 8 @f$</td></tr> 1873 * <tr><td>Range</td><td>@f$[a, b]@f$</td></tr> 1874 * </table> 1875 */ 1876 template<typename _RealType = double> 1877 class 1878 arcsine_distribution 1879 { 1880 static_assert(std::is_floating_point<_RealType>::value, 1881 "template argument not a floating point type"); 1882 1883 public: 1884 /** The type of the range of the distribution. */ 1885 typedef _RealType result_type; 1886 /** Parameter type. */ 1887 struct param_type 1888 { 1889 typedef arcsine_distribution<result_type> distribution_type; 1890 1891 param_type(result_type __a = result_type(0), 1892 result_type __b = result_type(1)) 1893 : _M_a(__a), _M_b(__b) 1894 { 1895 _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b); 1896 } 1897 1898 result_type 1899 a() const 1900 { return _M_a; } 1901 1902 result_type 1903 b() const 1904 { return _M_b; } 1905 1906 friend bool 1907 operator==(const param_type& __p1, const param_type& __p2) 1908 { return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; } 1909 1910 private: 1911 void _M_initialize(); 1912 1913 result_type _M_a; 1914 result_type _M_b; 1915 }; 1916 1917 /** 1918 * @brief Constructors. 1919 */ 1920 explicit 1921 arcsine_distribution(result_type __a = result_type(0), 1922 result_type __b = result_type(1)) 1923 : _M_param(__a, __b), 1924 _M_ud(-1.5707963267948966192313216916397514L, 1925 +1.5707963267948966192313216916397514L) 1926 { } 1927 1928 explicit 1929 arcsine_distribution(const param_type& __p) 1930 : _M_param(__p), 1931 _M_ud(-1.5707963267948966192313216916397514L, 1932 +1.5707963267948966192313216916397514L) 1933 { } 1934 1935 /** 1936 * @brief Resets the distribution state. 1937 */ 1938 void 1939 reset() 1940 { _M_ud.reset(); } 1941 1942 /** 1943 * @brief Return the parameters of the distribution. 1944 */ 1945 result_type 1946 a() const 1947 { return _M_param.a(); } 1948 1949 result_type 1950 b() const 1951 { return _M_param.b(); } 1952 1953 /** 1954 * @brief Returns the parameter set of the distribution. 1955 */ 1956 param_type 1957 param() const 1958 { return _M_param; } 1959 1960 /** 1961 * @brief Sets the parameter set of the distribution. 1962 * @param __param The new parameter set of the distribution. 1963 */ 1964 void 1965 param(const param_type& __param) 1966 { _M_param = __param; } 1967 1968 /** 1969 * @brief Returns the greatest lower bound value of the distribution. 1970 */ 1971 result_type 1972 min() const 1973 { return this->a(); } 1974 1975 /** 1976 * @brief Returns the least upper bound value of the distribution. 1977 */ 1978 result_type 1979 max() const 1980 { return this->b(); } 1981 1982 /** 1983 * @brief Generating functions. 1984 */ 1985 template<typename _UniformRandomNumberGenerator> 1986 result_type 1987 operator()(_UniformRandomNumberGenerator& __urng) 1988 { 1989 result_type __x = std::sin(this->_M_ud(__urng)); 1990 return (__x * (this->b() - this->a()) 1991 + this->a() + this->b()) / result_type(2); 1992 } 1993 1994 template<typename _UniformRandomNumberGenerator> 1995 result_type 1996 operator()(_UniformRandomNumberGenerator& __urng, 1997 const param_type& __p) 1998 { 1999 result_type __x = std::sin(this->_M_ud(__urng)); 2000 return (__x * (__p.b() - __p.a()) 2001 + __p.a() + __p.b()) / result_type(2); 2002 } 2003 2004 template<typename _ForwardIterator, 2005 typename _UniformRandomNumberGenerator> 2006 void 2007 __generate(_ForwardIterator __f, _ForwardIterator __t, 2008 _UniformRandomNumberGenerator& __urng) 2009 { this->__generate(__f, __t, __urng, _M_param); } 2010 2011 template<typename _ForwardIterator, 2012 typename _UniformRandomNumberGenerator> 2013 void 2014 __generate(_ForwardIterator __f, _ForwardIterator __t, 2015 _UniformRandomNumberGenerator& __urng, 2016 const param_type& __p) 2017 { this->__generate_impl(__f, __t, __urng, __p); } 2018 2019 template<typename _UniformRandomNumberGenerator> 2020 void 2021 __generate(result_type* __f, result_type* __t, 2022 _UniformRandomNumberGenerator& __urng, 2023 const param_type& __p) 2024 { this->__generate_impl(__f, __t, __urng, __p); } 2025 2026 /** 2027 * @brief Return true if two arcsine distributions have 2028 * the same parameters and the sequences that would 2029 * be generated are equal. 2030 */ 2031 friend bool 2032 operator==(const arcsine_distribution& __d1, 2033 const arcsine_distribution& __d2) 2034 { return (__d1._M_param == __d2._M_param 2035 && __d1._M_ud == __d2._M_ud); } 2036 2037 /** 2038 * @brief Inserts a %arcsine_distribution random number distribution 2039 * @p __x into the output stream @p __os. 2040 * 2041 * @param __os An output stream. 2042 * @param __x A %arcsine_distribution random number distribution. 2043 * 2044 * @returns The output stream with the state of @p __x inserted or in 2045 * an error state. 2046 */ 2047 template<typename _RealType1, typename _CharT, typename _Traits> 2048 friend std::basic_ostream<_CharT, _Traits>& 2049 operator<<(std::basic_ostream<_CharT, _Traits>&, 2050 const arcsine_distribution<_RealType1>&); 2051 2052 /** 2053 * @brief Extracts a %arcsine_distribution random number distribution 2054 * @p __x from the input stream @p __is. 2055 * 2056 * @param __is An input stream. 2057 * @param __x A %arcsine_distribution random number 2058 * generator engine. 2059 * 2060 * @returns The input stream with @p __x extracted or in an error state. 2061 */ 2062 template<typename _RealType1, typename _CharT, typename _Traits> 2063 friend std::basic_istream<_CharT, _Traits>& 2064 operator>>(std::basic_istream<_CharT, _Traits>&, 2065 arcsine_distribution<_RealType1>&); 2066 2067 private: 2068 template<typename _ForwardIterator, 2069 typename _UniformRandomNumberGenerator> 2070 void 2071 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2072 _UniformRandomNumberGenerator& __urng, 2073 const param_type& __p); 2074 2075 param_type _M_param; 2076 2077 std::uniform_real_distribution<result_type> _M_ud; 2078 }; 2079 2080 /** 2081 * @brief Return true if two arcsine distributions are not equal. 2082 */ 2083 template<typename _RealType> 2084 inline bool 2085 operator!=(const arcsine_distribution<_RealType>& __d1, 2086 const arcsine_distribution<_RealType>& __d2) 2087 { return !(__d1 == __d2); } 2088 2089 2090 /** 2091 * @brief A Hoyt continuous distribution for random numbers. 2092 * 2093 * The formula for the Hoyt probability density function is 2094 * @f[ 2095 * p(x|q,\omega) = \frac{(1 + q^2)x}{q\omega} 2096 * \exp\left(-\frac{(1 + q^2)^2 x^2}{4 q^2 \omega}\right) 2097 * I_0\left(\frac{(1 - q^4) x^2}{4 q^2 \omega}\right) 2098 * @f] 2099 * where @f$I_0(z)@f$ is the modified Bessel function of the first kind 2100 * of order 0 and @f$0 < q < 1@f$. 2101 * 2102 * <table border=1 cellpadding=10 cellspacing=0> 2103 * <caption align=top>Distribution Statistics</caption> 2104 * <tr><td>Mean</td><td>@f$ \sqrt{\frac{2}{\pi}} \sqrt{\frac{\omega}{1 + q^2}} 2105 * E(1 - q^2) @f$</td></tr> 2106 * <tr><td>Variance</td><td>@f$ \omega \left(1 - \frac{2E^2(1 - q^2)} 2107 * {\pi (1 + q^2)}\right) @f$</td></tr> 2108 * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr> 2109 * </table> 2110 * where @f$E(x)@f$ is the elliptic function of the second kind. 2111 */ 2112 template<typename _RealType = double> 2113 class 2114 hoyt_distribution 2115 { 2116 static_assert(std::is_floating_point<_RealType>::value, 2117 "template argument not a floating point type"); 2118 2119 public: 2120 /** The type of the range of the distribution. */ 2121 typedef _RealType result_type; 2122 /** Parameter type. */ 2123 struct param_type 2124 { 2125 typedef hoyt_distribution<result_type> distribution_type; 2126 2127 param_type(result_type __q = result_type(0.5L), 2128 result_type __omega = result_type(1)) 2129 : _M_q(__q), _M_omega(__omega) 2130 { 2131 _GLIBCXX_DEBUG_ASSERT(_M_q > result_type(0)); 2132 _GLIBCXX_DEBUG_ASSERT(_M_q < result_type(1)); 2133 } 2134 2135 result_type 2136 q() const 2137 { return _M_q; } 2138 2139 result_type 2140 omega() const 2141 { return _M_omega; } 2142 2143 friend bool 2144 operator==(const param_type& __p1, const param_type& __p2) 2145 { return __p1._M_q == __p2._M_q 2146 && __p1._M_omega == __p2._M_omega; } 2147 2148 private: 2149 void _M_initialize(); 2150 2151 result_type _M_q; 2152 result_type _M_omega; 2153 }; 2154 2155 /** 2156 * @brief Constructors. 2157 */ 2158 explicit 2159 hoyt_distribution(result_type __q = result_type(0.5L), 2160 result_type __omega = result_type(1)) 2161 : _M_param(__q, __omega), 2162 _M_ad(result_type(0.5L) * (result_type(1) + __q * __q), 2163 result_type(0.5L) * (result_type(1) + __q * __q) 2164 / (__q * __q)), 2165 _M_ed(result_type(1)) 2166 { } 2167 2168 explicit 2169 hoyt_distribution(const param_type& __p) 2170 : _M_param(__p), 2171 _M_ad(result_type(0.5L) * (result_type(1) + __p.q() * __p.q()), 2172 result_type(0.5L) * (result_type(1) + __p.q() * __p.q()) 2173 / (__p.q() * __p.q())), 2174 _M_ed(result_type(1)) 2175 { } 2176 2177 /** 2178 * @brief Resets the distribution state. 2179 */ 2180 void 2181 reset() 2182 { 2183 _M_ad.reset(); 2184 _M_ed.reset(); 2185 } 2186 2187 /** 2188 * @brief Return the parameters of the distribution. 2189 */ 2190 result_type 2191 q() const 2192 { return _M_param.q(); } 2193 2194 result_type 2195 omega() const 2196 { return _M_param.omega(); } 2197 2198 /** 2199 * @brief Returns the parameter set of the distribution. 2200 */ 2201 param_type 2202 param() const 2203 { return _M_param; } 2204 2205 /** 2206 * @brief Sets the parameter set of the distribution. 2207 * @param __param The new parameter set of the distribution. 2208 */ 2209 void 2210 param(const param_type& __param) 2211 { _M_param = __param; } 2212 2213 /** 2214 * @brief Returns the greatest lower bound value of the distribution. 2215 */ 2216 result_type 2217 min() const 2218 { return result_type(0); } 2219 2220 /** 2221 * @brief Returns the least upper bound value of the distribution. 2222 */ 2223 result_type 2224 max() const 2225 { return std::numeric_limits<result_type>::max(); } 2226 2227 /** 2228 * @brief Generating functions. 2229 */ 2230 template<typename _UniformRandomNumberGenerator> 2231 result_type 2232 operator()(_UniformRandomNumberGenerator& __urng); 2233 2234 template<typename _UniformRandomNumberGenerator> 2235 result_type 2236 operator()(_UniformRandomNumberGenerator& __urng, 2237 const param_type& __p); 2238 2239 template<typename _ForwardIterator, 2240 typename _UniformRandomNumberGenerator> 2241 void 2242 __generate(_ForwardIterator __f, _ForwardIterator __t, 2243 _UniformRandomNumberGenerator& __urng) 2244 { this->__generate(__f, __t, __urng, _M_param); } 2245 2246 template<typename _ForwardIterator, 2247 typename _UniformRandomNumberGenerator> 2248 void 2249 __generate(_ForwardIterator __f, _ForwardIterator __t, 2250 _UniformRandomNumberGenerator& __urng, 2251 const param_type& __p) 2252 { this->__generate_impl(__f, __t, __urng, __p); } 2253 2254 template<typename _UniformRandomNumberGenerator> 2255 void 2256 __generate(result_type* __f, result_type* __t, 2257 _UniformRandomNumberGenerator& __urng, 2258 const param_type& __p) 2259 { this->__generate_impl(__f, __t, __urng, __p); } 2260 2261 /** 2262 * @brief Return true if two Hoyt distributions have 2263 * the same parameters and the sequences that would 2264 * be generated are equal. 2265 */ 2266 friend bool 2267 operator==(const hoyt_distribution& __d1, 2268 const hoyt_distribution& __d2) 2269 { return (__d1._M_param == __d2._M_param 2270 && __d1._M_ad == __d2._M_ad 2271 && __d1._M_ed == __d2._M_ed); } 2272 2273 /** 2274 * @brief Inserts a %hoyt_distribution random number distribution 2275 * @p __x into the output stream @p __os. 2276 * 2277 * @param __os An output stream. 2278 * @param __x A %hoyt_distribution random number distribution. 2279 * 2280 * @returns The output stream with the state of @p __x inserted or in 2281 * an error state. 2282 */ 2283 template<typename _RealType1, typename _CharT, typename _Traits> 2284 friend std::basic_ostream<_CharT, _Traits>& 2285 operator<<(std::basic_ostream<_CharT, _Traits>&, 2286 const hoyt_distribution<_RealType1>&); 2287 2288 /** 2289 * @brief Extracts a %hoyt_distribution random number distribution 2290 * @p __x from the input stream @p __is. 2291 * 2292 * @param __is An input stream. 2293 * @param __x A %hoyt_distribution random number 2294 * generator engine. 2295 * 2296 * @returns The input stream with @p __x extracted or in an error state. 2297 */ 2298 template<typename _RealType1, typename _CharT, typename _Traits> 2299 friend std::basic_istream<_CharT, _Traits>& 2300 operator>>(std::basic_istream<_CharT, _Traits>&, 2301 hoyt_distribution<_RealType1>&); 2302 2303 private: 2304 template<typename _ForwardIterator, 2305 typename _UniformRandomNumberGenerator> 2306 void 2307 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2308 _UniformRandomNumberGenerator& __urng, 2309 const param_type& __p); 2310 2311 param_type _M_param; 2312 2313 __gnu_cxx::arcsine_distribution<result_type> _M_ad; 2314 std::exponential_distribution<result_type> _M_ed; 2315 }; 2316 2317 /** 2318 * @brief Return true if two Hoyt distributions are not equal. 2319 */ 2320 template<typename _RealType> 2321 inline bool 2322 operator!=(const hoyt_distribution<_RealType>& __d1, 2323 const hoyt_distribution<_RealType>& __d2) 2324 { return !(__d1 == __d2); } 2325 2326 2327 /** 2328 * @brief A triangular distribution for random numbers. 2329 * 2330 * The formula for the triangular probability density function is 2331 * @f[ 2332 * / 0 for x < a 2333 * p(x|a,b,c) = | \frac{2(x-a)}{(c-a)(b-a)} for a <= x <= b 2334 * | \frac{2(c-x)}{(c-a)(c-b)} for b < x <= c 2335 * \ 0 for c < x 2336 * @f] 2337 * 2338 * <table border=1 cellpadding=10 cellspacing=0> 2339 * <caption align=top>Distribution Statistics</caption> 2340 * <tr><td>Mean</td><td>@f$ \frac{a+b+c}{2} @f$</td></tr> 2341 * <tr><td>Variance</td><td>@f$ \frac{a^2+b^2+c^2-ab-ac-bc} 2342 * {18}@f$</td></tr> 2343 * <tr><td>Range</td><td>@f$[a, c]@f$</td></tr> 2344 * </table> 2345 */ 2346 template<typename _RealType = double> 2347 class triangular_distribution 2348 { 2349 static_assert(std::is_floating_point<_RealType>::value, 2350 "template argument not a floating point type"); 2351 2352 public: 2353 /** The type of the range of the distribution. */ 2354 typedef _RealType result_type; 2355 /** Parameter type. */ 2356 struct param_type 2357 { 2358 friend class triangular_distribution<_RealType>; 2359 2360 explicit 2361 param_type(_RealType __a = _RealType(0), 2362 _RealType __b = _RealType(0.5), 2363 _RealType __c = _RealType(1)) 2364 : _M_a(__a), _M_b(__b), _M_c(__c) 2365 { 2366 _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b); 2367 _GLIBCXX_DEBUG_ASSERT(_M_b <= _M_c); 2368 _GLIBCXX_DEBUG_ASSERT(_M_a < _M_c); 2369 2370 _M_r_ab = (_M_b - _M_a) / (_M_c - _M_a); 2371 _M_f_ab_ac = (_M_b - _M_a) * (_M_c - _M_a); 2372 _M_f_bc_ac = (_M_c - _M_b) * (_M_c - _M_a); 2373 } 2374 2375 _RealType 2376 a() const 2377 { return _M_a; } 2378 2379 _RealType 2380 b() const 2381 { return _M_b; } 2382 2383 _RealType 2384 c() const 2385 { return _M_c; } 2386 2387 friend bool 2388 operator==(const param_type& __p1, const param_type& __p2) 2389 { return (__p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b 2390 && __p1._M_c == __p2._M_c); } 2391 2392 private: 2393 2394 _RealType _M_a; 2395 _RealType _M_b; 2396 _RealType _M_c; 2397 _RealType _M_r_ab; 2398 _RealType _M_f_ab_ac; 2399 _RealType _M_f_bc_ac; 2400 }; 2401 2402 /** 2403 * @brief Constructs a triangle distribution with parameters 2404 * @f$ a @f$, @f$ b @f$ and @f$ c @f$. 2405 */ 2406 explicit 2407 triangular_distribution(result_type __a = result_type(0), 2408 result_type __b = result_type(0.5), 2409 result_type __c = result_type(1)) 2410 : _M_param(__a, __b, __c) 2411 { } 2412 2413 explicit 2414 triangular_distribution(const param_type& __p) 2415 : _M_param(__p) 2416 { } 2417 2418 /** 2419 * @brief Resets the distribution state. 2420 */ 2421 void 2422 reset() 2423 { } 2424 2425 /** 2426 * @brief Returns the @f$ a @f$ of the distribution. 2427 */ 2428 result_type 2429 a() const 2430 { return _M_param.a(); } 2431 2432 /** 2433 * @brief Returns the @f$ b @f$ of the distribution. 2434 */ 2435 result_type 2436 b() const 2437 { return _M_param.b(); } 2438 2439 /** 2440 * @brief Returns the @f$ c @f$ of the distribution. 2441 */ 2442 result_type 2443 c() const 2444 { return _M_param.c(); } 2445 2446 /** 2447 * @brief Returns the parameter set of the distribution. 2448 */ 2449 param_type 2450 param() const 2451 { return _M_param; } 2452 2453 /** 2454 * @brief Sets the parameter set of the distribution. 2455 * @param __param The new parameter set of the distribution. 2456 */ 2457 void 2458 param(const param_type& __param) 2459 { _M_param = __param; } 2460 2461 /** 2462 * @brief Returns the greatest lower bound value of the distribution. 2463 */ 2464 result_type 2465 min() const 2466 { return _M_param._M_a; } 2467 2468 /** 2469 * @brief Returns the least upper bound value of the distribution. 2470 */ 2471 result_type 2472 max() const 2473 { return _M_param._M_c; } 2474 2475 /** 2476 * @brief Generating functions. 2477 */ 2478 template<typename _UniformRandomNumberGenerator> 2479 result_type 2480 operator()(_UniformRandomNumberGenerator& __urng) 2481 { return this->operator()(__urng, _M_param); } 2482 2483 template<typename _UniformRandomNumberGenerator> 2484 result_type 2485 operator()(_UniformRandomNumberGenerator& __urng, 2486 const param_type& __p) 2487 { 2488 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2489 __aurng(__urng); 2490 result_type __rnd = __aurng(); 2491 if (__rnd <= __p._M_r_ab) 2492 return __p.a() + std::sqrt(__rnd * __p._M_f_ab_ac); 2493 else 2494 return __p.c() - std::sqrt((result_type(1) - __rnd) 2495 * __p._M_f_bc_ac); 2496 } 2497 2498 template<typename _ForwardIterator, 2499 typename _UniformRandomNumberGenerator> 2500 void 2501 __generate(_ForwardIterator __f, _ForwardIterator __t, 2502 _UniformRandomNumberGenerator& __urng) 2503 { this->__generate(__f, __t, __urng, _M_param); } 2504 2505 template<typename _ForwardIterator, 2506 typename _UniformRandomNumberGenerator> 2507 void 2508 __generate(_ForwardIterator __f, _ForwardIterator __t, 2509 _UniformRandomNumberGenerator& __urng, 2510 const param_type& __p) 2511 { this->__generate_impl(__f, __t, __urng, __p); } 2512 2513 template<typename _UniformRandomNumberGenerator> 2514 void 2515 __generate(result_type* __f, result_type* __t, 2516 _UniformRandomNumberGenerator& __urng, 2517 const param_type& __p) 2518 { this->__generate_impl(__f, __t, __urng, __p); } 2519 2520 /** 2521 * @brief Return true if two triangle distributions have the same 2522 * parameters and the sequences that would be generated 2523 * are equal. 2524 */ 2525 friend bool 2526 operator==(const triangular_distribution& __d1, 2527 const triangular_distribution& __d2) 2528 { return __d1._M_param == __d2._M_param; } 2529 2530 /** 2531 * @brief Inserts a %triangular_distribution random number distribution 2532 * @p __x into the output stream @p __os. 2533 * 2534 * @param __os An output stream. 2535 * @param __x A %triangular_distribution random number distribution. 2536 * 2537 * @returns The output stream with the state of @p __x inserted or in 2538 * an error state. 2539 */ 2540 template<typename _RealType1, typename _CharT, typename _Traits> 2541 friend std::basic_ostream<_CharT, _Traits>& 2542 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2543 const __gnu_cxx::triangular_distribution<_RealType1>& __x); 2544 2545 /** 2546 * @brief Extracts a %triangular_distribution random number distribution 2547 * @p __x from the input stream @p __is. 2548 * 2549 * @param __is An input stream. 2550 * @param __x A %triangular_distribution random number generator engine. 2551 * 2552 * @returns The input stream with @p __x extracted or in an error state. 2553 */ 2554 template<typename _RealType1, typename _CharT, typename _Traits> 2555 friend std::basic_istream<_CharT, _Traits>& 2556 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2557 __gnu_cxx::triangular_distribution<_RealType1>& __x); 2558 2559 private: 2560 template<typename _ForwardIterator, 2561 typename _UniformRandomNumberGenerator> 2562 void 2563 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2564 _UniformRandomNumberGenerator& __urng, 2565 const param_type& __p); 2566 2567 param_type _M_param; 2568 }; 2569 2570 /** 2571 * @brief Return true if two triangle distributions are different. 2572 */ 2573 template<typename _RealType> 2574 inline bool 2575 operator!=(const __gnu_cxx::triangular_distribution<_RealType>& __d1, 2576 const __gnu_cxx::triangular_distribution<_RealType>& __d2) 2577 { return !(__d1 == __d2); } 2578 2579 2580 /** 2581 * @brief A von Mises distribution for random numbers. 2582 * 2583 * The formula for the von Mises probability density function is 2584 * @f[ 2585 * p(x|\mu,\kappa) = \frac{e^{\kappa \cos(x-\mu)}} 2586 * {2\pi I_0(\kappa)} 2587 * @f] 2588 * 2589 * The generating functions use the method according to: 2590 * 2591 * D. J. Best and N. I. Fisher, 1979. "Efficient Simulation of the 2592 * von Mises Distribution", Journal of the Royal Statistical Society. 2593 * Series C (Applied Statistics), Vol. 28, No. 2, pp. 152-157. 2594 * 2595 * <table border=1 cellpadding=10 cellspacing=0> 2596 * <caption align=top>Distribution Statistics</caption> 2597 * <tr><td>Mean</td><td>@f$ \mu @f$</td></tr> 2598 * <tr><td>Variance</td><td>@f$ 1-I_1(\kappa)/I_0(\kappa) @f$</td></tr> 2599 * <tr><td>Range</td><td>@f$[-\pi, \pi]@f$</td></tr> 2600 * </table> 2601 */ 2602 template<typename _RealType = double> 2603 class von_mises_distribution 2604 { 2605 static_assert(std::is_floating_point<_RealType>::value, 2606 "template argument not a floating point type"); 2607 2608 public: 2609 /** The type of the range of the distribution. */ 2610 typedef _RealType result_type; 2611 /** Parameter type. */ 2612 struct param_type 2613 { 2614 friend class von_mises_distribution<_RealType>; 2615 2616 explicit 2617 param_type(_RealType __mu = _RealType(0), 2618 _RealType __kappa = _RealType(1)) 2619 : _M_mu(__mu), _M_kappa(__kappa) 2620 { 2621 const _RealType __pi = __gnu_cxx::__math_constants<_RealType>::__pi; 2622 _GLIBCXX_DEBUG_ASSERT(_M_mu >= -__pi && _M_mu <= __pi); 2623 _GLIBCXX_DEBUG_ASSERT(_M_kappa >= _RealType(0)); 2624 2625 auto __tau = std::sqrt(_RealType(4) * _M_kappa * _M_kappa 2626 + _RealType(1)) + _RealType(1); 2627 auto __rho = ((__tau - std::sqrt(_RealType(2) * __tau)) 2628 / (_RealType(2) * _M_kappa)); 2629 _M_r = (_RealType(1) + __rho * __rho) / (_RealType(2) * __rho); 2630 } 2631 2632 _RealType 2633 mu() const 2634 { return _M_mu; } 2635 2636 _RealType 2637 kappa() const 2638 { return _M_kappa; } 2639 2640 friend bool 2641 operator==(const param_type& __p1, const param_type& __p2) 2642 { return (__p1._M_mu == __p2._M_mu 2643 && __p1._M_kappa == __p2._M_kappa); } 2644 2645 private: 2646 _RealType _M_mu; 2647 _RealType _M_kappa; 2648 _RealType _M_r; 2649 }; 2650 2651 /** 2652 * @brief Constructs a von Mises distribution with parameters 2653 * @f$\mu@f$ and @f$\kappa@f$. 2654 */ 2655 explicit 2656 von_mises_distribution(result_type __mu = result_type(0), 2657 result_type __kappa = result_type(1)) 2658 : _M_param(__mu, __kappa) 2659 { } 2660 2661 explicit 2662 von_mises_distribution(const param_type& __p) 2663 : _M_param(__p) 2664 { } 2665 2666 /** 2667 * @brief Resets the distribution state. 2668 */ 2669 void 2670 reset() 2671 { } 2672 2673 /** 2674 * @brief Returns the @f$ \mu @f$ of the distribution. 2675 */ 2676 result_type 2677 mu() const 2678 { return _M_param.mu(); } 2679 2680 /** 2681 * @brief Returns the @f$ \kappa @f$ of the distribution. 2682 */ 2683 result_type 2684 kappa() const 2685 { return _M_param.kappa(); } 2686 2687 /** 2688 * @brief Returns the parameter set of the distribution. 2689 */ 2690 param_type 2691 param() const 2692 { return _M_param; } 2693 2694 /** 2695 * @brief Sets the parameter set of the distribution. 2696 * @param __param The new parameter set of the distribution. 2697 */ 2698 void 2699 param(const param_type& __param) 2700 { _M_param = __param; } 2701 2702 /** 2703 * @brief Returns the greatest lower bound value of the distribution. 2704 */ 2705 result_type 2706 min() const 2707 { 2708 return -__gnu_cxx::__math_constants<result_type>::__pi; 2709 } 2710 2711 /** 2712 * @brief Returns the least upper bound value of the distribution. 2713 */ 2714 result_type 2715 max() const 2716 { 2717 return __gnu_cxx::__math_constants<result_type>::__pi; 2718 } 2719 2720 /** 2721 * @brief Generating functions. 2722 */ 2723 template<typename _UniformRandomNumberGenerator> 2724 result_type 2725 operator()(_UniformRandomNumberGenerator& __urng) 2726 { return this->operator()(__urng, _M_param); } 2727 2728 template<typename _UniformRandomNumberGenerator> 2729 result_type 2730 operator()(_UniformRandomNumberGenerator& __urng, 2731 const param_type& __p) 2732 { 2733 const result_type __pi 2734 = __gnu_cxx::__math_constants<result_type>::__pi; 2735 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2736 __aurng(__urng); 2737 2738 result_type __f; 2739 while (1) 2740 { 2741 result_type __rnd = std::cos(__pi * __aurng()); 2742 __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd); 2743 result_type __c = __p._M_kappa * (__p._M_r - __f); 2744 2745 result_type __rnd2 = __aurng(); 2746 if (__c * (result_type(2) - __c) > __rnd2) 2747 break; 2748 if (std::log(__c / __rnd2) >= __c - result_type(1)) 2749 break; 2750 } 2751 2752 result_type __res = std::acos(__f); 2753#if _GLIBCXX_USE_C99_MATH_TR1 2754 __res = std::copysign(__res, __aurng() - result_type(0.5)); 2755#else 2756 if (__aurng() < result_type(0.5)) 2757 __res = -__res; 2758#endif 2759 __res += __p._M_mu; 2760 if (__res > __pi) 2761 __res -= result_type(2) * __pi; 2762 else if (__res < -__pi) 2763 __res += result_type(2) * __pi; 2764 return __res; 2765 } 2766 2767 template<typename _ForwardIterator, 2768 typename _UniformRandomNumberGenerator> 2769 void 2770 __generate(_ForwardIterator __f, _ForwardIterator __t, 2771 _UniformRandomNumberGenerator& __urng) 2772 { this->__generate(__f, __t, __urng, _M_param); } 2773 2774 template<typename _ForwardIterator, 2775 typename _UniformRandomNumberGenerator> 2776 void 2777 __generate(_ForwardIterator __f, _ForwardIterator __t, 2778 _UniformRandomNumberGenerator& __urng, 2779 const param_type& __p) 2780 { this->__generate_impl(__f, __t, __urng, __p); } 2781 2782 template<typename _UniformRandomNumberGenerator> 2783 void 2784 __generate(result_type* __f, result_type* __t, 2785 _UniformRandomNumberGenerator& __urng, 2786 const param_type& __p) 2787 { this->__generate_impl(__f, __t, __urng, __p); } 2788 2789 /** 2790 * @brief Return true if two von Mises distributions have the same 2791 * parameters and the sequences that would be generated 2792 * are equal. 2793 */ 2794 friend bool 2795 operator==(const von_mises_distribution& __d1, 2796 const von_mises_distribution& __d2) 2797 { return __d1._M_param == __d2._M_param; } 2798 2799 /** 2800 * @brief Inserts a %von_mises_distribution random number distribution 2801 * @p __x into the output stream @p __os. 2802 * 2803 * @param __os An output stream. 2804 * @param __x A %von_mises_distribution random number distribution. 2805 * 2806 * @returns The output stream with the state of @p __x inserted or in 2807 * an error state. 2808 */ 2809 template<typename _RealType1, typename _CharT, typename _Traits> 2810 friend std::basic_ostream<_CharT, _Traits>& 2811 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2812 const __gnu_cxx::von_mises_distribution<_RealType1>& __x); 2813 2814 /** 2815 * @brief Extracts a %von_mises_distribution random number distribution 2816 * @p __x from the input stream @p __is. 2817 * 2818 * @param __is An input stream. 2819 * @param __x A %von_mises_distribution random number generator engine. 2820 * 2821 * @returns The input stream with @p __x extracted or in an error state. 2822 */ 2823 template<typename _RealType1, typename _CharT, typename _Traits> 2824 friend std::basic_istream<_CharT, _Traits>& 2825 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2826 __gnu_cxx::von_mises_distribution<_RealType1>& __x); 2827 2828 private: 2829 template<typename _ForwardIterator, 2830 typename _UniformRandomNumberGenerator> 2831 void 2832 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2833 _UniformRandomNumberGenerator& __urng, 2834 const param_type& __p); 2835 2836 param_type _M_param; 2837 }; 2838 2839 /** 2840 * @brief Return true if two von Mises distributions are different. 2841 */ 2842 template<typename _RealType> 2843 inline bool 2844 operator!=(const __gnu_cxx::von_mises_distribution<_RealType>& __d1, 2845 const __gnu_cxx::von_mises_distribution<_RealType>& __d2) 2846 { return !(__d1 == __d2); } 2847 2848 2849 /** 2850 * @brief A discrete hypergeometric random number distribution. 2851 * 2852 * The hypergeometric distribution is a discrete probability distribution 2853 * that describes the probability of @p k successes in @p n draws @a without 2854 * replacement from a finite population of size @p N containing exactly @p K 2855 * successes. 2856 * 2857 * The formula for the hypergeometric probability density function is 2858 * @f[ 2859 * p(k|N,K,n) = \frac{\binom{K}{k} \binom{N-K}{n-k}}{\binom{N}{n}} 2860 * @f] 2861 * where @f$N@f$ is the total population of the distribution, 2862 * @f$K@f$ is the total population of the distribution. 2863 * 2864 * <table border=1 cellpadding=10 cellspacing=0> 2865 * <caption align=top>Distribution Statistics</caption> 2866 * <tr><td>Mean</td><td>@f$ n\frac{K}{N} @f$</td></tr> 2867 * <tr><td>Variance</td><td>@f$ n\frac{K}{N}\frac{N-K}{N}\frac{N-n}{N-1} 2868 * @f$</td></tr> 2869 * <tr><td>Range</td><td>@f$[max(0, n+K-N), min(K, n)]@f$</td></tr> 2870 * </table> 2871 */ 2872 template<typename _UIntType = unsigned int> 2873 class hypergeometric_distribution 2874 { 2875 static_assert(std::is_unsigned<_UIntType>::value, "template argument " 2876 "substituting _UIntType not an unsigned integral type"); 2877 2878 public: 2879 /** The type of the range of the distribution. */ 2880 typedef _UIntType result_type; 2881 2882 /** Parameter type. */ 2883 struct param_type 2884 { 2885 typedef hypergeometric_distribution<_UIntType> distribution_type; 2886 friend class hypergeometric_distribution<_UIntType>; 2887 2888 explicit 2889 param_type(result_type __N = 10, result_type __K = 5, 2890 result_type __n = 1) 2891 : _M_N{__N}, _M_K{__K}, _M_n{__n} 2892 { 2893 _GLIBCXX_DEBUG_ASSERT(_M_N >= _M_K); 2894 _GLIBCXX_DEBUG_ASSERT(_M_N >= _M_n); 2895 } 2896 2897 result_type 2898 total_size() const 2899 { return _M_N; } 2900 2901 result_type 2902 successful_size() const 2903 { return _M_K; } 2904 2905 result_type 2906 unsuccessful_size() const 2907 { return _M_N - _M_K; } 2908 2909 result_type 2910 total_draws() const 2911 { return _M_n; } 2912 2913 friend bool 2914 operator==(const param_type& __p1, const param_type& __p2) 2915 { return (__p1._M_N == __p2._M_N) 2916 && (__p1._M_K == __p2._M_K) 2917 && (__p1._M_n == __p2._M_n); } 2918 2919 private: 2920 2921 result_type _M_N; 2922 result_type _M_K; 2923 result_type _M_n; 2924 }; 2925 2926 // constructors and member function 2927 explicit 2928 hypergeometric_distribution(result_type __N = 10, result_type __K = 5, 2929 result_type __n = 1) 2930 : _M_param{__N, __K, __n} 2931 { } 2932 2933 explicit 2934 hypergeometric_distribution(const param_type& __p) 2935 : _M_param{__p} 2936 { } 2937 2938 /** 2939 * @brief Resets the distribution state. 2940 */ 2941 void 2942 reset() 2943 { } 2944 2945 /** 2946 * @brief Returns the distribution parameter @p N, 2947 * the total number of items. 2948 */ 2949 result_type 2950 total_size() const 2951 { return this->_M_param.total_size(); } 2952 2953 /** 2954 * @brief Returns the distribution parameter @p K, 2955 * the total number of successful items. 2956 */ 2957 result_type 2958 successful_size() const 2959 { return this->_M_param.successful_size(); } 2960 2961 /** 2962 * @brief Returns the total number of unsuccessful items @f$ N - K @f$. 2963 */ 2964 result_type 2965 unsuccessful_size() const 2966 { return this->_M_param.unsuccessful_size(); } 2967 2968 /** 2969 * @brief Returns the distribution parameter @p n, 2970 * the total number of draws. 2971 */ 2972 result_type 2973 total_draws() const 2974 { return this->_M_param.total_draws(); } 2975 2976 /** 2977 * @brief Returns the parameter set of the distribution. 2978 */ 2979 param_type 2980 param() const 2981 { return this->_M_param; } 2982 2983 /** 2984 * @brief Sets the parameter set of the distribution. 2985 * @param __param The new parameter set of the distribution. 2986 */ 2987 void 2988 param(const param_type& __param) 2989 { this->_M_param = __param; } 2990 2991 /** 2992 * @brief Returns the greatest lower bound value of the distribution. 2993 */ 2994 result_type 2995 min() const 2996 { 2997 using _IntType = typename std::make_signed<result_type>::type; 2998 return static_cast<result_type>(std::max(static_cast<_IntType>(0), 2999 static_cast<_IntType>(this->total_draws() 3000 - this->unsuccessful_size()))); 3001 } 3002 3003 /** 3004 * @brief Returns the least upper bound value of the distribution. 3005 */ 3006 result_type 3007 max() const 3008 { return std::min(this->successful_size(), this->total_draws()); } 3009 3010 /** 3011 * @brief Generating functions. 3012 */ 3013 template<typename _UniformRandomNumberGenerator> 3014 result_type 3015 operator()(_UniformRandomNumberGenerator& __urng) 3016 { return this->operator()(__urng, this->_M_param); } 3017 3018 template<typename _UniformRandomNumberGenerator> 3019 result_type 3020 operator()(_UniformRandomNumberGenerator& __urng, 3021 const param_type& __p); 3022 3023 template<typename _ForwardIterator, 3024 typename _UniformRandomNumberGenerator> 3025 void 3026 __generate(_ForwardIterator __f, _ForwardIterator __t, 3027 _UniformRandomNumberGenerator& __urng) 3028 { this->__generate(__f, __t, __urng, this->_M_param); } 3029 3030 template<typename _ForwardIterator, 3031 typename _UniformRandomNumberGenerator> 3032 void 3033 __generate(_ForwardIterator __f, _ForwardIterator __t, 3034 _UniformRandomNumberGenerator& __urng, 3035 const param_type& __p) 3036 { this->__generate_impl(__f, __t, __urng, __p); } 3037 3038 template<typename _UniformRandomNumberGenerator> 3039 void 3040 __generate(result_type* __f, result_type* __t, 3041 _UniformRandomNumberGenerator& __urng, 3042 const param_type& __p) 3043 { this->__generate_impl(__f, __t, __urng, __p); } 3044 3045 /** 3046 * @brief Return true if two hypergeometric distributions have the same 3047 * parameters and the sequences that would be generated 3048 * are equal. 3049 */ 3050 friend bool 3051 operator==(const hypergeometric_distribution& __d1, 3052 const hypergeometric_distribution& __d2) 3053 { return __d1._M_param == __d2._M_param; } 3054 3055 /** 3056 * @brief Inserts a %hypergeometric_distribution random number 3057 * distribution @p __x into the output stream @p __os. 3058 * 3059 * @param __os An output stream. 3060 * @param __x A %hypergeometric_distribution random number 3061 * distribution. 3062 * 3063 * @returns The output stream with the state of @p __x inserted or in 3064 * an error state. 3065 */ 3066 template<typename _UIntType1, typename _CharT, typename _Traits> 3067 friend std::basic_ostream<_CharT, _Traits>& 3068 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 3069 const __gnu_cxx::hypergeometric_distribution<_UIntType1>& 3070 __x); 3071 3072 /** 3073 * @brief Extracts a %hypergeometric_distribution random number 3074 * distribution @p __x from the input stream @p __is. 3075 * 3076 * @param __is An input stream. 3077 * @param __x A %hypergeometric_distribution random number generator 3078 * distribution. 3079 * 3080 * @returns The input stream with @p __x extracted or in an error 3081 * state. 3082 */ 3083 template<typename _UIntType1, typename _CharT, typename _Traits> 3084 friend std::basic_istream<_CharT, _Traits>& 3085 operator>>(std::basic_istream<_CharT, _Traits>& __is, 3086 __gnu_cxx::hypergeometric_distribution<_UIntType1>& __x); 3087 3088 private: 3089 3090 template<typename _ForwardIterator, 3091 typename _UniformRandomNumberGenerator> 3092 void 3093 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 3094 _UniformRandomNumberGenerator& __urng, 3095 const param_type& __p); 3096 3097 param_type _M_param; 3098 }; 3099 3100 /** 3101 * @brief Return true if two hypergeometric distributions are different. 3102 */ 3103 template<typename _UIntType> 3104 inline bool 3105 operator!=(const __gnu_cxx::hypergeometric_distribution<_UIntType>& __d1, 3106 const __gnu_cxx::hypergeometric_distribution<_UIntType>& __d2) 3107 { return !(__d1 == __d2); } 3108 3109_GLIBCXX_END_NAMESPACE_VERSION 3110} // namespace __gnu_cxx 3111 3112#include "ext/opt_random.h" 3113#include "random.tcc" 3114 3115#endif // _GLIBCXX_USE_C99_STDINT_TR1 3116 3117#endif // C++11 3118 3119#endif // _EXT_RANDOM 3120