Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. // random number generation (out of line) -*- C++ -*-
  2.  
  3. // Copyright (C) 2009-2015 Free Software Foundation, Inc.
  4. //
  5. // This file is part of the GNU ISO C++ Library.  This library is free
  6. // software; you can redistribute it and/or modify it under the
  7. // terms of the GNU General Public License as published by the
  8. // Free Software Foundation; either version 3, or (at your option)
  9. // any later version.
  10.  
  11. // This library is distributed in the hope that it will be useful,
  12. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14. // GNU General Public License for more details.
  15.  
  16. // Under Section 7 of GPL version 3, you are granted additional
  17. // permissions described in the GCC Runtime Library Exception, version
  18. // 3.1, as published by the Free Software Foundation.
  19.  
  20. // You should have received a copy of the GNU General Public License and
  21. // a copy of the GCC Runtime Library Exception along with this program;
  22. // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  23. // <http://www.gnu.org/licenses/>.
  24.  
  25.  
  26. /** @file tr1/random.tcc
  27.  *  This is an internal header file, included by other library headers.
  28.  *  Do not attempt to use it directly. @headername{tr1/random}
  29.  */
  30.  
  31. #ifndef _GLIBCXX_TR1_RANDOM_TCC
  32. #define _GLIBCXX_TR1_RANDOM_TCC 1
  33.  
  34. namespace std _GLIBCXX_VISIBILITY(default)
  35. {
  36. namespace tr1
  37. {
  38.   /*
  39.    * (Further) implementation-space details.
  40.    */
  41.   namespace __detail
  42.   {
  43.   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  44.  
  45.     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
  46.     // integer overflow.
  47.     //
  48.     // Because a and c are compile-time integral constants the compiler kindly
  49.     // elides any unreachable paths.
  50.     //
  51.     // Preconditions:  a > 0, m > 0.
  52.     //
  53.     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
  54.       struct _Mod
  55.       {
  56.         static _Tp
  57.         __calc(_Tp __x)
  58.         {
  59.           if (__a == 1)
  60.             __x %= __m;
  61.           else
  62.             {
  63.               static const _Tp __q = __m / __a;
  64.               static const _Tp __r = __m % __a;
  65.              
  66.               _Tp __t1 = __a * (__x % __q);
  67.               _Tp __t2 = __r * (__x / __q);
  68.               if (__t1 >= __t2)
  69.                 __x = __t1 - __t2;
  70.               else
  71.                 __x = __m - __t2 + __t1;
  72.             }
  73.  
  74.           if (__c != 0)
  75.             {
  76.               const _Tp __d = __m - __x;
  77.               if (__d > __c)
  78.                 __x += __c;
  79.               else
  80.                 __x = __c - __d;
  81.             }
  82.           return __x;
  83.         }
  84.       };
  85.  
  86.     // Special case for m == 0 -- use unsigned integer overflow as modulo
  87.     // operator.
  88.     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
  89.       struct _Mod<_Tp, __a, __c, __m, true>
  90.       {
  91.         static _Tp
  92.         __calc(_Tp __x)
  93.         { return __a * __x + __c; }
  94.       };
  95.   _GLIBCXX_END_NAMESPACE_VERSION
  96.   } // namespace __detail
  97.  
  98. _GLIBCXX_BEGIN_NAMESPACE_VERSION
  99.  
  100.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  101.     const _UIntType
  102.     linear_congruential<_UIntType, __a, __c, __m>::multiplier;
  103.  
  104.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  105.     const _UIntType
  106.     linear_congruential<_UIntType, __a, __c, __m>::increment;
  107.  
  108.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  109.     const _UIntType
  110.     linear_congruential<_UIntType, __a, __c, __m>::modulus;
  111.  
  112.   /**
  113.    * Seeds the LCR with integral value @p __x0, adjusted so that the
  114.    * ring identity is never a member of the convergence set.
  115.    */
  116.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  117.     void
  118.     linear_congruential<_UIntType, __a, __c, __m>::
  119.     seed(unsigned long __x0)
  120.     {
  121.       if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
  122.           && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
  123.         _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
  124.       else
  125.         _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
  126.     }
  127.  
  128.   /**
  129.    * Seeds the LCR engine with a value generated by @p __g.
  130.    */
  131.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  132.     template<class _Gen>
  133.       void
  134.       linear_congruential<_UIntType, __a, __c, __m>::
  135.       seed(_Gen& __g, false_type)
  136.       {
  137.         _UIntType __x0 = __g();
  138.         if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
  139.             && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
  140.           _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
  141.         else
  142.           _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
  143.       }
  144.  
  145.   /**
  146.    * Gets the next generated value in sequence.
  147.    */
  148.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  149.     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
  150.     linear_congruential<_UIntType, __a, __c, __m>::
  151.     operator()()
  152.     {
  153.       _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
  154.       return _M_x;
  155.     }
  156.  
  157.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
  158.            typename _CharT, typename _Traits>
  159.     std::basic_ostream<_CharT, _Traits>&
  160.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  161.                const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
  162.     {
  163.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  164.       typedef typename __ostream_type::ios_base    __ios_base;
  165.  
  166.       const typename __ios_base::fmtflags __flags = __os.flags();
  167.       const _CharT __fill = __os.fill();
  168.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
  169.       __os.fill(__os.widen(' '));
  170.  
  171.       __os << __lcr._M_x;
  172.  
  173.       __os.flags(__flags);
  174.       __os.fill(__fill);
  175.       return __os;
  176.     }
  177.  
  178.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
  179.            typename _CharT, typename _Traits>
  180.     std::basic_istream<_CharT, _Traits>&
  181.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  182.                linear_congruential<_UIntType, __a, __c, __m>& __lcr)
  183.     {
  184.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  185.       typedef typename __istream_type::ios_base    __ios_base;
  186.  
  187.       const typename __ios_base::fmtflags __flags = __is.flags();
  188.       __is.flags(__ios_base::dec);
  189.  
  190.       __is >> __lcr._M_x;
  191.  
  192.       __is.flags(__flags);
  193.       return __is;
  194.     }
  195.  
  196.  
  197.   template<class _UIntType, int __w, int __n, int __m, int __r,
  198.            _UIntType __a, int __u, int __s,
  199.            _UIntType __b, int __t, _UIntType __c, int __l>
  200.     const int
  201.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  202.                      __b, __t, __c, __l>::word_size;
  203.  
  204.   template<class _UIntType, int __w, int __n, int __m, int __r,
  205.            _UIntType __a, int __u, int __s,
  206.            _UIntType __b, int __t, _UIntType __c, int __l>
  207.     const int
  208.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  209.                      __b, __t, __c, __l>::state_size;
  210.    
  211.   template<class _UIntType, int __w, int __n, int __m, int __r,
  212.            _UIntType __a, int __u, int __s,
  213.            _UIntType __b, int __t, _UIntType __c, int __l>
  214.     const int
  215.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  216.                      __b, __t, __c, __l>::shift_size;
  217.  
  218.   template<class _UIntType, int __w, int __n, int __m, int __r,
  219.            _UIntType __a, int __u, int __s,
  220.            _UIntType __b, int __t, _UIntType __c, int __l>
  221.     const int
  222.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  223.                      __b, __t, __c, __l>::mask_bits;
  224.  
  225.   template<class _UIntType, int __w, int __n, int __m, int __r,
  226.            _UIntType __a, int __u, int __s,
  227.            _UIntType __b, int __t, _UIntType __c, int __l>
  228.     const _UIntType
  229.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  230.                      __b, __t, __c, __l>::parameter_a;
  231.  
  232.   template<class _UIntType, int __w, int __n, int __m, int __r,
  233.            _UIntType __a, int __u, int __s,
  234.            _UIntType __b, int __t, _UIntType __c, int __l>
  235.     const int
  236.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  237.                      __b, __t, __c, __l>::output_u;
  238.  
  239.   template<class _UIntType, int __w, int __n, int __m, int __r,
  240.            _UIntType __a, int __u, int __s,
  241.            _UIntType __b, int __t, _UIntType __c, int __l>
  242.     const int
  243.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  244.                      __b, __t, __c, __l>::output_s;
  245.  
  246.   template<class _UIntType, int __w, int __n, int __m, int __r,
  247.            _UIntType __a, int __u, int __s,
  248.            _UIntType __b, int __t, _UIntType __c, int __l>
  249.     const _UIntType
  250.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  251.                      __b, __t, __c, __l>::output_b;
  252.  
  253.   template<class _UIntType, int __w, int __n, int __m, int __r,
  254.            _UIntType __a, int __u, int __s,
  255.            _UIntType __b, int __t, _UIntType __c, int __l>
  256.     const int
  257.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  258.                      __b, __t, __c, __l>::output_t;
  259.  
  260.   template<class _UIntType, int __w, int __n, int __m, int __r,
  261.            _UIntType __a, int __u, int __s,
  262.            _UIntType __b, int __t, _UIntType __c, int __l>
  263.     const _UIntType
  264.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  265.                      __b, __t, __c, __l>::output_c;
  266.  
  267.   template<class _UIntType, int __w, int __n, int __m, int __r,
  268.            _UIntType __a, int __u, int __s,
  269.            _UIntType __b, int __t, _UIntType __c, int __l>
  270.     const int
  271.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  272.                      __b, __t, __c, __l>::output_l;
  273.  
  274.   template<class _UIntType, int __w, int __n, int __m, int __r,
  275.            _UIntType __a, int __u, int __s,
  276.            _UIntType __b, int __t, _UIntType __c, int __l>
  277.     void
  278.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  279.                      __b, __t, __c, __l>::
  280.     seed(unsigned long __value)
  281.     {
  282.       _M_x[0] = __detail::__mod<_UIntType, 1, 0,
  283.         __detail::_Shift<_UIntType, __w>::__value>(__value);
  284.  
  285.       for (int __i = 1; __i < state_size; ++__i)
  286.         {
  287.           _UIntType __x = _M_x[__i - 1];
  288.           __x ^= __x >> (__w - 2);
  289.           __x *= 1812433253ul;
  290.           __x += __i;
  291.           _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
  292.             __detail::_Shift<_UIntType, __w>::__value>(__x);     
  293.         }
  294.       _M_p = state_size;
  295.     }
  296.  
  297.   template<class _UIntType, int __w, int __n, int __m, int __r,
  298.            _UIntType __a, int __u, int __s,
  299.            _UIntType __b, int __t, _UIntType __c, int __l>
  300.     template<class _Gen>
  301.       void
  302.       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  303.                        __b, __t, __c, __l>::
  304.       seed(_Gen& __gen, false_type)
  305.       {
  306.         for (int __i = 0; __i < state_size; ++__i)
  307.           _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
  308.             __detail::_Shift<_UIntType, __w>::__value>(__gen());
  309.         _M_p = state_size;
  310.       }
  311.  
  312.   template<class _UIntType, int __w, int __n, int __m, int __r,
  313.            _UIntType __a, int __u, int __s,
  314.            _UIntType __b, int __t, _UIntType __c, int __l>
  315.     typename
  316.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  317.                      __b, __t, __c, __l>::result_type
  318.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
  319.                      __b, __t, __c, __l>::
  320.     operator()()
  321.     {
  322.       // Reload the vector - cost is O(n) amortized over n calls.
  323.       if (_M_p >= state_size)
  324.         {
  325.           const _UIntType __upper_mask = (~_UIntType()) << __r;
  326.           const _UIntType __lower_mask = ~__upper_mask;
  327.  
  328.           for (int __k = 0; __k < (__n - __m); ++__k)
  329.             {
  330.               _UIntType __y = ((_M_x[__k] & __upper_mask)
  331.                                | (_M_x[__k + 1] & __lower_mask));
  332.               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
  333.                            ^ ((__y & 0x01) ? __a : 0));
  334.             }
  335.  
  336.           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
  337.             {
  338.               _UIntType __y = ((_M_x[__k] & __upper_mask)
  339.                                | (_M_x[__k + 1] & __lower_mask));
  340.               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
  341.                            ^ ((__y & 0x01) ? __a : 0));
  342.             }
  343.  
  344.           _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
  345.                            | (_M_x[0] & __lower_mask));
  346.           _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
  347.                            ^ ((__y & 0x01) ? __a : 0));
  348.           _M_p = 0;
  349.         }
  350.  
  351.       // Calculate o(x(i)).
  352.       result_type __z = _M_x[_M_p++];
  353.       __z ^= (__z >> __u);
  354.       __z ^= (__z << __s) & __b;
  355.       __z ^= (__z << __t) & __c;
  356.       __z ^= (__z >> __l);
  357.  
  358.       return __z;
  359.     }
  360.  
  361.   template<class _UIntType, int __w, int __n, int __m, int __r,
  362.            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
  363.            _UIntType __c, int __l,
  364.            typename _CharT, typename _Traits>
  365.     std::basic_ostream<_CharT, _Traits>&
  366.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  367.                const mersenne_twister<_UIntType, __w, __n, __m,
  368.                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
  369.     {
  370.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  371.       typedef typename __ostream_type::ios_base    __ios_base;
  372.  
  373.       const typename __ios_base::fmtflags __flags = __os.flags();
  374.       const _CharT __fill = __os.fill();
  375.       const _CharT __space = __os.widen(' ');
  376.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
  377.       __os.fill(__space);
  378.  
  379.       for (int __i = 0; __i < __n - 1; ++__i)
  380.         __os << __x._M_x[__i] << __space;
  381.       __os << __x._M_x[__n - 1];
  382.  
  383.       __os.flags(__flags);
  384.       __os.fill(__fill);
  385.       return __os;
  386.     }
  387.  
  388.   template<class _UIntType, int __w, int __n, int __m, int __r,
  389.            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
  390.            _UIntType __c, int __l,
  391.            typename _CharT, typename _Traits>
  392.     std::basic_istream<_CharT, _Traits>&
  393.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  394.                mersenne_twister<_UIntType, __w, __n, __m,
  395.                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
  396.     {
  397.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  398.       typedef typename __istream_type::ios_base    __ios_base;
  399.  
  400.       const typename __ios_base::fmtflags __flags = __is.flags();
  401.       __is.flags(__ios_base::dec | __ios_base::skipws);
  402.  
  403.       for (int __i = 0; __i < __n; ++__i)
  404.         __is >> __x._M_x[__i];
  405.  
  406.       __is.flags(__flags);
  407.       return __is;
  408.     }
  409.  
  410.  
  411.   template<typename _IntType, _IntType __m, int __s, int __r>
  412.     const _IntType
  413.     subtract_with_carry<_IntType, __m, __s, __r>::modulus;
  414.  
  415.   template<typename _IntType, _IntType __m, int __s, int __r>
  416.     const int
  417.     subtract_with_carry<_IntType, __m, __s, __r>::long_lag;
  418.  
  419.   template<typename _IntType, _IntType __m, int __s, int __r>
  420.     const int
  421.     subtract_with_carry<_IntType, __m, __s, __r>::short_lag;
  422.  
  423.   template<typename _IntType, _IntType __m, int __s, int __r>
  424.     void
  425.     subtract_with_carry<_IntType, __m, __s, __r>::
  426.     seed(unsigned long __value)
  427.     {
  428.       if (__value == 0)
  429.         __value = 19780503;
  430.  
  431.       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
  432.         __lcg(__value);
  433.  
  434.       for (int __i = 0; __i < long_lag; ++__i)
  435.         _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
  436.  
  437.       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
  438.       _M_p = 0;
  439.     }
  440.  
  441.   template<typename _IntType, _IntType __m, int __s, int __r>
  442.     template<class _Gen>
  443.       void
  444.       subtract_with_carry<_IntType, __m, __s, __r>::
  445.       seed(_Gen& __gen, false_type)
  446.       {
  447.         const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
  448.  
  449.         for (int __i = 0; __i < long_lag; ++__i)
  450.           {
  451.             _UIntType __tmp = 0;
  452.             _UIntType __factor = 1;
  453.             for (int __j = 0; __j < __n; ++__j)
  454.               {
  455.                 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
  456.                          (__gen()) * __factor;
  457.                 __factor *= __detail::_Shift<_UIntType, 32>::__value;
  458.               }
  459.             _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
  460.           }
  461.         _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
  462.         _M_p = 0;
  463.       }
  464.  
  465.   template<typename _IntType, _IntType __m, int __s, int __r>
  466.     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
  467.     subtract_with_carry<_IntType, __m, __s, __r>::
  468.     operator()()
  469.     {
  470.       // Derive short lag index from current index.
  471.       int __ps = _M_p - short_lag;
  472.       if (__ps < 0)
  473.         __ps += long_lag;
  474.  
  475.       // Calculate new x(i) without overflow or division.
  476.       // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
  477.       // cannot overflow.
  478.       _UIntType __xi;
  479.       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
  480.         {
  481.           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
  482.           _M_carry = 0;
  483.         }
  484.       else
  485.         {
  486.           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
  487.           _M_carry = 1;
  488.         }
  489.       _M_x[_M_p] = __xi;
  490.  
  491.       // Adjust current index to loop around in ring buffer.
  492.       if (++_M_p >= long_lag)
  493.         _M_p = 0;
  494.  
  495.       return __xi;
  496.     }
  497.  
  498.   template<typename _IntType, _IntType __m, int __s, int __r,
  499.            typename _CharT, typename _Traits>
  500.     std::basic_ostream<_CharT, _Traits>&
  501.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  502.                const subtract_with_carry<_IntType, __m, __s, __r>& __x)
  503.     {
  504.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  505.       typedef typename __ostream_type::ios_base    __ios_base;
  506.  
  507.       const typename __ios_base::fmtflags __flags = __os.flags();
  508.       const _CharT __fill = __os.fill();
  509.       const _CharT __space = __os.widen(' ');
  510.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
  511.       __os.fill(__space);
  512.  
  513.       for (int __i = 0; __i < __r; ++__i)
  514.         __os << __x._M_x[__i] << __space;
  515.       __os << __x._M_carry;
  516.  
  517.       __os.flags(__flags);
  518.       __os.fill(__fill);
  519.       return __os;
  520.     }
  521.  
  522.   template<typename _IntType, _IntType __m, int __s, int __r,
  523.            typename _CharT, typename _Traits>
  524.     std::basic_istream<_CharT, _Traits>&
  525.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  526.                subtract_with_carry<_IntType, __m, __s, __r>& __x)
  527.     {
  528.       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
  529.       typedef typename __istream_type::ios_base    __ios_base;
  530.  
  531.       const typename __ios_base::fmtflags __flags = __is.flags();
  532.       __is.flags(__ios_base::dec | __ios_base::skipws);
  533.  
  534.       for (int __i = 0; __i < __r; ++__i)
  535.         __is >> __x._M_x[__i];
  536.       __is >> __x._M_carry;
  537.  
  538.       __is.flags(__flags);
  539.       return __is;
  540.     }
  541.  
  542.  
  543.   template<typename _RealType, int __w, int __s, int __r>
  544.     const int
  545.     subtract_with_carry_01<_RealType, __w, __s, __r>::word_size;
  546.  
  547.   template<typename _RealType, int __w, int __s, int __r>
  548.     const int
  549.     subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag;
  550.  
  551.   template<typename _RealType, int __w, int __s, int __r>
  552.     const int
  553.     subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag;
  554.  
  555.   template<typename _RealType, int __w, int __s, int __r>
  556.     void
  557.     subtract_with_carry_01<_RealType, __w, __s, __r>::
  558.     _M_initialize_npows()
  559.     {
  560.       for (int __j = 0; __j < __n; ++__j)
  561. #if _GLIBCXX_USE_C99_MATH_TR1
  562.         _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
  563. #else
  564.         _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
  565. #endif
  566.     }
  567.  
  568.   template<typename _RealType, int __w, int __s, int __r>
  569.     void
  570.     subtract_with_carry_01<_RealType, __w, __s, __r>::
  571.     seed(unsigned long __value)
  572.     {
  573.       if (__value == 0)
  574.         __value = 19780503;
  575.  
  576.       // _GLIBCXX_RESOLVE_LIB_DEFECTS
  577.       // 512. Seeding subtract_with_carry_01 from a single unsigned long.
  578.       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
  579.         __lcg(__value);
  580.  
  581.       this->seed(__lcg);
  582.     }
  583.  
  584.   template<typename _RealType, int __w, int __s, int __r>
  585.     template<class _Gen>
  586.       void
  587.       subtract_with_carry_01<_RealType, __w, __s, __r>::
  588.       seed(_Gen& __gen, false_type)
  589.       {
  590.         for (int __i = 0; __i < long_lag; ++__i)
  591.           {
  592.             for (int __j = 0; __j < __n - 1; ++__j)
  593.               _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
  594.             _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
  595.               __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
  596.           }
  597.  
  598.         _M_carry = 1;
  599.         for (int __j = 0; __j < __n; ++__j)
  600.           if (_M_x[long_lag - 1][__j] != 0)
  601.             {
  602.               _M_carry = 0;
  603.               break;
  604.             }
  605.  
  606.         _M_p = 0;
  607.       }
  608.  
  609.   template<typename _RealType, int __w, int __s, int __r>
  610.     typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
  611.     subtract_with_carry_01<_RealType, __w, __s, __r>::
  612.     operator()()
  613.     {
  614.       // Derive short lag index from current index.
  615.       int __ps = _M_p - short_lag;
  616.       if (__ps < 0)
  617.         __ps += long_lag;
  618.  
  619.       _UInt32Type __new_carry;
  620.       for (int __j = 0; __j < __n - 1; ++__j)
  621.         {
  622.           if (_M_x[__ps][__j] > _M_x[_M_p][__j]
  623.               || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
  624.             __new_carry = 0;
  625.           else
  626.             __new_carry = 1;
  627.  
  628.           _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
  629.           _M_carry = __new_carry;
  630.         }
  631.  
  632.       if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
  633.           || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
  634.         __new_carry = 0;
  635.       else
  636.         __new_carry = 1;
  637.      
  638.       _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
  639.         __detail::_Shift<_UInt32Type, __w % 32>::__value>
  640.         (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
  641.       _M_carry = __new_carry;
  642.  
  643.       result_type __ret = 0.0;
  644.       for (int __j = 0; __j < __n; ++__j)
  645.         __ret += _M_x[_M_p][__j] * _M_npows[__j];
  646.  
  647.       // Adjust current index to loop around in ring buffer.
  648.       if (++_M_p >= long_lag)
  649.         _M_p = 0;
  650.  
  651.       return __ret;
  652.     }
  653.  
  654.   template<typename _RealType, int __w, int __s, int __r,
  655.            typename _CharT, typename _Traits>
  656.     std::basic_ostream<_CharT, _Traits>&
  657.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  658.                const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
  659.     {
  660.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  661.       typedef typename __ostream_type::ios_base    __ios_base;
  662.  
  663.       const typename __ios_base::fmtflags __flags = __os.flags();
  664.       const _CharT __fill = __os.fill();
  665.       const _CharT __space = __os.widen(' ');
  666.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
  667.       __os.fill(__space);
  668.  
  669.       for (int __i = 0; __i < __r; ++__i)
  670.         for (int __j = 0; __j < __x.__n; ++__j)
  671.           __os << __x._M_x[__i][__j] << __space;
  672.       __os << __x._M_carry;
  673.  
  674.       __os.flags(__flags);
  675.       __os.fill(__fill);
  676.       return __os;
  677.     }
  678.  
  679.   template<typename _RealType, int __w, int __s, int __r,
  680.            typename _CharT, typename _Traits>
  681.     std::basic_istream<_CharT, _Traits>&
  682.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  683.                subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
  684.     {
  685.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  686.       typedef typename __istream_type::ios_base    __ios_base;
  687.  
  688.       const typename __ios_base::fmtflags __flags = __is.flags();
  689.       __is.flags(__ios_base::dec | __ios_base::skipws);
  690.  
  691.       for (int __i = 0; __i < __r; ++__i)
  692.         for (int __j = 0; __j < __x.__n; ++__j)
  693.           __is >> __x._M_x[__i][__j];
  694.       __is >> __x._M_carry;
  695.  
  696.       __is.flags(__flags);
  697.       return __is;
  698.     }
  699.  
  700.   template<class _UniformRandomNumberGenerator, int __p, int __r>
  701.     const int
  702.     discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size;
  703.  
  704.   template<class _UniformRandomNumberGenerator, int __p, int __r>
  705.     const int
  706.     discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block;
  707.  
  708.   template<class _UniformRandomNumberGenerator, int __p, int __r>
  709.     typename discard_block<_UniformRandomNumberGenerator,
  710.                            __p, __r>::result_type
  711.     discard_block<_UniformRandomNumberGenerator, __p, __r>::
  712.     operator()()
  713.     {
  714.       if (_M_n >= used_block)
  715.         {
  716.           while (_M_n < block_size)
  717.             {
  718.               _M_b();
  719.               ++_M_n;
  720.             }
  721.           _M_n = 0;
  722.         }
  723.       ++_M_n;
  724.       return _M_b();
  725.     }
  726.  
  727.   template<class _UniformRandomNumberGenerator, int __p, int __r,
  728.            typename _CharT, typename _Traits>
  729.     std::basic_ostream<_CharT, _Traits>&
  730.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  731.                const discard_block<_UniformRandomNumberGenerator,
  732.                __p, __r>& __x)
  733.     {
  734.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  735.       typedef typename __ostream_type::ios_base    __ios_base;
  736.  
  737.       const typename __ios_base::fmtflags __flags = __os.flags();
  738.       const _CharT __fill = __os.fill();
  739.       const _CharT __space = __os.widen(' ');
  740.       __os.flags(__ios_base::dec | __ios_base::fixed
  741.                  | __ios_base::left);
  742.       __os.fill(__space);
  743.  
  744.       __os << __x._M_b << __space << __x._M_n;
  745.  
  746.       __os.flags(__flags);
  747.       __os.fill(__fill);
  748.       return __os;
  749.     }
  750.  
  751.   template<class _UniformRandomNumberGenerator, int __p, int __r,
  752.            typename _CharT, typename _Traits>
  753.     std::basic_istream<_CharT, _Traits>&
  754.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  755.                discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
  756.     {
  757.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  758.       typedef typename __istream_type::ios_base    __ios_base;
  759.  
  760.       const typename __ios_base::fmtflags __flags = __is.flags();
  761.       __is.flags(__ios_base::dec | __ios_base::skipws);
  762.  
  763.       __is >> __x._M_b >> __x._M_n;
  764.  
  765.       __is.flags(__flags);
  766.       return __is;
  767.     }
  768.  
  769.  
  770.   template<class _UniformRandomNumberGenerator1, int __s1,
  771.            class _UniformRandomNumberGenerator2, int __s2>
  772.     const int
  773.     xor_combine<_UniformRandomNumberGenerator1, __s1,
  774.                 _UniformRandomNumberGenerator2, __s2>::shift1;
  775.      
  776.   template<class _UniformRandomNumberGenerator1, int __s1,
  777.            class _UniformRandomNumberGenerator2, int __s2>
  778.     const int
  779.     xor_combine<_UniformRandomNumberGenerator1, __s1,
  780.                 _UniformRandomNumberGenerator2, __s2>::shift2;
  781.  
  782.   template<class _UniformRandomNumberGenerator1, int __s1,
  783.            class _UniformRandomNumberGenerator2, int __s2>
  784.     void
  785.     xor_combine<_UniformRandomNumberGenerator1, __s1,
  786.                 _UniformRandomNumberGenerator2, __s2>::
  787.     _M_initialize_max()
  788.     {
  789.       const int __w = std::numeric_limits<result_type>::digits;
  790.  
  791.       const result_type __m1 =
  792.         std::min(result_type(_M_b1.max() - _M_b1.min()),
  793.                  __detail::_Shift<result_type, __w - __s1>::__value - 1);
  794.  
  795.       const result_type __m2 =
  796.         std::min(result_type(_M_b2.max() - _M_b2.min()),
  797.                  __detail::_Shift<result_type, __w - __s2>::__value - 1);
  798.  
  799.       // NB: In TR1 s1 is not required to be >= s2.
  800.       if (__s1 < __s2)
  801.         _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
  802.       else
  803.         _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
  804.     }
  805.  
  806.   template<class _UniformRandomNumberGenerator1, int __s1,
  807.            class _UniformRandomNumberGenerator2, int __s2>
  808.     typename xor_combine<_UniformRandomNumberGenerator1, __s1,
  809.                          _UniformRandomNumberGenerator2, __s2>::result_type
  810.     xor_combine<_UniformRandomNumberGenerator1, __s1,
  811.                 _UniformRandomNumberGenerator2, __s2>::
  812.     _M_initialize_max_aux(result_type __a, result_type __b, int __d)
  813.     {
  814.       const result_type __two2d = result_type(1) << __d;
  815.       const result_type __c = __a * __two2d;
  816.  
  817.       if (__a == 0 || __b < __two2d)
  818.         return __c + __b;
  819.  
  820.       const result_type __t = std::max(__c, __b);
  821.       const result_type __u = std::min(__c, __b);
  822.  
  823.       result_type __ub = __u;
  824.       result_type __p;
  825.       for (__p = 0; __ub != 1; __ub >>= 1)
  826.         ++__p;
  827.  
  828.       const result_type __two2p = result_type(1) << __p;
  829.       const result_type __k = __t / __two2p;
  830.  
  831.       if (__k & 1)
  832.         return (__k + 1) * __two2p - 1;
  833.  
  834.       if (__c >= __b)
  835.         return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
  836.                                                            / __two2d,
  837.                                                            __u % __two2p, __d);
  838.       else
  839.         return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
  840.                                                            / __two2d,
  841.                                                            __t % __two2p, __d);
  842.     }
  843.  
  844.   template<class _UniformRandomNumberGenerator1, int __s1,
  845.            class _UniformRandomNumberGenerator2, int __s2,
  846.            typename _CharT, typename _Traits>
  847.     std::basic_ostream<_CharT, _Traits>&
  848.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  849.                const xor_combine<_UniformRandomNumberGenerator1, __s1,
  850.                _UniformRandomNumberGenerator2, __s2>& __x)
  851.     {
  852.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  853.       typedef typename __ostream_type::ios_base    __ios_base;
  854.  
  855.       const typename __ios_base::fmtflags __flags = __os.flags();
  856.       const _CharT __fill = __os.fill();
  857.       const _CharT __space = __os.widen(' ');
  858.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
  859.       __os.fill(__space);
  860.  
  861.       __os << __x.base1() << __space << __x.base2();
  862.  
  863.       __os.flags(__flags);
  864.       __os.fill(__fill);
  865.       return __os;
  866.     }
  867.  
  868.   template<class _UniformRandomNumberGenerator1, int __s1,
  869.            class _UniformRandomNumberGenerator2, int __s2,
  870.            typename _CharT, typename _Traits>
  871.     std::basic_istream<_CharT, _Traits>&
  872.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  873.                xor_combine<_UniformRandomNumberGenerator1, __s1,
  874.                _UniformRandomNumberGenerator2, __s2>& __x)
  875.     {
  876.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  877.       typedef typename __istream_type::ios_base    __ios_base;
  878.  
  879.       const typename __ios_base::fmtflags __flags = __is.flags();
  880.       __is.flags(__ios_base::skipws);
  881.  
  882.       __is >> __x._M_b1 >> __x._M_b2;
  883.  
  884.       __is.flags(__flags);
  885.       return __is;
  886.     }
  887.  
  888.  
  889.   template<typename _IntType>
  890.     template<typename _UniformRandomNumberGenerator>
  891.       typename uniform_int<_IntType>::result_type
  892.       uniform_int<_IntType>::
  893.       _M_call(_UniformRandomNumberGenerator& __urng,
  894.               result_type __min, result_type __max, true_type)
  895.       {
  896.         // XXX Must be fixed to work well for *arbitrary* __urng.max(),
  897.         // __urng.min(), __max, __min.  Currently works fine only in the
  898.         // most common case __urng.max() - __urng.min() >= __max - __min,
  899.         // with __urng.max() > __urng.min() >= 0.
  900.         typedef typename __gnu_cxx::__add_unsigned<typename
  901.           _UniformRandomNumberGenerator::result_type>::__type __urntype;
  902.         typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
  903.                                                               __utype;
  904.         typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
  905.                                                         > sizeof(__utype)),
  906.           __urntype, __utype>::__type                         __uctype;
  907.  
  908.         result_type __ret;
  909.  
  910.         const __urntype __urnmin = __urng.min();
  911.         const __urntype __urnmax = __urng.max();
  912.         const __urntype __urnrange = __urnmax - __urnmin;
  913.         const __uctype __urange = __max - __min;
  914.         const __uctype __udenom = (__urnrange <= __urange
  915.                                    ? 1 : __urnrange / (__urange + 1));
  916.         do
  917.           __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
  918.         while (__ret > __max - __min);
  919.  
  920.         return __ret + __min;
  921.       }
  922.  
  923.   template<typename _IntType, typename _CharT, typename _Traits>
  924.     std::basic_ostream<_CharT, _Traits>&
  925.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  926.                const uniform_int<_IntType>& __x)
  927.     {
  928.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  929.       typedef typename __ostream_type::ios_base    __ios_base;
  930.  
  931.       const typename __ios_base::fmtflags __flags = __os.flags();
  932.       const _CharT __fill = __os.fill();
  933.       const _CharT __space = __os.widen(' ');
  934.       __os.flags(__ios_base::scientific | __ios_base::left);
  935.       __os.fill(__space);
  936.  
  937.       __os << __x.min() << __space << __x.max();
  938.  
  939.       __os.flags(__flags);
  940.       __os.fill(__fill);
  941.       return __os;
  942.     }
  943.  
  944.   template<typename _IntType, typename _CharT, typename _Traits>
  945.     std::basic_istream<_CharT, _Traits>&
  946.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  947.                uniform_int<_IntType>& __x)
  948.     {
  949.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  950.       typedef typename __istream_type::ios_base    __ios_base;
  951.  
  952.       const typename __ios_base::fmtflags __flags = __is.flags();
  953.       __is.flags(__ios_base::dec | __ios_base::skipws);
  954.  
  955.       __is >> __x._M_min >> __x._M_max;
  956.  
  957.       __is.flags(__flags);
  958.       return __is;
  959.     }
  960.  
  961.  
  962.   template<typename _CharT, typename _Traits>
  963.     std::basic_ostream<_CharT, _Traits>&
  964.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  965.                const bernoulli_distribution& __x)
  966.     {
  967.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  968.       typedef typename __ostream_type::ios_base    __ios_base;
  969.  
  970.       const typename __ios_base::fmtflags __flags = __os.flags();
  971.       const _CharT __fill = __os.fill();
  972.       const std::streamsize __precision = __os.precision();
  973.       __os.flags(__ios_base::scientific | __ios_base::left);
  974.       __os.fill(__os.widen(' '));
  975.       __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
  976.  
  977.       __os << __x.p();
  978.  
  979.       __os.flags(__flags);
  980.       __os.fill(__fill);
  981.       __os.precision(__precision);
  982.       return __os;
  983.     }
  984.  
  985.  
  986.   template<typename _IntType, typename _RealType>
  987.     template<class _UniformRandomNumberGenerator>
  988.       typename geometric_distribution<_IntType, _RealType>::result_type
  989.       geometric_distribution<_IntType, _RealType>::
  990.       operator()(_UniformRandomNumberGenerator& __urng)
  991.       {
  992.         // About the epsilon thing see this thread:
  993.         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
  994.         const _RealType __naf =
  995.           (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
  996.         // The largest _RealType convertible to _IntType.
  997.         const _RealType __thr =
  998.           std::numeric_limits<_IntType>::max() + __naf;
  999.  
  1000.         _RealType __cand;
  1001.         do
  1002.           __cand = std::ceil(std::log(__urng()) / _M_log_p);
  1003.         while (__cand >= __thr);
  1004.  
  1005.         return result_type(__cand + __naf);
  1006.       }
  1007.  
  1008.   template<typename _IntType, typename _RealType,
  1009.            typename _CharT, typename _Traits>
  1010.     std::basic_ostream<_CharT, _Traits>&
  1011.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1012.                const geometric_distribution<_IntType, _RealType>& __x)
  1013.     {
  1014.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1015.       typedef typename __ostream_type::ios_base    __ios_base;
  1016.  
  1017.       const typename __ios_base::fmtflags __flags = __os.flags();
  1018.       const _CharT __fill = __os.fill();
  1019.       const std::streamsize __precision = __os.precision();
  1020.       __os.flags(__ios_base::scientific | __ios_base::left);
  1021.       __os.fill(__os.widen(' '));
  1022.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1023.  
  1024.       __os << __x.p();
  1025.  
  1026.       __os.flags(__flags);
  1027.       __os.fill(__fill);
  1028.       __os.precision(__precision);
  1029.       return __os;
  1030.     }
  1031.  
  1032.  
  1033.   template<typename _IntType, typename _RealType>
  1034.     void
  1035.     poisson_distribution<_IntType, _RealType>::
  1036.     _M_initialize()
  1037.     {
  1038. #if _GLIBCXX_USE_C99_MATH_TR1
  1039.       if (_M_mean >= 12)
  1040.         {
  1041.           const _RealType __m = std::floor(_M_mean);
  1042.           _M_lm_thr = std::log(_M_mean);
  1043.           _M_lfm = std::tr1::lgamma(__m + 1);
  1044.           _M_sm = std::sqrt(__m);
  1045.  
  1046.           const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
  1047.           const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
  1048.                                                               / __pi_4));
  1049.           _M_d = std::tr1::round(std::max(_RealType(6),
  1050.                                           std::min(__m, __dx)));
  1051.           const _RealType __cx = 2 * __m + _M_d;
  1052.           _M_scx = std::sqrt(__cx / 2);
  1053.           _M_1cx = 1 / __cx;
  1054.  
  1055.           _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
  1056.           _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
  1057.         }
  1058.       else
  1059. #endif
  1060.         _M_lm_thr = std::exp(-_M_mean);
  1061.       }
  1062.  
  1063.   /**
  1064.    * A rejection algorithm when mean >= 12 and a simple method based
  1065.    * upon the multiplication of uniform random variates otherwise.
  1066.    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
  1067.    * is defined.
  1068.    *
  1069.    * Reference:
  1070.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
  1071.    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
  1072.    */
  1073.   template<typename _IntType, typename _RealType>
  1074.     template<class _UniformRandomNumberGenerator>
  1075.       typename poisson_distribution<_IntType, _RealType>::result_type
  1076.       poisson_distribution<_IntType, _RealType>::
  1077.       operator()(_UniformRandomNumberGenerator& __urng)
  1078.       {
  1079. #if _GLIBCXX_USE_C99_MATH_TR1
  1080.         if (_M_mean >= 12)
  1081.           {
  1082.             _RealType __x;
  1083.  
  1084.             // See comments above...
  1085.             const _RealType __naf =
  1086.               (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
  1087.             const _RealType __thr =
  1088.               std::numeric_limits<_IntType>::max() + __naf;
  1089.  
  1090.             const _RealType __m = std::floor(_M_mean);
  1091.             // sqrt(pi / 2)
  1092.             const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
  1093.             const _RealType __c1 = _M_sm * __spi_2;
  1094.             const _RealType __c2 = _M_c2b + __c1;
  1095.             const _RealType __c3 = __c2 + 1;
  1096.             const _RealType __c4 = __c3 + 1;
  1097.             // e^(1 / 78)
  1098.             const _RealType __e178 = 1.0129030479320018583185514777512983L;
  1099.             const _RealType __c5 = __c4 + __e178;
  1100.             const _RealType __c = _M_cb + __c5;
  1101.             const _RealType __2cx = 2 * (2 * __m + _M_d);
  1102.  
  1103.             bool __reject = true;
  1104.             do
  1105.               {
  1106.                 const _RealType __u = __c * __urng();
  1107.                 const _RealType __e = -std::log(__urng());
  1108.  
  1109.                 _RealType __w = 0.0;
  1110.                
  1111.                 if (__u <= __c1)
  1112.                   {
  1113.                     const _RealType __n = _M_nd(__urng);
  1114.                     const _RealType __y = -std::abs(__n) * _M_sm - 1;
  1115.                     __x = std::floor(__y);
  1116.                     __w = -__n * __n / 2;
  1117.                     if (__x < -__m)
  1118.                       continue;
  1119.                   }
  1120.                 else if (__u <= __c2)
  1121.                   {
  1122.                     const _RealType __n = _M_nd(__urng);
  1123.                     const _RealType __y = 1 + std::abs(__n) * _M_scx;
  1124.                     __x = std::ceil(__y);
  1125.                     __w = __y * (2 - __y) * _M_1cx;
  1126.                     if (__x > _M_d)
  1127.                       continue;
  1128.                   }
  1129.                 else if (__u <= __c3)
  1130.                   // NB: This case not in the book, nor in the Errata,
  1131.                   // but should be ok...
  1132.                   __x = -1;
  1133.                 else if (__u <= __c4)
  1134.                   __x = 0;
  1135.                 else if (__u <= __c5)
  1136.                   __x = 1;
  1137.                 else
  1138.                   {
  1139.                     const _RealType __v = -std::log(__urng());
  1140.                     const _RealType __y = _M_d + __v * __2cx / _M_d;
  1141.                     __x = std::ceil(__y);
  1142.                     __w = -_M_d * _M_1cx * (1 + __y / 2);
  1143.                   }
  1144.  
  1145.                 __reject = (__w - __e - __x * _M_lm_thr
  1146.                             > _M_lfm - std::tr1::lgamma(__x + __m + 1));
  1147.  
  1148.                 __reject |= __x + __m >= __thr;
  1149.  
  1150.               } while (__reject);
  1151.  
  1152.             return result_type(__x + __m + __naf);
  1153.           }
  1154.         else
  1155. #endif
  1156.           {
  1157.             _IntType     __x = 0;
  1158.             _RealType __prod = 1.0;
  1159.  
  1160.             do
  1161.               {
  1162.                 __prod *= __urng();
  1163.                 __x += 1;
  1164.               }
  1165.             while (__prod > _M_lm_thr);
  1166.  
  1167.             return __x - 1;
  1168.           }
  1169.       }
  1170.  
  1171.   template<typename _IntType, typename _RealType,
  1172.            typename _CharT, typename _Traits>
  1173.     std::basic_ostream<_CharT, _Traits>&
  1174.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1175.                const poisson_distribution<_IntType, _RealType>& __x)
  1176.     {
  1177.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1178.       typedef typename __ostream_type::ios_base    __ios_base;
  1179.  
  1180.       const typename __ios_base::fmtflags __flags = __os.flags();
  1181.       const _CharT __fill = __os.fill();
  1182.       const std::streamsize __precision = __os.precision();
  1183.       const _CharT __space = __os.widen(' ');
  1184.       __os.flags(__ios_base::scientific | __ios_base::left);
  1185.       __os.fill(__space);
  1186.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1187.  
  1188.       __os << __x.mean() << __space << __x._M_nd;
  1189.  
  1190.       __os.flags(__flags);
  1191.       __os.fill(__fill);
  1192.       __os.precision(__precision);
  1193.       return __os;
  1194.     }
  1195.  
  1196.   template<typename _IntType, typename _RealType,
  1197.            typename _CharT, typename _Traits>
  1198.     std::basic_istream<_CharT, _Traits>&
  1199.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  1200.                poisson_distribution<_IntType, _RealType>& __x)
  1201.     {
  1202.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  1203.       typedef typename __istream_type::ios_base    __ios_base;
  1204.  
  1205.       const typename __ios_base::fmtflags __flags = __is.flags();
  1206.       __is.flags(__ios_base::skipws);
  1207.  
  1208.       __is >> __x._M_mean >> __x._M_nd;
  1209.       __x._M_initialize();
  1210.  
  1211.       __is.flags(__flags);
  1212.       return __is;
  1213.     }
  1214.  
  1215.  
  1216.   template<typename _IntType, typename _RealType>
  1217.     void
  1218.     binomial_distribution<_IntType, _RealType>::
  1219.     _M_initialize()
  1220.     {
  1221.       const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
  1222.  
  1223.       _M_easy = true;
  1224.  
  1225. #if _GLIBCXX_USE_C99_MATH_TR1
  1226.       if (_M_t * __p12 >= 8)
  1227.         {
  1228.           _M_easy = false;
  1229.           const _RealType __np = std::floor(_M_t * __p12);
  1230.           const _RealType __pa = __np / _M_t;
  1231.           const _RealType __1p = 1 - __pa;
  1232.          
  1233.           const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
  1234.           const _RealType __d1x =
  1235.             std::sqrt(__np * __1p * std::log(32 * __np
  1236.                                              / (81 * __pi_4 * __1p)));
  1237.           _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
  1238.           const _RealType __d2x =
  1239.             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
  1240.                                              / (__pi_4 * __pa)));
  1241.           _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
  1242.          
  1243.           // sqrt(pi / 2)
  1244.           const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
  1245.           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
  1246.           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
  1247.           _M_c = 2 * _M_d1 / __np;
  1248.           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
  1249.           const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
  1250.           const _RealType __s1s = _M_s1 * _M_s1;
  1251.           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
  1252.                              * 2 * __s1s / _M_d1
  1253.                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
  1254.           const _RealType __s2s = _M_s2 * _M_s2;
  1255.           _M_s = (_M_a123 + 2 * __s2s / _M_d2
  1256.                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
  1257.           _M_lf = (std::tr1::lgamma(__np + 1)
  1258.                    + std::tr1::lgamma(_M_t - __np + 1));
  1259.           _M_lp1p = std::log(__pa / __1p);
  1260.  
  1261.           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
  1262.         }
  1263.       else
  1264. #endif
  1265.         _M_q = -std::log(1 - __p12);
  1266.     }
  1267.  
  1268.   template<typename _IntType, typename _RealType>
  1269.     template<class _UniformRandomNumberGenerator>
  1270.       typename binomial_distribution<_IntType, _RealType>::result_type
  1271.       binomial_distribution<_IntType, _RealType>::
  1272.       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
  1273.       {
  1274.         _IntType    __x = 0;
  1275.         _RealType __sum = 0;
  1276.  
  1277.         do
  1278.           {
  1279.             const _RealType __e = -std::log(__urng());
  1280.             __sum += __e / (__t - __x);
  1281.             __x += 1;
  1282.           }
  1283.         while (__sum <= _M_q);
  1284.  
  1285.         return __x - 1;
  1286.       }
  1287.  
  1288.   /**
  1289.    * A rejection algorithm when t * p >= 8 and a simple waiting time
  1290.    * method - the second in the referenced book - otherwise.
  1291.    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
  1292.    * is defined.
  1293.    *
  1294.    * Reference:
  1295.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
  1296.    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
  1297.    */
  1298.   template<typename _IntType, typename _RealType>
  1299.     template<class _UniformRandomNumberGenerator>
  1300.       typename binomial_distribution<_IntType, _RealType>::result_type
  1301.       binomial_distribution<_IntType, _RealType>::
  1302.       operator()(_UniformRandomNumberGenerator& __urng)
  1303.       {
  1304.         result_type __ret;
  1305.         const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
  1306.  
  1307. #if _GLIBCXX_USE_C99_MATH_TR1
  1308.         if (!_M_easy)
  1309.           {
  1310.             _RealType __x;
  1311.  
  1312.             // See comments above...
  1313.             const _RealType __naf =
  1314.               (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
  1315.             const _RealType __thr =
  1316.               std::numeric_limits<_IntType>::max() + __naf;
  1317.  
  1318.             const _RealType __np = std::floor(_M_t * __p12);
  1319.             const _RealType __pa = __np / _M_t;
  1320.  
  1321.             // sqrt(pi / 2)
  1322.             const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
  1323.             const _RealType __a1 = _M_a1;
  1324.             const _RealType __a12 = __a1 + _M_s2 * __spi_2;
  1325.             const _RealType __a123 = _M_a123;
  1326.             const _RealType __s1s = _M_s1 * _M_s1;
  1327.             const _RealType __s2s = _M_s2 * _M_s2;
  1328.  
  1329.             bool __reject;
  1330.             do
  1331.               {
  1332.                 const _RealType __u = _M_s * __urng();
  1333.  
  1334.                 _RealType __v;
  1335.  
  1336.                 if (__u <= __a1)
  1337.                   {
  1338.                     const _RealType __n = _M_nd(__urng);
  1339.                     const _RealType __y = _M_s1 * std::abs(__n);
  1340.                     __reject = __y >= _M_d1;
  1341.                     if (!__reject)
  1342.                       {
  1343.                         const _RealType __e = -std::log(__urng());
  1344.                         __x = std::floor(__y);
  1345.                         __v = -__e - __n * __n / 2 + _M_c;
  1346.                       }
  1347.                   }
  1348.                 else if (__u <= __a12)
  1349.                   {
  1350.                     const _RealType __n = _M_nd(__urng);
  1351.                     const _RealType __y = _M_s2 * std::abs(__n);
  1352.                     __reject = __y >= _M_d2;
  1353.                     if (!__reject)
  1354.                       {
  1355.                         const _RealType __e = -std::log(__urng());
  1356.                         __x = std::floor(-__y);
  1357.                         __v = -__e - __n * __n / 2;
  1358.                       }
  1359.                   }
  1360.                 else if (__u <= __a123)
  1361.                   {
  1362.                     const _RealType __e1 = -std::log(__urng());            
  1363.                     const _RealType __e2 = -std::log(__urng());
  1364.  
  1365.                     const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
  1366.                     __x = std::floor(__y);
  1367.                     __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
  1368.                                             -__y / (2 * __s1s)));
  1369.                     __reject = false;
  1370.                   }
  1371.                 else
  1372.                   {
  1373.                     const _RealType __e1 = -std::log(__urng());            
  1374.                     const _RealType __e2 = -std::log(__urng());
  1375.  
  1376.                     const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
  1377.                     __x = std::floor(-__y);
  1378.                     __v = -__e2 - _M_d2 * __y / (2 * __s2s);
  1379.                     __reject = false;
  1380.                   }
  1381.  
  1382.                 __reject = __reject || __x < -__np || __x > _M_t - __np;
  1383.                 if (!__reject)
  1384.                   {
  1385.                     const _RealType __lfx =
  1386.                       std::tr1::lgamma(__np + __x + 1)
  1387.                       + std::tr1::lgamma(_M_t - (__np + __x) + 1);
  1388.                     __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
  1389.                   }
  1390.  
  1391.                 __reject |= __x + __np >= __thr;
  1392.               }
  1393.             while (__reject);
  1394.  
  1395.             __x += __np + __naf;
  1396.  
  1397.             const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
  1398.             __ret = _IntType(__x) + __z;
  1399.           }
  1400.         else
  1401. #endif
  1402.           __ret = _M_waiting(__urng, _M_t);
  1403.  
  1404.         if (__p12 != _M_p)
  1405.           __ret = _M_t - __ret;
  1406.         return __ret;
  1407.       }
  1408.  
  1409.   template<typename _IntType, typename _RealType,
  1410.            typename _CharT, typename _Traits>
  1411.     std::basic_ostream<_CharT, _Traits>&
  1412.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1413.                const binomial_distribution<_IntType, _RealType>& __x)
  1414.     {
  1415.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1416.       typedef typename __ostream_type::ios_base    __ios_base;
  1417.  
  1418.       const typename __ios_base::fmtflags __flags = __os.flags();
  1419.       const _CharT __fill = __os.fill();
  1420.       const std::streamsize __precision = __os.precision();
  1421.       const _CharT __space = __os.widen(' ');
  1422.       __os.flags(__ios_base::scientific | __ios_base::left);
  1423.       __os.fill(__space);
  1424.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1425.  
  1426.       __os << __x.t() << __space << __x.p()
  1427.            << __space << __x._M_nd;
  1428.  
  1429.       __os.flags(__flags);
  1430.       __os.fill(__fill);
  1431.       __os.precision(__precision);
  1432.       return __os;
  1433.     }
  1434.  
  1435.   template<typename _IntType, typename _RealType,
  1436.            typename _CharT, typename _Traits>
  1437.     std::basic_istream<_CharT, _Traits>&
  1438.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  1439.                binomial_distribution<_IntType, _RealType>& __x)
  1440.     {
  1441.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  1442.       typedef typename __istream_type::ios_base    __ios_base;
  1443.  
  1444.       const typename __ios_base::fmtflags __flags = __is.flags();
  1445.       __is.flags(__ios_base::dec | __ios_base::skipws);
  1446.  
  1447.       __is >> __x._M_t >> __x._M_p >> __x._M_nd;
  1448.       __x._M_initialize();
  1449.  
  1450.       __is.flags(__flags);
  1451.       return __is;
  1452.     }
  1453.  
  1454.  
  1455.   template<typename _RealType, typename _CharT, typename _Traits>
  1456.     std::basic_ostream<_CharT, _Traits>&
  1457.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1458.                const uniform_real<_RealType>& __x)
  1459.     {
  1460.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1461.       typedef typename __ostream_type::ios_base    __ios_base;
  1462.  
  1463.       const typename __ios_base::fmtflags __flags = __os.flags();
  1464.       const _CharT __fill = __os.fill();
  1465.       const std::streamsize __precision = __os.precision();
  1466.       const _CharT __space = __os.widen(' ');
  1467.       __os.flags(__ios_base::scientific | __ios_base::left);
  1468.       __os.fill(__space);
  1469.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1470.  
  1471.       __os << __x.min() << __space << __x.max();
  1472.  
  1473.       __os.flags(__flags);
  1474.       __os.fill(__fill);
  1475.       __os.precision(__precision);
  1476.       return __os;
  1477.     }
  1478.  
  1479.   template<typename _RealType, typename _CharT, typename _Traits>
  1480.     std::basic_istream<_CharT, _Traits>&
  1481.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  1482.                uniform_real<_RealType>& __x)
  1483.     {
  1484.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  1485.       typedef typename __istream_type::ios_base    __ios_base;
  1486.  
  1487.       const typename __ios_base::fmtflags __flags = __is.flags();
  1488.       __is.flags(__ios_base::skipws);
  1489.  
  1490.       __is >> __x._M_min >> __x._M_max;
  1491.  
  1492.       __is.flags(__flags);
  1493.       return __is;
  1494.     }
  1495.  
  1496.  
  1497.   template<typename _RealType, typename _CharT, typename _Traits>
  1498.     std::basic_ostream<_CharT, _Traits>&
  1499.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1500.                const exponential_distribution<_RealType>& __x)
  1501.     {
  1502.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1503.       typedef typename __ostream_type::ios_base    __ios_base;
  1504.  
  1505.       const typename __ios_base::fmtflags __flags = __os.flags();
  1506.       const _CharT __fill = __os.fill();
  1507.       const std::streamsize __precision = __os.precision();
  1508.       __os.flags(__ios_base::scientific | __ios_base::left);
  1509.       __os.fill(__os.widen(' '));
  1510.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1511.  
  1512.       __os << __x.lambda();
  1513.  
  1514.       __os.flags(__flags);
  1515.       __os.fill(__fill);
  1516.       __os.precision(__precision);
  1517.       return __os;
  1518.     }
  1519.  
  1520.  
  1521.   /**
  1522.    * Polar method due to Marsaglia.
  1523.    *
  1524.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
  1525.    * New York, 1986, Ch. V, Sect. 4.4.
  1526.    */
  1527.   template<typename _RealType>
  1528.     template<class _UniformRandomNumberGenerator>
  1529.       typename normal_distribution<_RealType>::result_type
  1530.       normal_distribution<_RealType>::
  1531.       operator()(_UniformRandomNumberGenerator& __urng)
  1532.       {
  1533.         result_type __ret;
  1534.  
  1535.         if (_M_saved_available)
  1536.           {
  1537.             _M_saved_available = false;
  1538.             __ret = _M_saved;
  1539.           }
  1540.         else
  1541.           {
  1542.             result_type __x, __y, __r2;
  1543.             do
  1544.               {
  1545.                 __x = result_type(2.0) * __urng() - 1.0;
  1546.                 __y = result_type(2.0) * __urng() - 1.0;
  1547.                 __r2 = __x * __x + __y * __y;
  1548.               }
  1549.             while (__r2 > 1.0 || __r2 == 0.0);
  1550.  
  1551.             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
  1552.             _M_saved = __x * __mult;
  1553.             _M_saved_available = true;
  1554.             __ret = __y * __mult;
  1555.           }
  1556.        
  1557.         __ret = __ret * _M_sigma + _M_mean;
  1558.         return __ret;
  1559.       }
  1560.  
  1561.   template<typename _RealType, typename _CharT, typename _Traits>
  1562.     std::basic_ostream<_CharT, _Traits>&
  1563.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1564.                const normal_distribution<_RealType>& __x)
  1565.     {
  1566.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1567.       typedef typename __ostream_type::ios_base    __ios_base;
  1568.  
  1569.       const typename __ios_base::fmtflags __flags = __os.flags();
  1570.       const _CharT __fill = __os.fill();
  1571.       const std::streamsize __precision = __os.precision();
  1572.       const _CharT __space = __os.widen(' ');
  1573.       __os.flags(__ios_base::scientific | __ios_base::left);
  1574.       __os.fill(__space);
  1575.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1576.  
  1577.       __os << __x._M_saved_available << __space
  1578.            << __x.mean() << __space
  1579.            << __x.sigma();
  1580.       if (__x._M_saved_available)
  1581.         __os << __space << __x._M_saved;
  1582.  
  1583.       __os.flags(__flags);
  1584.       __os.fill(__fill);
  1585.       __os.precision(__precision);
  1586.       return __os;
  1587.     }
  1588.  
  1589.   template<typename _RealType, typename _CharT, typename _Traits>
  1590.     std::basic_istream<_CharT, _Traits>&
  1591.     operator>>(std::basic_istream<_CharT, _Traits>& __is,
  1592.                normal_distribution<_RealType>& __x)
  1593.     {
  1594.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
  1595.       typedef typename __istream_type::ios_base    __ios_base;
  1596.  
  1597.       const typename __ios_base::fmtflags __flags = __is.flags();
  1598.       __is.flags(__ios_base::dec | __ios_base::skipws);
  1599.  
  1600.       __is >> __x._M_saved_available >> __x._M_mean
  1601.            >> __x._M_sigma;
  1602.       if (__x._M_saved_available)
  1603.         __is >> __x._M_saved;
  1604.  
  1605.       __is.flags(__flags);
  1606.       return __is;
  1607.     }
  1608.  
  1609.  
  1610.   template<typename _RealType>
  1611.     void
  1612.     gamma_distribution<_RealType>::
  1613.     _M_initialize()
  1614.     {
  1615.       if (_M_alpha >= 1)
  1616.         _M_l_d = std::sqrt(2 * _M_alpha - 1);
  1617.       else
  1618.         _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
  1619.                   * (1 - _M_alpha));
  1620.     }
  1621.  
  1622.   /**
  1623.    * Cheng's rejection algorithm GB for alpha >= 1 and a modification
  1624.    * of Vaduva's rejection from Weibull algorithm due to Devroye for
  1625.    * alpha < 1.
  1626.    *
  1627.    * References:
  1628.    * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral
  1629.    * Shape Parameter. Applied Statistics, 26, 71-75, 1977.
  1630.    *
  1631.    * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection
  1632.    * and Composition Procedures. Math. Operationsforschung and Statistik,
  1633.    * Series in Statistics, 8, 545-576, 1977.
  1634.    *
  1635.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
  1636.    * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
  1637.    */
  1638.   template<typename _RealType>
  1639.     template<class _UniformRandomNumberGenerator>
  1640.       typename gamma_distribution<_RealType>::result_type
  1641.       gamma_distribution<_RealType>::
  1642.       operator()(_UniformRandomNumberGenerator& __urng)
  1643.       {
  1644.         result_type __x;
  1645.  
  1646.         bool __reject;
  1647.         if (_M_alpha >= 1)
  1648.           {
  1649.             // alpha - log(4)
  1650.             const result_type __b = _M_alpha
  1651.               - result_type(1.3862943611198906188344642429163531L);
  1652.             const result_type __c = _M_alpha + _M_l_d;
  1653.             const result_type __1l = 1 / _M_l_d;
  1654.  
  1655.             // 1 + log(9 / 2)
  1656.             const result_type __k = 2.5040773967762740733732583523868748L;
  1657.  
  1658.             do
  1659.               {
  1660.                 const result_type __u = __urng();
  1661.                 const result_type __v = __urng();
  1662.  
  1663.                 const result_type __y = __1l * std::log(__v / (1 - __v));
  1664.                 __x = _M_alpha * std::exp(__y);
  1665.  
  1666.                 const result_type __z = __u * __v * __v;
  1667.                 const result_type __r = __b + __c * __y - __x;
  1668.  
  1669.                 __reject = __r < result_type(4.5) * __z - __k;
  1670.                 if (__reject)
  1671.                   __reject = __r < std::log(__z);
  1672.               }
  1673.             while (__reject);
  1674.           }
  1675.         else
  1676.           {
  1677.             const result_type __c = 1 / _M_alpha;
  1678.  
  1679.             do
  1680.               {
  1681.                 const result_type __z = -std::log(__urng());
  1682.                 const result_type __e = -std::log(__urng());
  1683.  
  1684.                 __x = std::pow(__z, __c);
  1685.  
  1686.                 __reject = __z + __e < _M_l_d + __x;
  1687.               }
  1688.             while (__reject);
  1689.           }
  1690.  
  1691.         return __x;
  1692.       }
  1693.  
  1694.   template<typename _RealType, typename _CharT, typename _Traits>
  1695.     std::basic_ostream<_CharT, _Traits>&
  1696.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
  1697.                const gamma_distribution<_RealType>& __x)
  1698.     {
  1699.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
  1700.       typedef typename __ostream_type::ios_base    __ios_base;
  1701.  
  1702.       const typename __ios_base::fmtflags __flags = __os.flags();
  1703.       const _CharT __fill = __os.fill();
  1704.       const std::streamsize __precision = __os.precision();
  1705.       __os.flags(__ios_base::scientific | __ios_base::left);
  1706.       __os.fill(__os.widen(' '));
  1707.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
  1708.  
  1709.       __os << __x.alpha();
  1710.  
  1711.       __os.flags(__flags);
  1712.       __os.fill(__fill);
  1713.       __os.precision(__precision);
  1714.       return __os;
  1715.     }
  1716.  
  1717. _GLIBCXX_END_NAMESPACE_VERSION
  1718. }
  1719. }
  1720.  
  1721. #endif
  1722.