FM / tr1

/usr/include/c++/8/tr1

Name	Size	Mode	Actions
array	6964	0644	edit dl rm
bessel_function.tcc	22473	0644	edit dl rm
beta_function.tcc	5995	0644	edit dl rm
ccomplex	1255	0644	edit dl rm
cctype	1412	0644	edit dl rm
cfenv	2004	0644	edit dl rm
cfloat	1380	0644	edit dl rm
cinttypes	2256	0644	edit dl rm
climits	1454	0644	edit dl rm
cmath	43806	0644	edit dl rm
complex	12384	0644	edit dl rm
complex.h	1261	0644	edit dl rm
cstdarg	1246	0644	edit dl rm
cstdbool	1344	0644	edit dl rm
cstdint	2623	0644	edit dl rm
cstdio	1482	0644	edit dl rm
cstdlib	1796	0644	edit dl rm
ctgmath	1248	0644	edit dl rm
ctime	1234	0644	edit dl rm
ctype.h	1209	0644	edit dl rm
cwchar	1718	0644	edit dl rm
cwctype	1459	0644	edit dl rm
ell_integral.tcc	27724	0644	edit dl rm
exp_integral.tcc	16009	0644	edit dl rm
fenv.h	1204	0644	edit dl rm
float.h	1209	0644	edit dl rm
functional	70545	0644	edit dl rm
functional_hash.h	6043	0644	edit dl rm
gamma.tcc	14682	0644	edit dl rm
hashtable.h	41537	0644	edit dl rm
hashtable_policy.h	25086	0644	edit dl rm
hypergeometric.tcc	28066	0644	edit dl rm
inttypes.h	1267	0644	edit dl rm
legendre_function.tcc	10909	0644	edit dl rm
limits.h	1214	0644	edit dl rm
math.h	4553	0644	edit dl rm
memory	1791	0644	edit dl rm
modified_bessel_func.tcc	16320	0644	edit dl rm
poly_hermite.tcc	3925	0644	edit dl rm
poly_laguerre.tcc	11676	0644	edit dl rm
random	1589	0644	edit dl rm
random.h	73123	0644	edit dl rm
random.tcc	53927	0644	edit dl rm
regex	92880	0644	edit dl rm
riemann_zeta.tcc	14063	0644	edit dl rm
shared_ptr.h	32608	0644	edit dl rm
special_function_util.h	5055	0644	edit dl rm
stdarg.h	1214	0644	edit dl rm
stdbool.h	1219	0644	edit dl rm
stdint.h	1214	0644	edit dl rm
stdio.h	1209	0644	edit dl rm
stdlib.h	1487	0644	edit dl rm
tgmath.h	1255	0644	edit dl rm
tuple	12119	0644	edit dl rm
type_traits	19019	0644	edit dl rm
unordered_map	1574	0644	edit dl rm
unordered_map.h	10216	0644	edit dl rm
unordered_set	1574	0644	edit dl rm
unordered_set.h	9540	0644	edit dl rm
utility	3225	0644	edit dl rm
wchar.h	1249	0644	edit dl rm
wctype.h	1255	0644	edit dl rm

Edit: /usr/include/c++/8/tr1/random.tcc (53927B)

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

// Copyright (C) 2009-2018 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library.  This library is free
// software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 3, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// Under Section 7 of GPL version 3, you are granted additional
// permissions described in the GCC Runtime Library Exception, version
// 3.1, as published by the Free Software Foundation.

// You should have received a copy of the GNU General Public License and
// a copy of the GCC Runtime Library Exception along with this program;
// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
// <http://www.gnu.org/licenses/>.

/** @file tr1/random.tcc
 *  This is an internal header file, included by other library headers.
 *  Do not attempt to use it directly. @headername{tr1/random}
 */

#ifndef _GLIBCXX_TR1_RANDOM_TCC
#define _GLIBCXX_TR1_RANDOM_TCC 1

namespace std _GLIBCXX_VISIBILITY(default)
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION

