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