Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

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