namespace tr1
{
  /*
   * (Further) implementation-space details.
   */
  namespace __detail
  {
    // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
    // integer overflow.
    //
    // Because a and c are compile-time integral constants the compiler kindly
    // elides any unreachable paths.
    //
    // Preconditions:  a > 0, m > 0.
    //
    template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
      struct _Mod
      {
	static _Tp
	__calc(_Tp __x)
	{
	  if (__a == 1)
	    __x %= __m;
	  else
	    {
	      static const _Tp __q = __m / __a;
	      static const _Tp __r = __m % __a;
	      
	      _Tp __t1 = __a * (__x % __q);
	      _Tp __t2 = __r * (__x / __q);
	      if (__t1 >= __t2)
		__x = __t1 - __t2;
	      else
		__x = __m - __t2 + __t1;
	    }

if (__c != 0)
	    {
	      const _Tp __d = __m - __x;
	      if (__d > __c)
		__x += __c;
	      else
		__x = __c - __d;
	    }
	  return __x;
	}
      };

// Special case for m == 0 -- use unsigned integer overflow as modulo
    // operator.
    template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
      struct _Mod<_Tp, __a, __c, __m, true>
      {
	static _Tp
	__calc(_Tp __x)
	{ return __a * __x + __c; }
      };
  } // namespace __detail

template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    const _UIntType
    linear_congruential<_UIntType, __a, __c, __m>::multiplier;

template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    const _UIntType
    linear_congruential<_UIntType, __a, __c, __m>::increment;

template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    const _UIntType
    linear_congruential<_UIntType, __a, __c, __m>::modulus;

/**
   * Seeds the LCR with integral value @p __x0, adjusted so that the 
   * ring identity is never a member of the convergence set.
   */
  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    void
    linear_congruential<_UIntType, __a, __c, __m>::
    seed(unsigned long __x0)
    {
      if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
	  && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
	_M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
      else
	_M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
    }

/**
   * Seeds the LCR engine with a value generated by @p __g.
   */
  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    template<class _Gen>
      void
      linear_congruential<_UIntType, __a, __c, __m>::
      seed(_Gen& __g, false_type)
      {
	_UIntType __x0 = __g();
	if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
	    && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
	  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
	else
	  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
      }

/**
   * Gets the next generated value in sequence.
   */
  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
    typename linear_congruential<_UIntType, __a, __c, __m>::result_type
    linear_congruential<_UIntType, __a, __c, __m>::
    operator()()
    {
      _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
      return _M_x;
    }

template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __os.flags();
      const _CharT __fill = __os.fill();
      __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
      __os.fill(__os.widen(' '));

__os << __lcr._M_x;

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       linear_congruential<_UIntType, __a, __c, __m>& __lcr)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

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

__is >> __lcr._M_x;

__is.flags(__flags);
      return __is;
    }

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    const int
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::state_size;
    
  template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    const int
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::shift_size;

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    const _UIntType
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::parameter_a;

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    void
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::
    seed(unsigned long __value)
    {
      _M_x[0] = __detail::__mod<_UIntType, 1, 0,
	__detail::_Shift<_UIntType, __w>::__value>(__value);

for (int __i = 1; __i < state_size; ++__i)
	{
	  _UIntType __x = _M_x[__i - 1];
	  __x ^= __x >> (__w - 2);
	  __x *= 1812433253ul;
	  __x += __i;
	  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
	    __detail::_Shift<_UIntType, __w>::__value>(__x);	  
	}
      _M_p = state_size;
    }

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    template<class _Gen>
      void
      mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		       __b, __t, __c, __l>::
      seed(_Gen& __gen, false_type)
      {
	for (int __i = 0; __i < state_size; ++__i)
	  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
	    __detail::_Shift<_UIntType, __w>::__value>(__gen());
	_M_p = state_size;
      }

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s,
	   _UIntType __b, int __t, _UIntType __c, int __l>
    typename
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::result_type
    mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
		     __b, __t, __c, __l>::
    operator()()
    {
      // Reload the vector - cost is O(n) amortized over n calls.
      if (_M_p >= state_size)
	{
	  const _UIntType __upper_mask = (~_UIntType()) << __r;
	  const _UIntType __lower_mask = ~__upper_mask;

for (int __k = 0; __k < (__n - __m); ++__k)
	    {
	      _UIntType __y = ((_M_x[__k] & __upper_mask)
			       | (_M_x[__k + 1] & __lower_mask));
	      _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
			   ^ ((__y & 0x01) ? __a : 0));
	    }

for (int __k = (__n - __m); __k < (__n - 1); ++__k)
	    {
	      _UIntType __y = ((_M_x[__k] & __upper_mask)
			       | (_M_x[__k + 1] & __lower_mask));
	      _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
			   ^ ((__y & 0x01) ? __a : 0));
	    }

_UIntType __y = ((_M_x[__n - 1] & __upper_mask)
			   | (_M_x[0] & __lower_mask));
	  _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
			   ^ ((__y & 0x01) ? __a : 0));
	  _M_p = 0;
	}

// Calculate o(x(i)).
      result_type __z = _M_x[_M_p++];
      __z ^= (__z >> __u);
      __z ^= (__z << __s) & __b;
      __z ^= (__z << __t) & __c;
      __z ^= (__z >> __l);

