visus
/
toolchain-gcc-arm-embedded
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2 Commits
1 Branch
Branch: main
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'main'
${ noResults }
toolchain-gcc-arm-embedded/arm-none-eabi/include/c++/10.2.1/tr1/random.tcc

1719 lines
53KB
Raw Permalink Blame History 

  1. // random number generation (out of line) -*- C++ -*-

  2. 

  3. // Copyright (C) 2009-2020 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. _GLIBCXX_BEGIN_NAMESPACE_VERSION

  37. 

  38. namespace tr1

  39. {

  40.   /*

  41.    * (Further) implementation-space details.

  42.    */

  43.   namespace __detail

  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.   } // namespace __detail

  96. 

  97.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  98.     const _UIntType

  99.     linear_congruential<_UIntType, __a, __c, __m>::multiplier;

  100. 

  101.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  102.     const _UIntType

  103.     linear_congruential<_UIntType, __a, __c, __m>::increment;

  104. 

  105.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  106.     const _UIntType

  107.     linear_congruential<_UIntType, __a, __c, __m>::modulus;

  108. 

  109.   /**

  110.    * Seeds the LCR with integral value @p __x0, adjusted so that the 

  111.    * ring identity is never a member of the convergence set.

  112.    */

  113.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  114.     void

  115.     linear_congruential<_UIntType, __a, __c, __m>::

  116.     seed(unsigned long __x0)

  117.     {

  118.       if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)

  119. 	  && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))

  120. 	_M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);

  121.       else

  122. 	_M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);

  123.     }

  124. 

  125.   /**

  126.    * Seeds the LCR engine with a value generated by @p __g.

  127.    */

  128.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  129.     template<class _Gen>

  130.       void

  131.       linear_congruential<_UIntType, __a, __c, __m>::

  132.       seed(_Gen& __g, false_type)

  133.       {

  134. 	_UIntType __x0 = __g();

  135. 	if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)

  136. 	    && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))

  137. 	  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);

  138. 	else

  139. 	  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);

  140.       }

  141. 

  142.   /**

  143.    * Gets the next generated value in sequence.

  144.    */

  145.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>

  146.     typename linear_congruential<_UIntType, __a, __c, __m>::result_type

  147.     linear_congruential<_UIntType, __a, __c, __m>::

  148.     operator()()

  149.     {

  150.       _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);

  151.       return _M_x;

  152.     }

  153. 

  154.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,

  155. 	   typename _CharT, typename _Traits>

  156.     std::basic_ostream<_CharT, _Traits>&

  157.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  158. 	       const linear_congruential<_UIntType, __a, __c, __m>& __lcr)

  159.     {

  160.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  161.       typedef typename __ostream_type::ios_base    __ios_base;

  162. 

  163.       const typename __ios_base::fmtflags __flags = __os.flags();

  164.       const _CharT __fill = __os.fill();

  165.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);

  166.       __os.fill(__os.widen(' '));

  167. 

  168.       __os << __lcr._M_x;

  169. 

  170.       __os.flags(__flags);

  171.       __os.fill(__fill);

  172.       return __os;

  173.     }

  174. 

  175.   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,

  176. 	   typename _CharT, typename _Traits>

  177.     std::basic_istream<_CharT, _Traits>&

  178.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  179. 	       linear_congruential<_UIntType, __a, __c, __m>& __lcr)

  180.     {

  181.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  182.       typedef typename __istream_type::ios_base    __ios_base;

  183. 

  184.       const typename __ios_base::fmtflags __flags = __is.flags();

  185.       __is.flags(__ios_base::dec);

  186. 

  187.       __is >> __lcr._M_x;

  188. 

  189.       __is.flags(__flags);

  190.       return __is;

  191.     } 

  192. 

  193. 

  194.   template<class _UIntType, int __w, int __n, int __m, int __r,

  195. 	   _UIntType __a, int __u, int __s,

  196. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  197.     const int

  198.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  199. 		     __b, __t, __c, __l>::word_size;

  200. 

  201.   template<class _UIntType, int __w, int __n, int __m, int __r,

  202. 	   _UIntType __a, int __u, int __s,

  203. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  204.     const int

  205.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  206. 		     __b, __t, __c, __l>::state_size;

  207.     

  208.   template<class _UIntType, int __w, int __n, int __m, int __r,

  209. 	   _UIntType __a, int __u, int __s,

  210. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  211.     const int

  212.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  213. 		     __b, __t, __c, __l>::shift_size;

  214. 

  215.   template<class _UIntType, int __w, int __n, int __m, int __r,

  216. 	   _UIntType __a, int __u, int __s,

  217. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  218.     const int

  219.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  220. 		     __b, __t, __c, __l>::mask_bits;

  221. 

  222.   template<class _UIntType, int __w, int __n, int __m, int __r,

  223. 	   _UIntType __a, int __u, int __s,

  224. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  225.     const _UIntType

  226.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  227. 		     __b, __t, __c, __l>::parameter_a;

  228. 

  229.   template<class _UIntType, int __w, int __n, int __m, int __r,

  230. 	   _UIntType __a, int __u, int __s,

  231. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  232.     const int

  233.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  234. 		     __b, __t, __c, __l>::output_u;

  235. 

  236.   template<class _UIntType, int __w, int __n, int __m, int __r,

  237. 	   _UIntType __a, int __u, int __s,

  238. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  239.     const int

  240.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  241. 		     __b, __t, __c, __l>::output_s;

  242. 

  243.   template<class _UIntType, int __w, int __n, int __m, int __r,

  244. 	   _UIntType __a, int __u, int __s,

  245. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  246.     const _UIntType

  247.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  248. 		     __b, __t, __c, __l>::output_b;

  249. 

  250.   template<class _UIntType, int __w, int __n, int __m, int __r,

  251. 	   _UIntType __a, int __u, int __s,

  252. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  253.     const int

  254.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  255. 		     __b, __t, __c, __l>::output_t;

  256. 

  257.   template<class _UIntType, int __w, int __n, int __m, int __r,

  258. 	   _UIntType __a, int __u, int __s,

  259. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  260.     const _UIntType

  261.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  262. 		     __b, __t, __c, __l>::output_c;

  263. 

  264.   template<class _UIntType, int __w, int __n, int __m, int __r,

  265. 	   _UIntType __a, int __u, int __s,

  266. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  267.     const int

  268.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  269. 		     __b, __t, __c, __l>::output_l;

  270. 

  271.   template<class _UIntType, int __w, int __n, int __m, int __r,

  272. 	   _UIntType __a, int __u, int __s,

  273. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  274.     void

  275.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  276. 		     __b, __t, __c, __l>::

  277.     seed(unsigned long __value)

  278.     {

  279.       _M_x[0] = __detail::__mod<_UIntType, 1, 0,

  280. 	__detail::_Shift<_UIntType, __w>::__value>(__value);

  281. 

  282.       for (int __i = 1; __i < state_size; ++__i)

  283. 	{

  284. 	  _UIntType __x = _M_x[__i - 1];

  285. 	  __x ^= __x >> (__w - 2);

  286. 	  __x *= 1812433253ul;

  287. 	  __x += __i;

  288. 	  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,

  289. 	    __detail::_Shift<_UIntType, __w>::__value>(__x);	  

  290. 	}

  291.       _M_p = state_size;

  292.     }

  293. 

  294.   template<class _UIntType, int __w, int __n, int __m, int __r,

  295. 	   _UIntType __a, int __u, int __s,

  296. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  297.     template<class _Gen>

  298.       void

  299.       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  300. 		       __b, __t, __c, __l>::

  301.       seed(_Gen& __gen, false_type)

  302.       {

  303. 	for (int __i = 0; __i < state_size; ++__i)

  304. 	  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,

  305. 	    __detail::_Shift<_UIntType, __w>::__value>(__gen());

  306. 	_M_p = state_size;

  307.       }

  308. 

  309.   template<class _UIntType, int __w, int __n, int __m, int __r,

  310. 	   _UIntType __a, int __u, int __s,

  311. 	   _UIntType __b, int __t, _UIntType __c, int __l>

  312.     typename

  313.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  314. 		     __b, __t, __c, __l>::result_type

  315.     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,

  316. 		     __b, __t, __c, __l>::

  317.     operator()()

  318.     {

  319.       // Reload the vector - cost is O(n) amortized over n calls.

  320.       if (_M_p >= state_size)

  321. 	{

  322. 	  const _UIntType __upper_mask = (~_UIntType()) << __r;

  323. 	  const _UIntType __lower_mask = ~__upper_mask;

  324. 

  325. 	  for (int __k = 0; __k < (__n - __m); ++__k)

  326. 	    {

  327. 	      _UIntType __y = ((_M_x[__k] & __upper_mask)

  328. 			       | (_M_x[__k + 1] & __lower_mask));

  329. 	      _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)

  330. 			   ^ ((__y & 0x01) ? __a : 0));

  331. 	    }

  332. 

  333. 	  for (int __k = (__n - __m); __k < (__n - 1); ++__k)

  334. 	    {

  335. 	      _UIntType __y = ((_M_x[__k] & __upper_mask)

  336. 			       | (_M_x[__k + 1] & __lower_mask));

  337. 	      _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)

  338. 			   ^ ((__y & 0x01) ? __a : 0));

  339. 	    }

  340. 

  341. 	  _UIntType __y = ((_M_x[__n - 1] & __upper_mask)

  342. 			   | (_M_x[0] & __lower_mask));

  343. 	  _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)

  344. 			   ^ ((__y & 0x01) ? __a : 0));

  345. 	  _M_p = 0;

  346. 	}

  347. 

  348.       // Calculate o(x(i)).

  349.       result_type __z = _M_x[_M_p++];

  350.       __z ^= (__z >> __u);

  351.       __z ^= (__z << __s) & __b;

  352.       __z ^= (__z << __t) & __c;

  353.       __z ^= (__z >> __l);

  354. 

  355.       return __z;

  356.     }

  357. 

  358.   template<class _UIntType, int __w, int __n, int __m, int __r,

  359. 	   _UIntType __a, int __u, int __s, _UIntType __b, int __t,

  360. 	   _UIntType __c, int __l,

  361. 	   typename _CharT, typename _Traits>

  362.     std::basic_ostream<_CharT, _Traits>&

  363.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  364. 	       const mersenne_twister<_UIntType, __w, __n, __m,

  365. 	       __r, __a, __u, __s, __b, __t, __c, __l>& __x)

  366.     {

  367.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  368.       typedef typename __ostream_type::ios_base    __ios_base;

  369. 

  370.       const typename __ios_base::fmtflags __flags = __os.flags();

  371.       const _CharT __fill = __os.fill();

  372.       const _CharT __space = __os.widen(' ');

  373.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);

  374.       __os.fill(__space);

  375. 

  376.       for (int __i = 0; __i < __n - 1; ++__i)

  377. 	__os << __x._M_x[__i] << __space;

  378.       __os << __x._M_x[__n - 1];

  379. 

  380.       __os.flags(__flags);

  381.       __os.fill(__fill);

  382.       return __os;

  383.     }

  384. 

  385.   template<class _UIntType, int __w, int __n, int __m, int __r,

  386. 	   _UIntType __a, int __u, int __s, _UIntType __b, int __t,

  387. 	   _UIntType __c, int __l,

  388. 	   typename _CharT, typename _Traits>

  389.     std::basic_istream<_CharT, _Traits>&

  390.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  391. 	       mersenne_twister<_UIntType, __w, __n, __m,

  392. 	       __r, __a, __u, __s, __b, __t, __c, __l>& __x)

  393.     {

  394.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  395.       typedef typename __istream_type::ios_base    __ios_base;

  396. 

  397.       const typename __ios_base::fmtflags __flags = __is.flags();

  398.       __is.flags(__ios_base::dec | __ios_base::skipws);

  399. 

  400.       for (int __i = 0; __i < __n; ++__i)

  401. 	__is >> __x._M_x[__i];

  402. 

  403.       __is.flags(__flags);

  404.       return __is;

  405.     }

  406. 

  407. 

  408.   template<typename _IntType, _IntType __m, int __s, int __r>

  409.     const _IntType

  410.     subtract_with_carry<_IntType, __m, __s, __r>::modulus;

  411. 

  412.   template<typename _IntType, _IntType __m, int __s, int __r>

  413.     const int

  414.     subtract_with_carry<_IntType, __m, __s, __r>::long_lag;

  415. 

  416.   template<typename _IntType, _IntType __m, int __s, int __r>

  417.     const int

  418.     subtract_with_carry<_IntType, __m, __s, __r>::short_lag;

  419. 

  420.   template<typename _IntType, _IntType __m, int __s, int __r>

  421.     void

  422.     subtract_with_carry<_IntType, __m, __s, __r>::

  423.     seed(unsigned long __value)

  424.     {

  425.       if (__value == 0)

  426. 	__value = 19780503;

  427. 

  428.       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>

  429. 	__lcg(__value);

  430. 

  431.       for (int __i = 0; __i < long_lag; ++__i)

  432. 	_M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());

  433. 

  434.       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;

  435.       _M_p = 0;

  436.     }

  437. 

  438.   template<typename _IntType, _IntType __m, int __s, int __r>

  439.     template<class _Gen>

  440.       void

  441.       subtract_with_carry<_IntType, __m, __s, __r>::

  442.       seed(_Gen& __gen, false_type)

  443.       {

  444. 	const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;

  445. 

  446. 	for (int __i = 0; __i < long_lag; ++__i)

  447. 	  {

  448. 	    _UIntType __tmp = 0;

  449. 	    _UIntType __factor = 1;

  450. 	    for (int __j = 0; __j < __n; ++__j)

  451. 	      {

  452. 		__tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>

  453. 		         (__gen()) * __factor;

  454. 		__factor *= __detail::_Shift<_UIntType, 32>::__value;

  455. 	      }

  456. 	    _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);

  457. 	  }

  458. 	_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;

  459. 	_M_p = 0;

  460.       }

  461. 

  462.   template<typename _IntType, _IntType __m, int __s, int __r>

  463.     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type

  464.     subtract_with_carry<_IntType, __m, __s, __r>::

  465.     operator()()

  466.     {

  467.       // Derive short lag index from current index.

  468.       int __ps = _M_p - short_lag;

  469.       if (__ps < 0)

  470. 	__ps += long_lag;

  471. 

  472.       // Calculate new x(i) without overflow or division.

  473.       // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry

  474.       // cannot overflow.

  475.       _UIntType __xi;

  476.       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)

  477. 	{

  478. 	  __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;

  479. 	  _M_carry = 0;

  480. 	}

  481.       else

  482. 	{

  483. 	  __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];

  484. 	  _M_carry = 1;

  485. 	}

  486.       _M_x[_M_p] = __xi;

  487. 

  488.       // Adjust current index to loop around in ring buffer.

  489.       if (++_M_p >= long_lag)

  490. 	_M_p = 0;

  491. 

  492.       return __xi;

  493.     }

  494. 

  495.   template<typename _IntType, _IntType __m, int __s, int __r,

  496. 	   typename _CharT, typename _Traits>

  497.     std::basic_ostream<_CharT, _Traits>&

  498.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  499. 	       const subtract_with_carry<_IntType, __m, __s, __r>& __x)

  500.     {

  501.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  502.       typedef typename __ostream_type::ios_base    __ios_base;

  503. 

  504.       const typename __ios_base::fmtflags __flags = __os.flags();

  505.       const _CharT __fill = __os.fill();

  506.       const _CharT __space = __os.widen(' ');

  507.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);

  508.       __os.fill(__space);

  509. 

  510.       for (int __i = 0; __i < __r; ++__i)

  511. 	__os << __x._M_x[__i] << __space;

  512.       __os << __x._M_carry;

  513. 

  514.       __os.flags(__flags);

  515.       __os.fill(__fill);

  516.       return __os;

  517.     }

  518. 

  519.   template<typename _IntType, _IntType __m, int __s, int __r,

  520. 	   typename _CharT, typename _Traits>

  521.     std::basic_istream<_CharT, _Traits>&

  522.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  523. 	       subtract_with_carry<_IntType, __m, __s, __r>& __x)

  524.     {

  525.       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;

  526.       typedef typename __istream_type::ios_base    __ios_base;

  527. 

  528.       const typename __ios_base::fmtflags __flags = __is.flags();

  529.       __is.flags(__ios_base::dec | __ios_base::skipws);

  530. 

  531.       for (int __i = 0; __i < __r; ++__i)

  532. 	__is >> __x._M_x[__i];

  533.       __is >> __x._M_carry;

  534. 

  535.       __is.flags(__flags);

  536.       return __is;

  537.     }

  538. 

  539. 

  540.   template<typename _RealType, int __w, int __s, int __r>

  541.     const int

  542.     subtract_with_carry_01<_RealType, __w, __s, __r>::word_size;

  543. 

  544.   template<typename _RealType, int __w, int __s, int __r>

  545.     const int

  546.     subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag;

  547. 

  548.   template<typename _RealType, int __w, int __s, int __r>

  549.     const int

  550.     subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag;

  551. 

  552.   template<typename _RealType, int __w, int __s, int __r>

  553.     void

  554.     subtract_with_carry_01<_RealType, __w, __s, __r>::

  555.     _M_initialize_npows()

  556.     {

  557.       for (int __j = 0; __j < __n; ++__j)

  558. #if _GLIBCXX_USE_C99_MATH_TR1

  559. 	_M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);

  560. #else

  561.         _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);

  562. #endif

  563.     }

  564. 

  565.   template<typename _RealType, int __w, int __s, int __r>

  566.     void

  567.     subtract_with_carry_01<_RealType, __w, __s, __r>::

  568.     seed(unsigned long __value)

  569.     {

  570.       if (__value == 0)

  571. 	__value = 19780503;

  572. 

  573.       // _GLIBCXX_RESOLVE_LIB_DEFECTS

  574.       // 512. Seeding subtract_with_carry_01 from a single unsigned long.

  575.       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>

  576. 	__lcg(__value);

  577. 

  578.       this->seed(__lcg);

  579.     }

  580. 

  581.   template<typename _RealType, int __w, int __s, int __r>

  582.     template<class _Gen>

  583.       void

  584.       subtract_with_carry_01<_RealType, __w, __s, __r>::

  585.       seed(_Gen& __gen, false_type)

  586.       {

  587. 	for (int __i = 0; __i < long_lag; ++__i)

  588. 	  {

  589. 	    for (int __j = 0; __j < __n - 1; ++__j)

  590. 	      _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());

  591. 	    _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,

  592. 	      __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());

  593. 	  }

  594. 

  595. 	_M_carry = 1;

  596. 	for (int __j = 0; __j < __n; ++__j)

  597. 	  if (_M_x[long_lag - 1][__j] != 0)

  598. 	    {

  599. 	      _M_carry = 0;

  600. 	      break;

  601. 	    }

  602. 

  603. 	_M_p = 0;

  604.       }

  605. 

  606.   template<typename _RealType, int __w, int __s, int __r>

  607.     typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type

  608.     subtract_with_carry_01<_RealType, __w, __s, __r>::

  609.     operator()()

  610.     {

  611.       // Derive short lag index from current index.

  612.       int __ps = _M_p - short_lag;

  613.       if (__ps < 0)

  614. 	__ps += long_lag;

  615. 

  616.       _UInt32Type __new_carry;

  617.       for (int __j = 0; __j < __n - 1; ++__j)

  618. 	{

  619. 	  if (_M_x[__ps][__j] > _M_x[_M_p][__j]

  620. 	      || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))

  621. 	    __new_carry = 0;

  622. 	  else

  623. 	    __new_carry = 1;

  624. 

  625. 	  _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;

  626. 	  _M_carry = __new_carry;

  627. 	}

  628. 

  629.       if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]

  630. 	  || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))

  631. 	__new_carry = 0;

  632.       else

  633. 	__new_carry = 1;

  634.       

  635.       _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,

  636. 	__detail::_Shift<_UInt32Type, __w % 32>::__value>

  637. 	(_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);

  638.       _M_carry = __new_carry;

  639. 

  640.       result_type __ret = 0.0;

  641.       for (int __j = 0; __j < __n; ++__j)

  642. 	__ret += _M_x[_M_p][__j] * _M_npows[__j];

  643. 

  644.       // Adjust current index to loop around in ring buffer.

  645.       if (++_M_p >= long_lag)

  646. 	_M_p = 0;

  647. 

  648.       return __ret;

  649.     }

  650. 

  651.   template<typename _RealType, int __w, int __s, int __r,

  652. 	   typename _CharT, typename _Traits>

  653.     std::basic_ostream<_CharT, _Traits>&

  654.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  655. 	       const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)

  656.     {

  657.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  658.       typedef typename __ostream_type::ios_base    __ios_base;

  659. 

  660.       const typename __ios_base::fmtflags __flags = __os.flags();

  661.       const _CharT __fill = __os.fill();

  662.       const _CharT __space = __os.widen(' ');

  663.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);

  664.       __os.fill(__space);

  665. 

  666.       for (int __i = 0; __i < __r; ++__i)

  667. 	for (int __j = 0; __j < __x.__n; ++__j)

  668. 	  __os << __x._M_x[__i][__j] << __space;

  669.       __os << __x._M_carry;

  670. 

  671.       __os.flags(__flags);

  672.       __os.fill(__fill);

  673.       return __os;

  674.     }

  675. 

  676.   template<typename _RealType, int __w, int __s, int __r,

  677. 	   typename _CharT, typename _Traits>

  678.     std::basic_istream<_CharT, _Traits>&

  679.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  680. 	       subtract_with_carry_01<_RealType, __w, __s, __r>& __x)

  681.     {

  682.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  683.       typedef typename __istream_type::ios_base    __ios_base;

  684. 

  685.       const typename __ios_base::fmtflags __flags = __is.flags();

  686.       __is.flags(__ios_base::dec | __ios_base::skipws);

  687. 

  688.       for (int __i = 0; __i < __r; ++__i)

  689. 	for (int __j = 0; __j < __x.__n; ++__j)

  690. 	  __is >> __x._M_x[__i][__j];

  691.       __is >> __x._M_carry;

  692. 

  693.       __is.flags(__flags);

  694.       return __is;

  695.     }

  696. 

  697.   template<class _UniformRandomNumberGenerator, int __p, int __r>

  698.     const int

  699.     discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size;

  700. 

  701.   template<class _UniformRandomNumberGenerator, int __p, int __r>

  702.     const int

  703.     discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block;

  704. 

  705.   template<class _UniformRandomNumberGenerator, int __p, int __r>

  706.     typename discard_block<_UniformRandomNumberGenerator,

  707. 			   __p, __r>::result_type

  708.     discard_block<_UniformRandomNumberGenerator, __p, __r>::

  709.     operator()()

  710.     {

  711.       if (_M_n >= used_block)

  712. 	{

  713. 	  while (_M_n < block_size)

  714. 	    {

  715. 	      _M_b();

  716. 	      ++_M_n;

  717. 	    }

  718. 	  _M_n = 0;

  719. 	}

  720.       ++_M_n;

  721.       return _M_b();

  722.     }

  723. 

  724.   template<class _UniformRandomNumberGenerator, int __p, int __r,

  725. 	   typename _CharT, typename _Traits>

  726.     std::basic_ostream<_CharT, _Traits>&

  727.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  728. 	       const discard_block<_UniformRandomNumberGenerator,

  729. 	       __p, __r>& __x)

  730.     {

  731.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  732.       typedef typename __ostream_type::ios_base    __ios_base;

  733. 

  734.       const typename __ios_base::fmtflags __flags = __os.flags();

  735.       const _CharT __fill = __os.fill();

  736.       const _CharT __space = __os.widen(' ');

  737.       __os.flags(__ios_base::dec | __ios_base::fixed

  738. 		 | __ios_base::left);

  739.       __os.fill(__space);

  740. 

  741.       __os << __x._M_b << __space << __x._M_n;

  742. 

  743.       __os.flags(__flags);

  744.       __os.fill(__fill);

  745.       return __os;

  746.     }

  747. 

  748.   template<class _UniformRandomNumberGenerator, int __p, int __r,

  749. 	   typename _CharT, typename _Traits>

  750.     std::basic_istream<_CharT, _Traits>&

  751.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  752. 	       discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)

  753.     {

  754.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  755.       typedef typename __istream_type::ios_base    __ios_base;

  756. 

  757.       const typename __ios_base::fmtflags __flags = __is.flags();

  758.       __is.flags(__ios_base::dec | __ios_base::skipws);

  759. 

  760.       __is >> __x._M_b >> __x._M_n;

  761. 

  762.       __is.flags(__flags);

  763.       return __is;

  764.     }

  765. 

  766. 

  767.   template<class _UniformRandomNumberGenerator1, int __s1,

  768. 	   class _UniformRandomNumberGenerator2, int __s2>

  769.     const int

  770.     xor_combine<_UniformRandomNumberGenerator1, __s1,

  771. 		_UniformRandomNumberGenerator2, __s2>::shift1;

  772.      

  773.   template<class _UniformRandomNumberGenerator1, int __s1,

  774. 	   class _UniformRandomNumberGenerator2, int __s2>

  775.     const int

  776.     xor_combine<_UniformRandomNumberGenerator1, __s1,

  777. 		_UniformRandomNumberGenerator2, __s2>::shift2;

  778. 

  779.   template<class _UniformRandomNumberGenerator1, int __s1,

  780. 	   class _UniformRandomNumberGenerator2, int __s2>

  781.     void

  782.     xor_combine<_UniformRandomNumberGenerator1, __s1,

  783. 		_UniformRandomNumberGenerator2, __s2>::

  784.     _M_initialize_max()

  785.     {

  786.       const int __w = std::numeric_limits<result_type>::digits;

  787. 

  788.       const result_type __m1 =

  789. 	std::min(result_type(_M_b1.max() - _M_b1.min()),

  790. 		 __detail::_Shift<result_type, __w - __s1>::__value - 1);

  791. 

  792.       const result_type __m2 =

  793. 	std::min(result_type(_M_b2.max() - _M_b2.min()),

  794. 		 __detail::_Shift<result_type, __w - __s2>::__value - 1);

  795. 

  796.       // NB: In TR1 s1 is not required to be >= s2.

  797.       if (__s1 < __s2)

  798. 	_M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;

  799.       else

  800. 	_M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;

  801.     }

  802. 

  803.   template<class _UniformRandomNumberGenerator1, int __s1,

  804. 	   class _UniformRandomNumberGenerator2, int __s2>

  805.     typename xor_combine<_UniformRandomNumberGenerator1, __s1,

  806. 			 _UniformRandomNumberGenerator2, __s2>::result_type

  807.     xor_combine<_UniformRandomNumberGenerator1, __s1,

  808. 		_UniformRandomNumberGenerator2, __s2>::

  809.     _M_initialize_max_aux(result_type __a, result_type __b, int __d)

  810.     {

  811.       const result_type __two2d = result_type(1) << __d;

  812.       const result_type __c = __a * __two2d;

  813. 

  814.       if (__a == 0 || __b < __two2d)

  815. 	return __c + __b;

  816. 

  817.       const result_type __t = std::max(__c, __b);

  818.       const result_type __u = std::min(__c, __b);

  819. 

  820.       result_type __ub = __u;

  821.       result_type __p;

  822.       for (__p = 0; __ub != 1; __ub >>= 1)

  823. 	++__p;

  824. 

  825.       const result_type __two2p = result_type(1) << __p;

  826.       const result_type __k = __t / __two2p;

  827. 

  828.       if (__k & 1)

  829. 	return (__k + 1) * __two2p - 1;

  830. 

  831.       if (__c >= __b)

  832. 	return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)

  833. 							   / __two2d,

  834. 							   __u % __two2p, __d);

  835.       else

  836. 	return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)

  837. 							   / __two2d,

  838. 							   __t % __two2p, __d);

  839.     }

  840. 

  841.   template<class _UniformRandomNumberGenerator1, int __s1,

  842. 	   class _UniformRandomNumberGenerator2, int __s2,

  843. 	   typename _CharT, typename _Traits>

  844.     std::basic_ostream<_CharT, _Traits>&

  845.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  846. 	       const xor_combine<_UniformRandomNumberGenerator1, __s1,

  847. 	       _UniformRandomNumberGenerator2, __s2>& __x)

  848.     {

  849.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  850.       typedef typename __ostream_type::ios_base    __ios_base;

  851. 

  852.       const typename __ios_base::fmtflags __flags = __os.flags();

  853.       const _CharT __fill = __os.fill();

  854.       const _CharT __space = __os.widen(' ');

  855.       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);

  856.       __os.fill(__space);

  857. 

  858.       __os << __x.base1() << __space << __x.base2();

  859. 

  860.       __os.flags(__flags);

  861.       __os.fill(__fill);

  862.       return __os; 

  863.     }

  864. 

  865.   template<class _UniformRandomNumberGenerator1, int __s1,

  866. 	   class _UniformRandomNumberGenerator2, int __s2,

  867. 	   typename _CharT, typename _Traits>

  868.     std::basic_istream<_CharT, _Traits>&

  869.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  870. 	       xor_combine<_UniformRandomNumberGenerator1, __s1,

  871. 	       _UniformRandomNumberGenerator2, __s2>& __x)

  872.     {

  873.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  874.       typedef typename __istream_type::ios_base    __ios_base;

  875. 

  876.       const typename __ios_base::fmtflags __flags = __is.flags();

  877.       __is.flags(__ios_base::skipws);

  878. 

  879.       __is >> __x._M_b1 >> __x._M_b2;

  880. 

  881.       __is.flags(__flags);

  882.       return __is;

  883.     }

  884. 

  885. 

  886.   template<typename _IntType>

  887.     template<typename _UniformRandomNumberGenerator>

  888.       typename uniform_int<_IntType>::result_type

  889.       uniform_int<_IntType>::

  890.       _M_call(_UniformRandomNumberGenerator& __urng,

  891. 	      result_type __min, result_type __max, true_type)

  892.       {

  893. 	// XXX Must be fixed to work well for *arbitrary* __urng.max(),

  894. 	// __urng.min(), __max, __min.  Currently works fine only in the

  895. 	// most common case __urng.max() - __urng.min() >= __max - __min,

  896. 	// with __urng.max() > __urng.min() >= 0.

  897. 	typedef typename __gnu_cxx::__add_unsigned<typename

  898. 	  _UniformRandomNumberGenerator::result_type>::__type __urntype;

  899. 	typedef typename __gnu_cxx::__add_unsigned<result_type>::__type

  900. 	                                                      __utype;

  901. 	typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)

  902. 							> sizeof(__utype)),

  903. 	  __urntype, __utype>::__type                         __uctype;

  904. 

  905. 	result_type __ret;

  906. 

  907. 	const __urntype __urnmin = __urng.min();

  908. 	const __urntype __urnmax = __urng.max();

  909. 	const __urntype __urnrange = __urnmax - __urnmin;

  910. 	const __uctype __urange = __max - __min;

  911. 	const __uctype __udenom = (__urnrange <= __urange

  912. 				   ? 1 : __urnrange / (__urange + 1));

  913. 	do

  914. 	  __ret = (__urntype(__urng()) -  __urnmin) / __udenom;

  915. 	while (__ret > __max - __min);

  916. 

  917. 	return __ret + __min;

  918.       }

  919. 

  920.   template<typename _IntType, typename _CharT, typename _Traits>

  921.     std::basic_ostream<_CharT, _Traits>&

  922.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  923. 	       const uniform_int<_IntType>& __x)

  924.     {

  925.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  926.       typedef typename __ostream_type::ios_base    __ios_base;

  927. 

  928.       const typename __ios_base::fmtflags __flags = __os.flags();

  929.       const _CharT __fill = __os.fill();

  930.       const _CharT __space = __os.widen(' ');

  931.       __os.flags(__ios_base::scientific | __ios_base::left);

  932.       __os.fill(__space);

  933. 

  934.       __os << __x.min() << __space << __x.max();

  935. 

  936.       __os.flags(__flags);

  937.       __os.fill(__fill);

  938.       return __os;

  939.     }

  940. 

  941.   template<typename _IntType, typename _CharT, typename _Traits>

  942.     std::basic_istream<_CharT, _Traits>&

  943.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  944. 	       uniform_int<_IntType>& __x)

  945.     {

  946.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  947.       typedef typename __istream_type::ios_base    __ios_base;

  948. 

  949.       const typename __ios_base::fmtflags __flags = __is.flags();

  950.       __is.flags(__ios_base::dec | __ios_base::skipws);

  951. 

  952.       __is >> __x._M_min >> __x._M_max;

  953. 

  954.       __is.flags(__flags);

  955.       return __is;

  956.     }

  957. 

  958.   

  959.   template<typename _CharT, typename _Traits>

  960.     std::basic_ostream<_CharT, _Traits>&

  961.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  962. 	       const bernoulli_distribution& __x)

  963.     {

  964.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  965.       typedef typename __ostream_type::ios_base    __ios_base;

  966. 

  967.       const typename __ios_base::fmtflags __flags = __os.flags();

  968.       const _CharT __fill = __os.fill();

  969.       const std::streamsize __precision = __os.precision();

  970.       __os.flags(__ios_base::scientific | __ios_base::left);

  971.       __os.fill(__os.widen(' '));

  972.       __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);

  973. 

  974.       __os << __x.p();

  975. 

  976.       __os.flags(__flags);

  977.       __os.fill(__fill);

  978.       __os.precision(__precision);

  979.       return __os;

  980.     }

  981. 

  982. 

  983.   template<typename _IntType, typename _RealType>

  984.     template<class _UniformRandomNumberGenerator>

  985.       typename geometric_distribution<_IntType, _RealType>::result_type

  986.       geometric_distribution<_IntType, _RealType>::

  987.       operator()(_UniformRandomNumberGenerator& __urng)

  988.       {

  989. 	// About the epsilon thing see this thread:

  990.         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html

  991. 	const _RealType __naf =

  992. 	  (1 - std::numeric_limits<_RealType>::epsilon()) / 2;

  993. 	// The largest _RealType convertible to _IntType.

  994. 	const _RealType __thr =

  995. 	  std::numeric_limits<_IntType>::max() + __naf;

  996. 

  997. 	_RealType __cand;

  998. 	do

  999. 	  __cand = std::ceil(std::log(__urng()) / _M_log_p);

  1000. 	while (__cand >= __thr);

  1001. 

  1002. 	return result_type(__cand + __naf);

  1003.       }

  1004. 

  1005.   template<typename _IntType, typename _RealType,

  1006. 	   typename _CharT, typename _Traits>

  1007.     std::basic_ostream<_CharT, _Traits>&

  1008.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1009. 	       const geometric_distribution<_IntType, _RealType>& __x)

  1010.     {

  1011.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1012.       typedef typename __ostream_type::ios_base    __ios_base;

  1013. 

  1014.       const typename __ios_base::fmtflags __flags = __os.flags();

  1015.       const _CharT __fill = __os.fill();

  1016.       const std::streamsize __precision = __os.precision();

  1017.       __os.flags(__ios_base::scientific | __ios_base::left);

  1018.       __os.fill(__os.widen(' '));

  1019.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1020. 

  1021.       __os << __x.p();

  1022. 

  1023.       __os.flags(__flags);

  1024.       __os.fill(__fill);

  1025.       __os.precision(__precision);

  1026.       return __os;

  1027.     }

  1028. 

  1029. 

  1030.   template<typename _IntType, typename _RealType>

  1031.     void

  1032.     poisson_distribution<_IntType, _RealType>::

  1033.     _M_initialize()

  1034.     {

  1035. #if _GLIBCXX_USE_C99_MATH_TR1

  1036.       if (_M_mean >= 12)

  1037. 	{

  1038. 	  const _RealType __m = std::floor(_M_mean);

  1039. 	  _M_lm_thr = std::log(_M_mean);

  1040. 	  _M_lfm = std::tr1::lgamma(__m + 1);

  1041. 	  _M_sm = std::sqrt(__m);

  1042. 

  1043. 	  const _RealType __pi_4 = 0.7853981633974483096156608458198757L;

  1044. 	  const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m

  1045. 							      / __pi_4));

  1046. 	  _M_d = std::tr1::round(std::max(_RealType(6),

  1047. 					  std::min(__m, __dx)));

  1048. 	  const _RealType __cx = 2 * __m + _M_d;

  1049. 	  _M_scx = std::sqrt(__cx / 2);

  1050. 	  _M_1cx = 1 / __cx;

  1051. 

  1052. 	  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);

  1053. 	  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;

  1054. 	}

  1055.       else

  1056. #endif

  1057. 	_M_lm_thr = std::exp(-_M_mean);

  1058.       }

  1059. 

  1060.   /**

  1061.    * A rejection algorithm when mean >= 12 and a simple method based

  1062.    * upon the multiplication of uniform random variates otherwise.

  1063.    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1

  1064.    * is defined.

  1065.    *

  1066.    * Reference:

  1067.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,

  1068.    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).

  1069.    */

  1070.   template<typename _IntType, typename _RealType>

  1071.     template<class _UniformRandomNumberGenerator>

  1072.       typename poisson_distribution<_IntType, _RealType>::result_type

  1073.       poisson_distribution<_IntType, _RealType>::

  1074.       operator()(_UniformRandomNumberGenerator& __urng)

  1075.       {

  1076. #if _GLIBCXX_USE_C99_MATH_TR1

  1077. 	if (_M_mean >= 12)

  1078. 	  {

  1079. 	    _RealType __x;

  1080. 

  1081. 	    // See comments above...

  1082. 	    const _RealType __naf =

  1083. 	      (1 - std::numeric_limits<_RealType>::epsilon()) / 2;

  1084. 	    const _RealType __thr =

  1085. 	      std::numeric_limits<_IntType>::max() + __naf;

  1086. 

  1087. 	    const _RealType __m = std::floor(_M_mean);

  1088. 	    // sqrt(pi / 2)

  1089. 	    const _RealType __spi_2 = 1.2533141373155002512078826424055226L;

  1090. 	    const _RealType __c1 = _M_sm * __spi_2;

  1091. 	    const _RealType __c2 = _M_c2b + __c1; 

  1092. 	    const _RealType __c3 = __c2 + 1;

  1093. 	    const _RealType __c4 = __c3 + 1;

  1094. 	    // e^(1 / 78)

  1095. 	    const _RealType __e178 = 1.0129030479320018583185514777512983L;

  1096. 	    const _RealType __c5 = __c4 + __e178;

  1097. 	    const _RealType __c = _M_cb + __c5;

  1098. 	    const _RealType __2cx = 2 * (2 * __m + _M_d);

  1099. 

  1100. 	    bool __reject = true;

  1101. 	    do

  1102. 	      {

  1103. 		const _RealType __u = __c * __urng();

  1104. 		const _RealType __e = -std::log(__urng());

  1105. 

  1106. 		_RealType __w = 0.0;

  1107. 		

  1108. 		if (__u <= __c1)

  1109. 		  {

  1110. 		    const _RealType __n = _M_nd(__urng);

  1111. 		    const _RealType __y = -std::abs(__n) * _M_sm - 1;

  1112. 		    __x = std::floor(__y);

  1113. 		    __w = -__n * __n / 2;

  1114. 		    if (__x < -__m)

  1115. 		      continue;

  1116. 		  }

  1117. 		else if (__u <= __c2)

  1118. 		  {

  1119. 		    const _RealType __n = _M_nd(__urng);

  1120. 		    const _RealType __y = 1 + std::abs(__n) * _M_scx;

  1121. 		    __x = std::ceil(__y);

  1122. 		    __w = __y * (2 - __y) * _M_1cx;

  1123. 		    if (__x > _M_d)

  1124. 		      continue;

  1125. 		  }

  1126. 		else if (__u <= __c3)

  1127. 		  // NB: This case not in the book, nor in the Errata,

  1128. 		  // but should be ok...

  1129. 		  __x = -1;

  1130. 		else if (__u <= __c4)

  1131. 		  __x = 0;

  1132. 		else if (__u <= __c5)

  1133. 		  __x = 1;

  1134. 		else

  1135. 		  {

  1136. 		    const _RealType __v = -std::log(__urng());

  1137. 		    const _RealType __y = _M_d + __v * __2cx / _M_d;

  1138. 		    __x = std::ceil(__y);

  1139. 		    __w = -_M_d * _M_1cx * (1 + __y / 2);

  1140. 		  }

  1141. 

  1142. 		__reject = (__w - __e - __x * _M_lm_thr

  1143. 			    > _M_lfm - std::tr1::lgamma(__x + __m + 1));

  1144. 

  1145. 		__reject |= __x + __m >= __thr;

  1146. 

  1147. 	      } while (__reject);

  1148. 

  1149. 	    return result_type(__x + __m + __naf);

  1150. 	  }

  1151. 	else

  1152. #endif

  1153. 	  {

  1154. 	    _IntType     __x = 0;

  1155. 	    _RealType __prod = 1.0;

  1156. 

  1157. 	    do

  1158. 	      {

  1159. 		__prod *= __urng();

  1160. 		__x += 1;

  1161. 	      }

  1162. 	    while (__prod > _M_lm_thr);

  1163. 

  1164. 	    return __x - 1;

  1165. 	  }

  1166.       }

  1167. 

  1168.   template<typename _IntType, typename _RealType,

  1169. 	   typename _CharT, typename _Traits>

  1170.     std::basic_ostream<_CharT, _Traits>&

  1171.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1172. 	       const poisson_distribution<_IntType, _RealType>& __x)

  1173.     {

  1174.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1175.       typedef typename __ostream_type::ios_base    __ios_base;

  1176. 

  1177.       const typename __ios_base::fmtflags __flags = __os.flags();

  1178.       const _CharT __fill = __os.fill();

  1179.       const std::streamsize __precision = __os.precision();

  1180.       const _CharT __space = __os.widen(' ');

  1181.       __os.flags(__ios_base::scientific | __ios_base::left);

  1182.       __os.fill(__space);

  1183.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1184. 

  1185.       __os << __x.mean() << __space << __x._M_nd;

  1186. 

  1187.       __os.flags(__flags);

  1188.       __os.fill(__fill);

  1189.       __os.precision(__precision);

  1190.       return __os;

  1191.     }

  1192. 

  1193.   template<typename _IntType, typename _RealType,

  1194. 	   typename _CharT, typename _Traits>

  1195.     std::basic_istream<_CharT, _Traits>&

  1196.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  1197. 	       poisson_distribution<_IntType, _RealType>& __x)

  1198.     {

  1199.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  1200.       typedef typename __istream_type::ios_base    __ios_base;

  1201. 

  1202.       const typename __ios_base::fmtflags __flags = __is.flags();

  1203.       __is.flags(__ios_base::skipws);

  1204. 

  1205.       __is >> __x._M_mean >> __x._M_nd;

  1206.       __x._M_initialize();

  1207. 

  1208.       __is.flags(__flags);

  1209.       return __is;

  1210.     }

  1211. 

  1212. 

  1213.   template<typename _IntType, typename _RealType>

  1214.     void

  1215.     binomial_distribution<_IntType, _RealType>::

  1216.     _M_initialize()

  1217.     {

  1218.       const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;

  1219. 

  1220.       _M_easy = true;

  1221. 

  1222. #if _GLIBCXX_USE_C99_MATH_TR1

  1223.       if (_M_t * __p12 >= 8)

  1224. 	{

  1225. 	  _M_easy = false;

  1226. 	  const _RealType __np = std::floor(_M_t * __p12);

  1227. 	  const _RealType __pa = __np / _M_t;

  1228. 	  const _RealType __1p = 1 - __pa;

  1229. 	  

  1230. 	  const _RealType __pi_4 = 0.7853981633974483096156608458198757L;

  1231. 	  const _RealType __d1x =

  1232. 	    std::sqrt(__np * __1p * std::log(32 * __np

  1233. 					     / (81 * __pi_4 * __1p)));

  1234. 	  _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));

  1235. 	  const _RealType __d2x =

  1236. 	    std::sqrt(__np * __1p * std::log(32 * _M_t * __1p

  1237. 					     / (__pi_4 * __pa)));

  1238. 	  _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));

  1239. 	  

  1240. 	  // sqrt(pi / 2)

  1241. 	  const _RealType __spi_2 = 1.2533141373155002512078826424055226L;

  1242. 	  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));

  1243. 	  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));

  1244. 	  _M_c = 2 * _M_d1 / __np;

  1245. 	  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;

  1246. 	  const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;

  1247. 	  const _RealType __s1s = _M_s1 * _M_s1;

  1248. 	  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))

  1249. 			     * 2 * __s1s / _M_d1

  1250. 			     * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));

  1251. 	  const _RealType __s2s = _M_s2 * _M_s2;

  1252. 	  _M_s = (_M_a123 + 2 * __s2s / _M_d2

  1253. 		  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));

  1254. 	  _M_lf = (std::tr1::lgamma(__np + 1)

  1255. 		   + std::tr1::lgamma(_M_t - __np + 1));

  1256. 	  _M_lp1p = std::log(__pa / __1p);

  1257. 

  1258. 	  _M_q = -std::log(1 - (__p12 - __pa) / __1p);

  1259. 	}

  1260.       else

  1261. #endif

  1262. 	_M_q = -std::log(1 - __p12);

  1263.     }

  1264. 

  1265.   template<typename _IntType, typename _RealType>

  1266.     template<class _UniformRandomNumberGenerator>

  1267.       typename binomial_distribution<_IntType, _RealType>::result_type

  1268.       binomial_distribution<_IntType, _RealType>::

  1269.       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)

  1270.       {

  1271. 	_IntType    __x = 0;

  1272. 	_RealType __sum = 0;

  1273. 

  1274. 	do

  1275. 	  {

  1276. 	    const _RealType __e = -std::log(__urng());

  1277. 	    __sum += __e / (__t - __x);

  1278. 	    __x += 1;

  1279. 	  }

  1280. 	while (__sum <= _M_q);

  1281. 

  1282. 	return __x - 1;

  1283.       }

  1284. 

  1285.   /**

  1286.    * A rejection algorithm when t * p >= 8 and a simple waiting time

  1287.    * method - the second in the referenced book - otherwise.

  1288.    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1

  1289.    * is defined.

  1290.    *

  1291.    * Reference:

  1292.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,

  1293.    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).

  1294.    */

  1295.   template<typename _IntType, typename _RealType>

  1296.     template<class _UniformRandomNumberGenerator>

  1297.       typename binomial_distribution<_IntType, _RealType>::result_type

  1298.       binomial_distribution<_IntType, _RealType>::

  1299.       operator()(_UniformRandomNumberGenerator& __urng)

  1300.       {

  1301. 	result_type __ret;

  1302. 	const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;

  1303. 

  1304. #if _GLIBCXX_USE_C99_MATH_TR1

  1305. 	if (!_M_easy)

  1306. 	  {

  1307. 	    _RealType __x;

  1308. 

  1309. 	    // See comments above...

  1310. 	    const _RealType __naf =

  1311. 	      (1 - std::numeric_limits<_RealType>::epsilon()) / 2;

  1312. 	    const _RealType __thr =

  1313. 	      std::numeric_limits<_IntType>::max() + __naf;

  1314. 

  1315. 	    const _RealType __np = std::floor(_M_t * __p12);

  1316. 	    const _RealType __pa = __np / _M_t;

  1317. 

  1318. 	    // sqrt(pi / 2)

  1319. 	    const _RealType __spi_2 = 1.2533141373155002512078826424055226L;

  1320. 	    const _RealType __a1 = _M_a1;

  1321. 	    const _RealType __a12 = __a1 + _M_s2 * __spi_2;

  1322. 	    const _RealType __a123 = _M_a123;

  1323. 	    const _RealType __s1s = _M_s1 * _M_s1;

  1324. 	    const _RealType __s2s = _M_s2 * _M_s2;

  1325. 

  1326. 	    bool __reject;

  1327. 	    do

  1328. 	      {

  1329. 		const _RealType __u = _M_s * __urng();

  1330. 

  1331. 		_RealType __v;

  1332. 

  1333. 		if (__u <= __a1)

  1334. 		  {

  1335. 		    const _RealType __n = _M_nd(__urng);

  1336. 		    const _RealType __y = _M_s1 * std::abs(__n);

  1337. 		    __reject = __y >= _M_d1;

  1338. 		    if (!__reject)

  1339. 		      {

  1340. 			const _RealType __e = -std::log(__urng());

  1341. 			__x = std::floor(__y);

  1342. 			__v = -__e - __n * __n / 2 + _M_c;

  1343. 		      }

  1344. 		  }

  1345. 		else if (__u <= __a12)

  1346. 		  {

  1347. 		    const _RealType __n = _M_nd(__urng);

  1348. 		    const _RealType __y = _M_s2 * std::abs(__n);

  1349. 		    __reject = __y >= _M_d2;

  1350. 		    if (!__reject)

  1351. 		      {

  1352. 			const _RealType __e = -std::log(__urng());

  1353. 			__x = std::floor(-__y);

  1354. 			__v = -__e - __n * __n / 2;

  1355. 		      }

  1356. 		  }

  1357. 		else if (__u <= __a123)

  1358. 		  {

  1359. 		    const _RealType __e1 = -std::log(__urng());		    

  1360. 		    const _RealType __e2 = -std::log(__urng());

  1361. 

  1362. 		    const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;

  1363. 		    __x = std::floor(__y);

  1364. 		    __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)

  1365. 					    -__y / (2 * __s1s)));

  1366. 		    __reject = false;

  1367. 		  }

  1368. 		else

  1369. 		  {

  1370. 		    const _RealType __e1 = -std::log(__urng());		    

  1371. 		    const _RealType __e2 = -std::log(__urng());

  1372. 

  1373. 		    const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;

  1374. 		    __x = std::floor(-__y);

  1375. 		    __v = -__e2 - _M_d2 * __y / (2 * __s2s);

  1376. 		    __reject = false;

  1377. 		  }

  1378. 

  1379. 		__reject = __reject || __x < -__np || __x > _M_t - __np;

  1380. 		if (!__reject)

  1381. 		  {

  1382. 		    const _RealType __lfx =

  1383. 		      std::tr1::lgamma(__np + __x + 1)

  1384. 		      + std::tr1::lgamma(_M_t - (__np + __x) + 1);

  1385. 		    __reject = __v > _M_lf - __lfx + __x * _M_lp1p;

  1386. 		  }

  1387. 

  1388. 		__reject |= __x + __np >= __thr;

  1389. 	      }

  1390. 	    while (__reject);

  1391. 

  1392. 	    __x += __np + __naf;

  1393. 

  1394. 	    const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 

  1395. 	    __ret = _IntType(__x) + __z;

  1396. 	  }

  1397. 	else

  1398. #endif

  1399. 	  __ret = _M_waiting(__urng, _M_t);

  1400. 

  1401. 	if (__p12 != _M_p)

  1402. 	  __ret = _M_t - __ret;

  1403. 	return __ret;

  1404.       }

  1405. 

  1406.   template<typename _IntType, typename _RealType,

  1407. 	   typename _CharT, typename _Traits>

  1408.     std::basic_ostream<_CharT, _Traits>&

  1409.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1410. 	       const binomial_distribution<_IntType, _RealType>& __x)

  1411.     {

  1412.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1413.       typedef typename __ostream_type::ios_base    __ios_base;

  1414. 

  1415.       const typename __ios_base::fmtflags __flags = __os.flags();

  1416.       const _CharT __fill = __os.fill();

  1417.       const std::streamsize __precision = __os.precision();

  1418.       const _CharT __space = __os.widen(' ');

  1419.       __os.flags(__ios_base::scientific | __ios_base::left);

  1420.       __os.fill(__space);

  1421.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1422. 

  1423.       __os << __x.t() << __space << __x.p() 

  1424. 	   << __space << __x._M_nd;

  1425. 

  1426.       __os.flags(__flags);

  1427.       __os.fill(__fill);

  1428.       __os.precision(__precision);

  1429.       return __os;

  1430.     }

  1431. 

  1432.   template<typename _IntType, typename _RealType,

  1433. 	   typename _CharT, typename _Traits>

  1434.     std::basic_istream<_CharT, _Traits>&

  1435.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  1436. 	       binomial_distribution<_IntType, _RealType>& __x)

  1437.     {

  1438.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  1439.       typedef typename __istream_type::ios_base    __ios_base;

  1440. 

  1441.       const typename __ios_base::fmtflags __flags = __is.flags();

  1442.       __is.flags(__ios_base::dec | __ios_base::skipws);

  1443. 

  1444.       __is >> __x._M_t >> __x._M_p >> __x._M_nd;

  1445.       __x._M_initialize();

  1446. 

  1447.       __is.flags(__flags);

  1448.       return __is;

  1449.     }

  1450. 

  1451. 

  1452.   template<typename _RealType, typename _CharT, typename _Traits>

  1453.     std::basic_ostream<_CharT, _Traits>&

  1454.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1455. 	       const uniform_real<_RealType>& __x)

  1456.     {

  1457.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1458.       typedef typename __ostream_type::ios_base    __ios_base;

  1459. 

  1460.       const typename __ios_base::fmtflags __flags = __os.flags();

  1461.       const _CharT __fill = __os.fill();

  1462.       const std::streamsize __precision = __os.precision();

  1463.       const _CharT __space = __os.widen(' ');

  1464.       __os.flags(__ios_base::scientific | __ios_base::left);

  1465.       __os.fill(__space);

  1466.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1467. 

  1468.       __os << __x.min() << __space << __x.max();

  1469. 

  1470.       __os.flags(__flags);

  1471.       __os.fill(__fill);

  1472.       __os.precision(__precision);

  1473.       return __os;

  1474.     }

  1475. 

  1476.   template<typename _RealType, typename _CharT, typename _Traits>

  1477.     std::basic_istream<_CharT, _Traits>&

  1478.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  1479. 	       uniform_real<_RealType>& __x)

  1480.     {

  1481.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  1482.       typedef typename __istream_type::ios_base    __ios_base;

  1483. 

  1484.       const typename __ios_base::fmtflags __flags = __is.flags();

  1485.       __is.flags(__ios_base::skipws);

  1486. 

  1487.       __is >> __x._M_min >> __x._M_max;

  1488. 

  1489.       __is.flags(__flags);

  1490.       return __is;

  1491.     }

  1492. 

  1493. 

  1494.   template<typename _RealType, typename _CharT, typename _Traits>

  1495.     std::basic_ostream<_CharT, _Traits>&

  1496.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1497. 	       const exponential_distribution<_RealType>& __x)

  1498.     {

  1499.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1500.       typedef typename __ostream_type::ios_base    __ios_base;

  1501. 

  1502.       const typename __ios_base::fmtflags __flags = __os.flags();

  1503.       const _CharT __fill = __os.fill();

  1504.       const std::streamsize __precision = __os.precision();

  1505.       __os.flags(__ios_base::scientific | __ios_base::left);

  1506.       __os.fill(__os.widen(' '));

  1507.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1508. 

  1509.       __os << __x.lambda();

  1510. 

  1511.       __os.flags(__flags);

  1512.       __os.fill(__fill);

  1513.       __os.precision(__precision);

  1514.       return __os;

  1515.     }

  1516. 

  1517. 

  1518.   /**

  1519.    * Polar method due to Marsaglia.

  1520.    *

  1521.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,

  1522.    * New York, 1986, Ch. V, Sect. 4.4.

  1523.    */

  1524.   template<typename _RealType>

  1525.     template<class _UniformRandomNumberGenerator>

  1526.       typename normal_distribution<_RealType>::result_type

  1527.       normal_distribution<_RealType>::

  1528.       operator()(_UniformRandomNumberGenerator& __urng)

  1529.       {

  1530. 	result_type __ret;

  1531. 

  1532. 	if (_M_saved_available)

  1533. 	  {

  1534. 	    _M_saved_available = false;

  1535. 	    __ret = _M_saved;

  1536. 	  }

  1537. 	else

  1538. 	  {

  1539. 	    result_type __x, __y, __r2;

  1540. 	    do

  1541. 	      {

  1542. 		__x = result_type(2.0) * __urng() - 1.0;

  1543. 		__y = result_type(2.0) * __urng() - 1.0;

  1544. 		__r2 = __x * __x + __y * __y;

  1545. 	      }

  1546. 	    while (__r2 > 1.0 || __r2 == 0.0);

  1547. 

  1548. 	    const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);

  1549. 	    _M_saved = __x * __mult;

  1550. 	    _M_saved_available = true;

  1551. 	    __ret = __y * __mult;

  1552. 	  }

  1553. 	

  1554. 	__ret = __ret * _M_sigma + _M_mean;

  1555. 	return __ret;

  1556.       }

  1557. 

  1558.   template<typename _RealType, typename _CharT, typename _Traits>

  1559.     std::basic_ostream<_CharT, _Traits>&

  1560.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1561. 	       const normal_distribution<_RealType>& __x)

  1562.     {

  1563.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1564.       typedef typename __ostream_type::ios_base    __ios_base;

  1565. 

  1566.       const typename __ios_base::fmtflags __flags = __os.flags();

  1567.       const _CharT __fill = __os.fill();

  1568.       const std::streamsize __precision = __os.precision();

  1569.       const _CharT __space = __os.widen(' ');

  1570.       __os.flags(__ios_base::scientific | __ios_base::left);

  1571.       __os.fill(__space);

  1572.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1573. 

  1574.       __os << __x._M_saved_available << __space

  1575. 	   << __x.mean() << __space

  1576. 	   << __x.sigma();

  1577.       if (__x._M_saved_available)

  1578. 	__os << __space << __x._M_saved;

  1579. 

  1580.       __os.flags(__flags);

  1581.       __os.fill(__fill);

  1582.       __os.precision(__precision);

  1583.       return __os;

  1584.     }

  1585. 

  1586.   template<typename _RealType, typename _CharT, typename _Traits>

  1587.     std::basic_istream<_CharT, _Traits>&

  1588.     operator>>(std::basic_istream<_CharT, _Traits>& __is,

  1589. 	       normal_distribution<_RealType>& __x)

  1590.     {

  1591.       typedef std::basic_istream<_CharT, _Traits>  __istream_type;

  1592.       typedef typename __istream_type::ios_base    __ios_base;

  1593. 

  1594.       const typename __ios_base::fmtflags __flags = __is.flags();

  1595.       __is.flags(__ios_base::dec | __ios_base::skipws);

  1596. 

  1597.       __is >> __x._M_saved_available >> __x._M_mean

  1598. 	   >> __x._M_sigma;

  1599.       if (__x._M_saved_available)

  1600. 	__is >> __x._M_saved;

  1601. 

  1602.       __is.flags(__flags);

  1603.       return __is;

  1604.     }

  1605. 

  1606. 

  1607.   template<typename _RealType>

  1608.     void

  1609.     gamma_distribution<_RealType>::

  1610.     _M_initialize()

  1611.     {

  1612.       if (_M_alpha >= 1)

  1613. 	_M_l_d = std::sqrt(2 * _M_alpha - 1);

  1614.       else

  1615. 	_M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))

  1616. 		  * (1 - _M_alpha));

  1617.     }

  1618. 

  1619.   /**

  1620.    * Cheng's rejection algorithm GB for alpha >= 1 and a modification

  1621.    * of Vaduva's rejection from Weibull algorithm due to Devroye for

  1622.    * alpha < 1.

  1623.    *

  1624.    * References:

  1625.    * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral

  1626.    * Shape Parameter. Applied Statistics, 26, 71-75, 1977.

  1627.    *

  1628.    * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection

  1629.    * and Composition Procedures. Math. Operationsforschung and Statistik,

  1630.    * Series in Statistics, 8, 545-576, 1977.

  1631.    *

  1632.    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,

  1633.    * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).

  1634.    */

  1635.   template<typename _RealType>

  1636.     template<class _UniformRandomNumberGenerator>

  1637.       typename gamma_distribution<_RealType>::result_type

  1638.       gamma_distribution<_RealType>::

  1639.       operator()(_UniformRandomNumberGenerator& __urng)

  1640.       {

  1641. 	result_type __x;

  1642. 

  1643. 	bool __reject;

  1644. 	if (_M_alpha >= 1)

  1645. 	  {

  1646. 	    // alpha - log(4)

  1647. 	    const result_type __b = _M_alpha

  1648. 	      - result_type(1.3862943611198906188344642429163531L);

  1649. 	    const result_type __c = _M_alpha + _M_l_d;

  1650. 	    const result_type __1l = 1 / _M_l_d;

  1651. 

  1652. 	    // 1 + log(9 / 2)

  1653. 	    const result_type __k = 2.5040773967762740733732583523868748L;

  1654. 

  1655. 	    do

  1656. 	      {

  1657. 		const result_type __u = __urng();

  1658. 		const result_type __v = __urng();

  1659. 

  1660. 		const result_type __y = __1l * std::log(__v / (1 - __v));

  1661. 		__x = _M_alpha * std::exp(__y);

  1662. 

  1663. 		const result_type __z = __u * __v * __v;

  1664. 		const result_type __r = __b + __c * __y - __x;

  1665. 

  1666. 		__reject = __r < result_type(4.5) * __z - __k;

  1667. 		if (__reject)

  1668. 		  __reject = __r < std::log(__z);

  1669. 	      }

  1670. 	    while (__reject);

  1671. 	  }

  1672. 	else

  1673. 	  {

  1674. 	    const result_type __c = 1 / _M_alpha;

  1675. 

  1676. 	    do

  1677. 	      {

  1678. 		const result_type __z = -std::log(__urng());

  1679. 		const result_type __e = -std::log(__urng());

  1680. 

  1681. 		__x = std::pow(__z, __c);

  1682. 

  1683. 		__reject = __z + __e < _M_l_d + __x;

  1684. 	      }

  1685. 	    while (__reject);

  1686. 	  }

  1687. 

  1688. 	return __x;

  1689.       }

  1690. 

  1691.   template<typename _RealType, typename _CharT, typename _Traits>

  1692.     std::basic_ostream<_CharT, _Traits>&

  1693.     operator<<(std::basic_ostream<_CharT, _Traits>& __os,

  1694. 	       const gamma_distribution<_RealType>& __x)

  1695.     {

  1696.       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;

  1697.       typedef typename __ostream_type::ios_base    __ios_base;

  1698. 

  1699.       const typename __ios_base::fmtflags __flags = __os.flags();

  1700.       const _CharT __fill = __os.fill();

  1701.       const std::streamsize __precision = __os.precision();

  1702.       __os.flags(__ios_base::scientific | __ios_base::left);

  1703.       __os.fill(__os.widen(' '));

  1704.       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

  1705. 

  1706.       __os << __x.alpha();

  1707. 

  1708.       __os.flags(__flags);

  1709.       __os.fill(__fill);

  1710.       __os.precision(__precision);

  1711.       return __os;

  1712.     }

  1713. }

  1714. 

  1715. _GLIBCXX_END_NAMESPACE_VERSION

  1716. }

  1717. 

  1718. #endif