libstdc++-v3/include/ext/random

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