return __z;
    }

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s, _UIntType __b, int __t,
	   _UIntType __c, int __l,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const mersenne_twister<_UIntType, __w, __n, __m,
	       __r, __a, __u, __s, __b, __t, __c, __l>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __os.flags();
      const _CharT __fill = __os.fill();
      const _CharT __space = __os.widen(' ');
      __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
      __os.fill(__space);

for (int __i = 0; __i < __n - 1; ++__i)
	__os << __x._M_x[__i] << __space;
      __os << __x._M_x[__n - 1];

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<class _UIntType, int __w, int __n, int __m, int __r,
	   _UIntType __a, int __u, int __s, _UIntType __b, int __t,
	   _UIntType __c, int __l,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       mersenne_twister<_UIntType, __w, __n, __m,
	       __r, __a, __u, __s, __b, __t, __c, __l>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

for (int __i = 0; __i < __n; ++__i)
	__is >> __x._M_x[__i];

__is.flags(__flags);
      return __is;
    }

template<typename _IntType, _IntType __m, int __s, int __r>
    const _IntType
    subtract_with_carry<_IntType, __m, __s, __r>::modulus;

template<typename _IntType, _IntType __m, int __s, int __r>
    const int
    subtract_with_carry<_IntType, __m, __s, __r>::long_lag;

template<typename _IntType, _IntType __m, int __s, int __r>
    const int
    subtract_with_carry<_IntType, __m, __s, __r>::short_lag;

template<typename _IntType, _IntType __m, int __s, int __r>
    void
    subtract_with_carry<_IntType, __m, __s, __r>::
    seed(unsigned long __value)
    {
      if (__value == 0)
	__value = 19780503;

std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
	__lcg(__value);

for (int __i = 0; __i < long_lag; ++__i)
	_M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());

_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
      _M_p = 0;
    }

template<typename _IntType, _IntType __m, int __s, int __r>
    template<class _Gen>
      void
      subtract_with_carry<_IntType, __m, __s, __r>::
      seed(_Gen& __gen, false_type)
      {
	const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;

for (int __i = 0; __i < long_lag; ++__i)
	  {
	    _UIntType __tmp = 0;
	    _UIntType __factor = 1;
	    for (int __j = 0; __j < __n; ++__j)
	      {
		__tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
		         (__gen()) * __factor;
		__factor *= __detail::_Shift<_UIntType, 32>::__value;
	      }
	    _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
	  }
	_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
	_M_p = 0;
      }

template<typename _IntType, _IntType __m, int __s, int __r>
    typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
    subtract_with_carry<_IntType, __m, __s, __r>::
    operator()()
    {
      // Derive short lag index from current index.
      int __ps = _M_p - short_lag;
      if (__ps < 0)
	__ps += long_lag;

// Calculate new x(i) without overflow or division.
      // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
      // cannot overflow.
      _UIntType __xi;
      if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
	{
	  __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
	  _M_carry = 0;
	}
      else
	{
	  __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
	  _M_carry = 1;
	}
      _M_x[_M_p] = __xi;

// Adjust current index to loop around in ring buffer.
      if (++_M_p >= long_lag)
	_M_p = 0;

return __xi;
    }

template<typename _IntType, _IntType __m, int __s, int __r,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const subtract_with_carry<_IntType, __m, __s, __r>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

for (int __i = 0; __i < __r; ++__i)
	__os << __x._M_x[__i] << __space;
      __os << __x._M_carry;

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<typename _IntType, _IntType __m, int __s, int __r,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       subtract_with_carry<_IntType, __m, __s, __r>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

for (int __i = 0; __i < __r; ++__i)
	__is >> __x._M_x[__i];
      __is >> __x._M_carry;

__is.flags(__flags);
      return __is;
    }

template<typename _RealType, int __w, int __s, int __r>
    const int
    subtract_with_carry_01<_RealType, __w, __s, __r>::word_size;

template<typename _RealType, int __w, int __s, int __r>
    const int
    subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag;

template<typename _RealType, int __w, int __s, int __r>
    const int
    subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag;

template<typename _RealType, int __w, int __s, int __r>
    void
    subtract_with_carry_01<_RealType, __w, __s, __r>::
    _M_initialize_npows()
    {
      for (int __j = 0; __j < __n; ++__j)
#if _GLIBCXX_USE_C99_MATH_TR1
	_M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
#else
        _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
#endif
    }

template<typename _RealType, int __w, int __s, int __r>
    void
    subtract_with_carry_01<_RealType, __w, __s, __r>::
    seed(unsigned long __value)
    {
      if (__value == 0)
	__value = 19780503;

// _GLIBCXX_RESOLVE_LIB_DEFECTS
      // 512. Seeding subtract_with_carry_01 from a single unsigned long.
      std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
	__lcg(__value);

this->seed(__lcg);
    }

