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