template<typename _RealType, int __w, int __s, int __r>
    template<class _Gen>
      void
      subtract_with_carry_01<_RealType, __w, __s, __r>::
      seed(_Gen& __gen, false_type)
      {
	for (int __i = 0; __i < long_lag; ++__i)
	  {
	    for (int __j = 0; __j < __n - 1; ++__j)
	      _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
	    _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
	      __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
	  }

_M_carry = 1;
	for (int __j = 0; __j < __n; ++__j)
	  if (_M_x[long_lag - 1][__j] != 0)
	    {
	      _M_carry = 0;
	      break;
	    }

_M_p = 0;
      }

template<typename _RealType, int __w, int __s, int __r>
    typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
    subtract_with_carry_01<_RealType, __w, __s, __r>::
    operator()()
    {
      // Derive short lag index from current index.
      int __ps = _M_p - short_lag;
      if (__ps < 0)
	__ps += long_lag;

_UInt32Type __new_carry;
      for (int __j = 0; __j < __n - 1; ++__j)
	{
	  if (_M_x[__ps][__j] > _M_x[_M_p][__j]
	      || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
	    __new_carry = 0;
	  else
	    __new_carry = 1;

_M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
	  _M_carry = __new_carry;
	}

if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
	  || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
	__new_carry = 0;
      else
	__new_carry = 1;
      
      _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
	__detail::_Shift<_UInt32Type, __w % 32>::__value>
	(_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
      _M_carry = __new_carry;

result_type __ret = 0.0;
      for (int __j = 0; __j < __n; ++__j)
	__ret += _M_x[_M_p][__j] * _M_npows[__j];

// Adjust current index to loop around in ring buffer.
      if (++_M_p >= long_lag)
	_M_p = 0;

return __ret;
    }

template<typename _RealType, int __w, int __s, int __r,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

for (int __i = 0; __i < __r; ++__i)
	for (int __j = 0; __j < __x.__n; ++__j)
	  __os << __x._M_x[__i][__j] << __space;
      __os << __x._M_carry;

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<typename _RealType, int __w, int __s, int __r,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

for (int __i = 0; __i < __r; ++__i)
	for (int __j = 0; __j < __x.__n; ++__j)
	  __is >> __x._M_x[__i][__j];
      __is >> __x._M_carry;

__is.flags(__flags);
      return __is;
    }

template<class _UniformRandomNumberGenerator, int __p, int __r>
    const int
    discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size;

template<class _UniformRandomNumberGenerator, int __p, int __r>
    const int
    discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block;

template<class _UniformRandomNumberGenerator, int __p, int __r>
    typename discard_block<_UniformRandomNumberGenerator,
			   __p, __r>::result_type
    discard_block<_UniformRandomNumberGenerator, __p, __r>::
    operator()()
    {
      if (_M_n >= used_block)
	{
	  while (_M_n < block_size)
	    {
	      _M_b();
	      ++_M_n;
	    }
	  _M_n = 0;
	}
      ++_M_n;
      return _M_b();
    }

template<class _UniformRandomNumberGenerator, int __p, int __r,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const discard_block<_UniformRandomNumberGenerator,
	       __p, __r>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

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

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<class _UniformRandomNumberGenerator, int __p, int __r,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

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

__is.flags(__flags);
      return __is;
    }

template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2>
    const int
    xor_combine<_UniformRandomNumberGenerator1, __s1,
		_UniformRandomNumberGenerator2, __s2>::shift1;
     
  template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2>
    const int
    xor_combine<_UniformRandomNumberGenerator1, __s1,
		_UniformRandomNumberGenerator2, __s2>::shift2;

template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2>
    void
    xor_combine<_UniformRandomNumberGenerator1, __s1,
		_UniformRandomNumberGenerator2, __s2>::
    _M_initialize_max()
    {
      const int __w = std::numeric_limits<result_type>::digits;

const result_type __m1 =
	std::min(result_type(_M_b1.max() - _M_b1.min()),
		 __detail::_Shift<result_type, __w - __s1>::__value - 1);

const result_type __m2 =
	std::min(result_type(_M_b2.max() - _M_b2.min()),
		 __detail::_Shift<result_type, __w - __s2>::__value - 1);

// NB: In TR1 s1 is not required to be >= s2.
      if (__s1 < __s2)
	_M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
      else
	_M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
    }

template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2>
    typename xor_combine<_UniformRandomNumberGenerator1, __s1,
			 _UniformRandomNumberGenerator2, __s2>::result_type
    xor_combine<_UniformRandomNumberGenerator1, __s1,
		_UniformRandomNumberGenerator2, __s2>::
    _M_initialize_max_aux(result_type __a, result_type __b, int __d)
    {
      const result_type __two2d = result_type(1) << __d;
      const result_type __c = __a * __two2d;

if (__a == 0 || __b < __two2d)
	return __c + __b;

const result_type __t = std::max(__c, __b);
      const result_type __u = std::min(__c, __b);

result_type __ub = __u;
      result_type __p;
      for (__p = 0; __ub != 1; __ub >>= 1)
	++__p;

const result_type __two2p = result_type(1) << __p;
      const result_type __k = __t / __two2p;

if (__k & 1)
	return (__k + 1) * __two2p - 1;

if (__c >= __b)
	return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
							   / __two2d,
							   __u % __two2p, __d);
      else
	return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
							   / __two2d,
							   __t % __two2p, __d);
    }

template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const xor_combine<_UniformRandomNumberGenerator1, __s1,
	       _UniformRandomNumberGenerator2, __s2>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

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

__os.flags(__flags);
      __os.fill(__fill);
      return __os; 
    }

template<class _UniformRandomNumberGenerator1, int __s1,
	   class _UniformRandomNumberGenerator2, int __s2,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       xor_combine<_UniformRandomNumberGenerator1, __s1,
	       _UniformRandomNumberGenerator2, __s2>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

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

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

__is.flags(__flags);
      return __is;
    }

template<typename _IntType>
    template<typename _UniformRandomNumberGenerator>
      typename uniform_int<_IntType>::result_type
      uniform_int<_IntType>::
      _M_call(_UniformRandomNumberGenerator& __urng,
	      result_type __min, result_type __max, true_type)
      {
	// XXX Must be fixed to work well for *arbitrary* __urng.max(),
	// __urng.min(), __max, __min.  Currently works fine only in the
	// most common case __urng.max() - __urng.min() >= __max - __min,
	// with __urng.max() > __urng.min() >= 0.
	typedef typename __gnu_cxx::__add_unsigned<typename
	  _UniformRandomNumberGenerator::result_type>::__type __urntype;
	typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
	                                                      __utype;
	typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
							> sizeof(__utype)),
	  __urntype, __utype>::__type                         __uctype;

result_type __ret;

const __urntype __urnmin = __urng.min();
	const __urntype __urnmax = __urng.max();
	const __urntype __urnrange = __urnmax - __urnmin;
	const __uctype __urange = __max - __min;
	const __uctype __udenom = (__urnrange <= __urange
				   ? 1 : __urnrange / (__urange + 1));
	do
	  __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
	while (__ret > __max - __min);

return __ret + __min;
      }

template<typename _IntType, typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const uniform_int<_IntType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __os.flags();
      const _CharT __fill = __os.fill();
      const _CharT __space = __os.widen(' ');
      __os.flags(__ios_base::scientific | __ios_base::left);
      __os.fill(__space);

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

__os.flags(__flags);
      __os.fill(__fill);
      return __os;
    }

template<typename _IntType, typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       uniform_int<_IntType>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

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

__is.flags(__flags);
      return __is;
    }

template<typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const bernoulli_distribution& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __os.flags();
      const _CharT __fill = __os.fill();
      const std::streamsize __precision = __os.precision();
      __os.flags(__ios_base::scientific | __ios_base::left);
      __os.fill(__os.widen(' '));
      __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);

__os << __x.p();

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _IntType, typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename geometric_distribution<_IntType, _RealType>::result_type
      geometric_distribution<_IntType, _RealType>::
      operator()(_UniformRandomNumberGenerator& __urng)
      {
	// About the epsilon thing see this thread:
        // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
	const _RealType __naf =
	  (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
	// The largest _RealType convertible to _IntType.
	const _RealType __thr =
	  std::numeric_limits<_IntType>::max() + __naf;

_RealType __cand;
	do
	  __cand = std::ceil(std::log(__urng()) / _M_log_p);
	while (__cand >= __thr);

return result_type(__cand + __naf);
      }

template<typename _IntType, typename _RealType,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const geometric_distribution<_IntType, _RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

__os << __x.p();

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _IntType, typename _RealType>
    void
    poisson_distribution<_IntType, _RealType>::
    _M_initialize()
    {
#if _GLIBCXX_USE_C99_MATH_TR1
      if (_M_mean >= 12)
	{
	  const _RealType __m = std::floor(_M_mean);
	  _M_lm_thr = std::log(_M_mean);
	  _M_lfm = std::tr1::lgamma(__m + 1);
	  _M_sm = std::sqrt(__m);

const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
	  const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
							      / __pi_4));
	  _M_d = std::tr1::round(std::max(_RealType(6),
					  std::min(__m, __dx)));
	  const _RealType __cx = 2 * __m + _M_d;
	  _M_scx = std::sqrt(__cx / 2);
	  _M_1cx = 1 / __cx;

_M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
	  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
	}
      else
#endif
	_M_lm_thr = std::exp(-_M_mean);
      }

/**
   * A rejection algorithm when mean >= 12 and a simple method based
   * upon the multiplication of uniform random variates otherwise.
   * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
   * is defined.
   *
   * Reference:
   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
   * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
   */
  template<typename _IntType, typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename poisson_distribution<_IntType, _RealType>::result_type
      poisson_distribution<_IntType, _RealType>::
      operator()(_UniformRandomNumberGenerator& __urng)
      {
#if _GLIBCXX_USE_C99_MATH_TR1
	if (_M_mean >= 12)
	  {
	    _RealType __x;

// See comments above...
	    const _RealType __naf =
	      (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
	    const _RealType __thr =
	      std::numeric_limits<_IntType>::max() + __naf;

const _RealType __m = std::floor(_M_mean);
	    // sqrt(pi / 2)
	    const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
	    const _RealType __c1 = _M_sm * __spi_2;
	    const _RealType __c2 = _M_c2b + __c1; 
	    const _RealType __c3 = __c2 + 1;
	    const _RealType __c4 = __c3 + 1;
	    // e^(1 / 78)
	    const _RealType __e178 = 1.0129030479320018583185514777512983L;
	    const _RealType __c5 = __c4 + __e178;
	    const _RealType __c = _M_cb + __c5;
	    const _RealType __2cx = 2 * (2 * __m + _M_d);

bool __reject = true;
	    do
	      {
		const _RealType __u = __c * __urng();
		const _RealType __e = -std::log(__urng());

_RealType __w = 0.0;
		
		if (__u <= __c1)
		  {
		    const _RealType __n = _M_nd(__urng);
		    const _RealType __y = -std::abs(__n) * _M_sm - 1;
		    __x = std::floor(__y);
		    __w = -__n * __n / 2;
		    if (__x < -__m)
		      continue;
		  }
		else if (__u <= __c2)
		  {
		    const _RealType __n = _M_nd(__urng);
		    const _RealType __y = 1 + std::abs(__n) * _M_scx;
		    __x = std::ceil(__y);
		    __w = __y * (2 - __y) * _M_1cx;
		    if (__x > _M_d)
		      continue;
		  }
		else if (__u <= __c3)
		  // NB: This case not in the book, nor in the Errata,
		  // but should be ok...
		  __x = -1;
		else if (__u <= __c4)
		  __x = 0;
		else if (__u <= __c5)
		  __x = 1;
		else
		  {
		    const _RealType __v = -std::log(__urng());
		    const _RealType __y = _M_d + __v * __2cx / _M_d;
		    __x = std::ceil(__y);
		    __w = -_M_d * _M_1cx * (1 + __y / 2);
		  }

__reject = (__w - __e - __x * _M_lm_thr
			    > _M_lfm - std::tr1::lgamma(__x + __m + 1));

__reject |= __x + __m >= __thr;

} while (__reject);

return result_type(__x + __m + __naf);
	  }
	else
#endif
	  {
	    _IntType     __x = 0;
	    _RealType __prod = 1.0;

do
	      {
		__prod *= __urng();
		__x += 1;
	      }
	    while (__prod > _M_lm_thr);

return __x - 1;
	  }
      }

template<typename _IntType, typename _RealType,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const poisson_distribution<_IntType, _RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __os.flags();
      const _CharT __fill = __os.fill();
      const std::streamsize __precision = __os.precision();
      const _CharT __space = __os.widen(' ');
      __os.flags(__ios_base::scientific | __ios_base::left);
      __os.fill(__space);
      __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);

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

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _IntType, typename _RealType,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       poisson_distribution<_IntType, _RealType>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

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

__is >> __x._M_mean >> __x._M_nd;
      __x._M_initialize();

__is.flags(__flags);
      return __is;
    }

template<typename _IntType, typename _RealType>
    void
    binomial_distribution<_IntType, _RealType>::
    _M_initialize()
    {
      const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;

_M_easy = true;

#if _GLIBCXX_USE_C99_MATH_TR1
      if (_M_t * __p12 >= 8)
	{
	  _M_easy = false;
	  const _RealType __np = std::floor(_M_t * __p12);
	  const _RealType __pa = __np / _M_t;
	  const _RealType __1p = 1 - __pa;
	  
	  const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
	  const _RealType __d1x =
	    std::sqrt(__np * __1p * std::log(32 * __np
					     / (81 * __pi_4 * __1p)));
	  _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
	  const _RealType __d2x =
	    std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
					     / (__pi_4 * __pa)));
	  _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
	  
	  // sqrt(pi / 2)
	  const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
	  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
	  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
	  _M_c = 2 * _M_d1 / __np;
	  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
	  const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
	  const _RealType __s1s = _M_s1 * _M_s1;
	  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
			     * 2 * __s1s / _M_d1
			     * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
	  const _RealType __s2s = _M_s2 * _M_s2;
	  _M_s = (_M_a123 + 2 * __s2s / _M_d2
		  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
	  _M_lf = (std::tr1::lgamma(__np + 1)
		   + std::tr1::lgamma(_M_t - __np + 1));
	  _M_lp1p = std::log(__pa / __1p);

_M_q = -std::log(1 - (__p12 - __pa) / __1p);
	}
      else
#endif
	_M_q = -std::log(1 - __p12);
    }

template<typename _IntType, typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename binomial_distribution<_IntType, _RealType>::result_type
      binomial_distribution<_IntType, _RealType>::
      _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
      {
	_IntType    __x = 0;
	_RealType __sum = 0;

do
	  {
	    const _RealType __e = -std::log(__urng());
	    __sum += __e / (__t - __x);
	    __x += 1;
	  }
	while (__sum <= _M_q);

return __x - 1;
      }

/**
   * A rejection algorithm when t * p >= 8 and a simple waiting time
   * method - the second in the referenced book - otherwise.
   * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
   * is defined.
   *
   * Reference:
   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
   * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
   */
  template<typename _IntType, typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename binomial_distribution<_IntType, _RealType>::result_type
      binomial_distribution<_IntType, _RealType>::
      operator()(_UniformRandomNumberGenerator& __urng)
      {
	result_type __ret;
	const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;

#if _GLIBCXX_USE_C99_MATH_TR1
	if (!_M_easy)
	  {
	    _RealType __x;

// See comments above...
	    const _RealType __naf =
	      (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
	    const _RealType __thr =
	      std::numeric_limits<_IntType>::max() + __naf;

const _RealType __np = std::floor(_M_t * __p12);
	    const _RealType __pa = __np / _M_t;

// sqrt(pi / 2)
	    const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
	    const _RealType __a1 = _M_a1;
	    const _RealType __a12 = __a1 + _M_s2 * __spi_2;
	    const _RealType __a123 = _M_a123;
	    const _RealType __s1s = _M_s1 * _M_s1;
	    const _RealType __s2s = _M_s2 * _M_s2;

bool __reject;
	    do
	      {
		const _RealType __u = _M_s * __urng();

_RealType __v;

if (__u <= __a1)
		  {
		    const _RealType __n = _M_nd(__urng);
		    const _RealType __y = _M_s1 * std::abs(__n);
		    __reject = __y >= _M_d1;
		    if (!__reject)
		      {
			const _RealType __e = -std::log(__urng());
			__x = std::floor(__y);
			__v = -__e - __n * __n / 2 + _M_c;
		      }
		  }
		else if (__u <= __a12)
		  {
		    const _RealType __n = _M_nd(__urng);
		    const _RealType __y = _M_s2 * std::abs(__n);
		    __reject = __y >= _M_d2;
		    if (!__reject)
		      {
			const _RealType __e = -std::log(__urng());
			__x = std::floor(-__y);
			__v = -__e - __n * __n / 2;
		      }
		  }
		else if (__u <= __a123)
		  {
		    const _RealType __e1 = -std::log(__urng());		    
		    const _RealType __e2 = -std::log(__urng());

const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
		    __x = std::floor(__y);
		    __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
					    -__y / (2 * __s1s)));
		    __reject = false;
		  }
		else
		  {
		    const _RealType __e1 = -std::log(__urng());		    
		    const _RealType __e2 = -std::log(__urng());

const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
		    __x = std::floor(-__y);
		    __v = -__e2 - _M_d2 * __y / (2 * __s2s);
		    __reject = false;
		  }

__reject = __reject || __x < -__np || __x > _M_t - __np;
		if (!__reject)
		  {
		    const _RealType __lfx =
		      std::tr1::lgamma(__np + __x + 1)
		      + std::tr1::lgamma(_M_t - (__np + __x) + 1);
		    __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
		  }

__reject |= __x + __np >= __thr;
	      }
	    while (__reject);

__x += __np + __naf;

const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 
	    __ret = _IntType(__x) + __z;
	  }
	else
#endif
	  __ret = _M_waiting(__urng, _M_t);

if (__p12 != _M_p)
	  __ret = _M_t - __ret;
	return __ret;
      }

template<typename _IntType, typename _RealType,
	   typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const binomial_distribution<_IntType, _RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

__os << __x.t() << __space << __x.p() 
	   << __space << __x._M_nd;

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _IntType, typename _RealType,
	   typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       binomial_distribution<_IntType, _RealType>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

__is >> __x._M_t >> __x._M_p >> __x._M_nd;
      __x._M_initialize();

__is.flags(__flags);
      return __is;
    }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const uniform_real<_RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

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

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       uniform_real<_RealType>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

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

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

__is.flags(__flags);
      return __is;
    }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const exponential_distribution<_RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

__os << __x.lambda();

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

/**
   * Polar method due to Marsaglia.
   *
   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
   * New York, 1986, Ch. V, Sect. 4.4.
   */
  template<typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename normal_distribution<_RealType>::result_type
      normal_distribution<_RealType>::
      operator()(_UniformRandomNumberGenerator& __urng)
      {
	result_type __ret;

if (_M_saved_available)
	  {
	    _M_saved_available = false;
	    __ret = _M_saved;
	  }
	else
	  {
	    result_type __x, __y, __r2;
	    do
	      {
		__x = result_type(2.0) * __urng() - 1.0;
		__y = result_type(2.0) * __urng() - 1.0;
		__r2 = __x * __x + __y * __y;
	      }
	    while (__r2 > 1.0 || __r2 == 0.0);

const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
	    _M_saved = __x * __mult;
	    _M_saved_available = true;
	    __ret = __y * __mult;
	  }
	
	__ret = __ret * _M_sigma + _M_mean;
	return __ret;
      }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const normal_distribution<_RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

__os << __x._M_saved_available << __space
	   << __x.mean() << __space
	   << __x.sigma();
      if (__x._M_saved_available)
	__os << __space << __x._M_saved;

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_istream<_CharT, _Traits>&
    operator>>(std::basic_istream<_CharT, _Traits>& __is,
	       normal_distribution<_RealType>& __x)
    {
      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
      typedef typename __istream_type::ios_base    __ios_base;

const typename __ios_base::fmtflags __flags = __is.flags();
      __is.flags(__ios_base::dec | __ios_base::skipws);

__is >> __x._M_saved_available >> __x._M_mean
	   >> __x._M_sigma;
      if (__x._M_saved_available)
	__is >> __x._M_saved;

__is.flags(__flags);
      return __is;
    }

template<typename _RealType>
    void
    gamma_distribution<_RealType>::
    _M_initialize()
    {
      if (_M_alpha >= 1)
	_M_l_d = std::sqrt(2 * _M_alpha - 1);
      else
	_M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
		  * (1 - _M_alpha));
    }

/**
   * Cheng's rejection algorithm GB for alpha >= 1 and a modification
   * of Vaduva's rejection from Weibull algorithm due to Devroye for
   * alpha < 1.
   *
   * References:
   * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral
   * Shape Parameter. Applied Statistics, 26, 71-75, 1977.
   *
   * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection
   * and Composition Procedures. Math. Operationsforschung and Statistik,
   * Series in Statistics, 8, 545-576, 1977.
   *
   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
   * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
   */
  template<typename _RealType>
    template<class _UniformRandomNumberGenerator>
      typename gamma_distribution<_RealType>::result_type
      gamma_distribution<_RealType>::
      operator()(_UniformRandomNumberGenerator& __urng)
      {
	result_type __x;

bool __reject;
	if (_M_alpha >= 1)
	  {
	    // alpha - log(4)
	    const result_type __b = _M_alpha
	      - result_type(1.3862943611198906188344642429163531L);
	    const result_type __c = _M_alpha + _M_l_d;
	    const result_type __1l = 1 / _M_l_d;

// 1 + log(9 / 2)
	    const result_type __k = 2.5040773967762740733732583523868748L;

do
	      {
		const result_type __u = __urng();
		const result_type __v = __urng();

const result_type __y = __1l * std::log(__v / (1 - __v));
		__x = _M_alpha * std::exp(__y);

const result_type __z = __u * __v * __v;
		const result_type __r = __b + __c * __y - __x;

__reject = __r < result_type(4.5) * __z - __k;
		if (__reject)
		  __reject = __r < std::log(__z);
	      }
	    while (__reject);
	  }
	else
	  {
	    const result_type __c = 1 / _M_alpha;

do
	      {
		const result_type __z = -std::log(__urng());
		const result_type __e = -std::log(__urng());

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

__reject = __z + __e < _M_l_d + __x;
	      }
	    while (__reject);
	  }

return __x;
      }

template<typename _RealType, typename _CharT, typename _Traits>
    std::basic_ostream<_CharT, _Traits>&
    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
	       const gamma_distribution<_RealType>& __x)
    {
      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
      typedef typename __ostream_type::ios_base    __ios_base;

__os << __x.alpha();

__os.flags(__flags);
      __os.fill(__fill);
      __os.precision(__precision);
      return __os;
    }
}

_GLIBCXX_END_NAMESPACE_VERSION
}

#endif

cmd: