// random number generation -*- C++ -*-

// Copyright (C) 2005 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 2, 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.

// You should have received a copy of the GNU General Public License along
// with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the GNU General Public License.

#ifndef _STD_TR1_RANDOM
#define _STD_TR1_RANDOM 1

/**
 * @file 
 * This is a TR1 C++ Library header. 
 */

#include <algorithm>
#include <bits/concept_check.h>
#include <bits/cpp_type_traits.h>
#include <cmath>
#include <debug/debug.h>
#include <iterator>
#include <iosfwd>
#include <limits>
#include <tr1/type_traits>


namespace std {
namespace tr1 {
	// [5.1] Random number generation

	/**
	 * @addtogroup tr1_random Random Number Generation
	 * A facility for generating random numbers on selected distributions.
	 * @{
	 */

	/*
	 * Implementation-space details.
	 */
	namespace _Private
	{
		// Type selectors -- are these already implemeted elsewhere?
		template<bool, typename _TpTrue, typename _TpFalse>
			struct _Select
			{
				typedef _TpTrue Type;
			};

		template<typename _TpTrue, typename _TpFalse>
			struct _Select<false, _TpTrue, _TpFalse>
			{
				typedef _TpFalse Type;
			};

	/*
	 * An adaptor class for converting the output of any Generator into
	 * the input for a specific Distribution.
	 */
	template<typename _Generator, typename _Distribution>
		struct _Adaptor
		{ 
			typedef typename _Generator::result_type   generated_type;
			typedef typename _Distribution::input_type result_type;

		public:
			_Adaptor(const _Generator& __g)
			: _M_g(__g)
			{ }

			result_type
			operator()();

		private:
			_Generator _M_g;
		};

	/*
	 * Converts a value generated by the adapted random number genereator into a
	 * value in the input domain for the dependent random number distribution.
	 *
	 * Because the type traits are compile time constants only the appropriate
	 * clause of the if statements will actually be emitted by the compiler.
	 */
	template<typename _Generator, typename _Distribution>
		typename _Adaptor<_Generator, _Distribution>::result_type
		_Adaptor<_Generator, _Distribution>::
		operator()()
		{
			result_type __return_value = 0;
			if (is_integral<generated_type>::value
			 && is_integral<result_type>::value)
			{
				__return_value = _M_g();
			}
			else if (is_integral<generated_type>::value
					 && !is_integral<result_type>::value)
			{
				__return_value = result_type(_M_g())
					               / result_type(_M_g.max() - _M_g.min() + 1);
			}
			else if (!is_integral<generated_type>::value
					  && !is_integral<result_type>::value)
			{
				__return_value = result_type(_M_g())
					               / result_type(_M_g.max() - _M_g.min());
			}
			return __return_value;
		}
	} // namespace std::tr1::_Private


	/**
	 * Produces random numbers on a given disribution function using a un uniform
	 * random number generation engine.
	 *
	 * @todo the engine_value_type needs to be studied more carefully.
	 */
	template<typename _Generator, typename _Dist>
		class variate_generator
		{
      // Concept requirements.
      __glibcxx_class_requires(_Generator, _CopyConstructibleConcept)
//      __glibcxx_class_requires(_Generator, _GeneratorConcept)
//      __glibcxx_class_requires(_Dist,      _GeneratorConcept)

		public:
			typedef _Generator                             engine_type;
			typedef _Private::_Adaptor<_Generator, _Dist>  engine_value_type;
			typedef _Dist                                  distribution_type;
			typedef typename _Dist::result_type            result_type;
			
			// tr1:5.1.1 table 5.1 requirement
			typedef typename std::__enable_if<
			                         result_type,
			                         is_arithmetic<result_type>::value
			                      >::__type _IsValidType;

		public:
			/**
			 * Constructs a variate generator with the uniform random number
			 * generator @p eng for the random distribution @p d.
			 *
			 * @throws Any exceptions which may thrown by the copy constructors of the
			 * @p _Generator or @p _Dist objects.
			 */
			variate_generator(engine_type __eng, distribution_type __dist)
			: _M_engine(__eng), _M_dist(__dist)
			{ }

			/**
			 * Gets the next generated value on the distribution.
			 */
			result_type
			operator()();

			template<typename _Tp>
				result_type
				operator()(_Tp __value);

			/**
			 * Gets a reference to the underlying uniform random number generator
			 * object.
			 */
			engine_value_type&
			engine()
			{ return _M_engine; }

			/**
			 * Gets a const reference to the underlying uniform random number
			 * generator object.
			 */
			const engine_value_type&
			engine() const
			{ return _M_engine; }

			/**
			 * Gets a reference to the underlying random distribution.
			 */
			distribution_type&
			distribution()
			{ return _M_dist; }

			/**
			 * Gets a const reference to the underlying random distribution.
			 */
			const distribution_type&
			distribution() const
			{ return _M_dist; }

			/**
			 * Gets the closed lower bound of the distribution interval.
			 */
			result_type
			min() const
			{ return this->distribution().min(); }

			/**
			 * Gets the closed upper bound of the distribution interval.
			 */
			result_type
			max() const
			{ return this->distribution().max(); }

		private:
			engine_value_type _M_engine;
			distribution_type _M_dist;
		};

	/**
	 * Gets the next random value on the given distribution.
	 */
	template<typename _Generator, typename _Dist>
		typename variate_generator<_Generator, _Dist>::result_type
		variate_generator<_Generator, _Dist>::
		operator()()
		{ return _M_dist(_M_engine); }

	/**
	 * WTF?
	 */
	template<typename _Generator, typename _Dist>
		template<typename _Tp>
			typename variate_generator<_Generator, _Dist>::result_type
			variate_generator<_Generator, _Dist>::
			operator()(_Tp __value)
			{ return _M_dist(_M_engine, __value); }


	/**
	 * @addtogroup tr1_random_generators Random Number Generators
	 * @ingroup tr1_random
	 *
	 * These classes define objects which provide random or pseudorandom numbers,
	 * either from a discrete or a continuous interval.  The random number
	 * generator supplied as a part of this library are all uniform random number
	 * generators which provide a sequence of random number uniformly distributed
	 * over their range.
	 *
	 * A number generator is a function object with an operator() that takes zero
	 * arguments and returns a number.
	 *
	 * A compliant random number generator must satisy the following requirements.
	 * <table border=1 cellpadding=10 cellspacing=0>
	 * <caption align=top>Random Number Generator Requirements</caption>
	 * <tr><td>To be documented.</td></tr>
	 * </table>
	 * 
	 * @{
	 */

	/**
	 * @brief A model of a linear congruential random number generator.
	 *
	 * A random number generator that produces pseudorandom numbers using the
	 * linear function @f$x_{i+1}\leftarrow(ax_{i} + c) \bmod m @f$.
	 *
	 * The template parameter @p UIntType must be an unsigned integral type large
	 * enough to store values up to (m-1). If the template parameter @p m is 0,
	 * the modulus @p m used is std::numeric_limits<UIntType>::max() plus 1.
	 * Otherwise, the template parameters @p a and @p c must be less than @p m.
	 *
	 * The size of the state is @f$ 1 @f$.
	 */
	template<class UIntType, UIntType a, UIntType c, UIntType m>
		class linear_congruential
		{
			__glibcxx_class_requires(UIntType, _UnsignedIntegerConcept);
//			__glibcpp_class_requires(a < m && c < m)

		public:
			/** The type of the generated random value. */
			typedef UIntType result_type;

			/** The multiplier. */
			static const UIntType multiplier = a;
			/** An increment. */
			static const UIntType increment = c;
			/** The modulus. */
			static const UIntType modulus = m;

			/**
			 * Constructs a %linear_congruential random number generator engine with
			 * seed @p s.  The default seed value is 1.
			 *
			 * @param s The initial seed value.
			 */
			explicit linear_congruential(unsigned long s = 1);

			/**
			 * Constructs a %linear_congruential random number generator engine
			 * seeded from the generator function @g.
			 *
			 * @param g The seed generator function.
			 */
			template<class Gen>
				linear_congruential(Gen& g);

			/**
			 * Resets the %linear_congruential random number generator engine sequence
			 * to @p.
			 *
			 * @param s The new seed.
			 */
			void
			seed(unsigned long s = 1);

			/**
			 * Resets the %linear_congruential random number generator engine sequence
			 * using a vlaue from the generator function @g.
			 *
			 * @param g the seed generator function.
			 */
			template<class Gen>
				void
				seed(Gen& g)
				{ seed(g, typename is_fundamental<Gen>::type()); }

			/**
			 * Gets the smallest possible value in the output range.
			 */
			result_type
			min() const;

			/**
			 * Gets the largest possible value in the output range.
			 */
			result_type
			max() const;

			/**
			 * Gets the next random number in the sequence.
			 */
			result_type
			operator()();

		 /**
		  * Compares two linear congruential random number generator objects of the
		  * same type for equality.
			*  
		  * @param lhs A linear congruential random number generator object.
		  * @param rhs Another linear congruential random number generator object.
		  *
		  * @returns true if the two objects are equal, false otherwise.
		  */
			friend bool
			operator==(const linear_congruential& lhs,
								 const linear_congruential& rhs)
			{ return lhs.m_x == rhs.m_x; }

			/**
			 * Compares two linear congruential random number generator objects of the
			 * same type for inequality.
			 *
			 * @param lhs A linear congruential random number generator object.
			 * @param rhs Another linear congruential random number generator object.
			 *
			 * @returns true if the two objects are not equal, false otherwise.
			 */
			friend bool
			operator!=(const linear_congruential& lhs,
								 const linear_congruential& rhs)
			{ return !(lhs == rhs); }

			/**
			 * Writes the textual representation of the state x(i) of x to @p os.
			 *
			 * @param os  The output stream.
			 * @param lcr A linear_congruential random number generator.
			 * @returns os.
			 */
			template<typename CharT, typename Traits>
				friend std::basic_ostream<CharT, Traits>&
				operator<<(std::basic_ostream<CharT, Traits>& os,
									 const linear_congruential& lcr)
				{ return os << lcr.m_x; }

			/**
			 * Sets the state of the engine by reading its textual
			 * representation from @p is.
			 *
			 * The textual representation must have been previously written using an
			 * output stream whose imbued locale and whose type's template
			 * specialization arguments CharT and Traits were the same as those of
			 * @p is.
			 *
			 * @param is  The input stream.
			 * @param lcr A linear_congruential random number generator.
			 * @returns os.
			 */
			template<typename CharT, typename Traits>
				friend std::basic_istream<CharT, Traits>&
				operator>>(std::basic_istream<CharT, Traits>& is,
									 linear_congruential& lcr)
				{ return is >> lcr.m_x; }

		private:
			template<class Gen>
				void
				seed(Gen& g, true_type)
				{ return seed(static_cast<unsigned long>(g)); }

			template<class Gen>
				void
				seed(Gen& g, false_type);

		private:
			UIntType m_x;
		};

	/**
	 * The classic Minimum Standard rand0 of Lewis, Goodman, and Miller.
	 */
	typedef linear_congruential<unsigned int, 16807, 0, 2147483647> minstd_rand0;

	/**
	 * An alternative LCR (Lehmer Generator function) .
	 */
	typedef linear_congruential<unsigned int, 48271, 0, 2147483647> minstd_rand;


	/**
	 * A generalized feedback shift register discrete random number generator.
	 *
	 * This algorithm avoind multiplication and division and is designed to be
	 * friendly to a pipelined architecture.  If the parameters are chosen
	 * correctly, this generator will produce numbers with a very long period and
	 * fairly good apparent entropy, although still not cryptographically strong.
	 *
	 * The best way to use this generator is with the predefined mt19937 class.
	 *
	 * This algorithm was originally invented by Makoto Matsumoto and
	 * Takuji Nishimura.
	 *
	 * @var word_size   The number of bits in each element of the state vector.
	 * @var state_size  The degree of recursion.
	 * @var shift_size  The period parameter.
	 * @var mask_bits   The separation point bit index.
	 * @var parameter_a The last row of the twist matrix.
	 * @var output_u    The first right-shift tempering matrix parameter.
	 * @var output_s    The first left-shift tempering matrix parameter.
	 * @var output_b    The first left-shift tempering matrix mask.
	 * @var output_t    The second left-shift tempering matrix parameter.
	 * @var output_c    The second left-shift tempering matrix mask.
	 * @var output_l    The second right-shift tempering matrix parameter.
	 */
	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>
		class mersenne_twister
		{
			__glibcxx_class_requires(UIntType, _UnsignedIntegerConcept);

		public:
			// types
			typedef UIntType result_type ;

			// parameter values
			static const int      word_size   = w;
			static const int      state_size  = n;
			static const int      shift_size  = m;
			static const int      mask_bits   = r;
			static const UIntType parameter_a = a;
			static const int      output_u    = u;
			static const int      output_s    = s;
			static const UIntType output_b    = b;
			static const int      output_t    = t;
			static const UIntType output_c    = c;
			static const int      output_l    = l;
			
			// constructors and member function
			mersenne_twister()
			{ seed(); }
			
			explicit
			mersenne_twister(unsigned long value)
			{ seed(value); }

			template<class Gen>
				mersenne_twister(Gen& g)
				{ seed(g); }

			void
			seed()
			{ seed(0UL); }

			void
			seed(unsigned long value);

			template<class Gen>
				void
				seed(Gen& g)
				{ seed(g, typename is_fundamental<Gen>::type()); }

			result_type
			min() const
			{ return 0; };

			result_type
			max() const;

			result_type
			operator()();

		private:
			template<class Gen>
				void
				seed(Gen& g, true_type)
				{ return seed(static_cast<unsigned long>(g)); }

			template<class Gen>
				void
				seed(Gen& g, false_type);

		private:
			UIntType _M_x[state_size];
			int      _M_p;
		};

	/**
	 * The classic Mersenne Twister.
	 *
	 * Reference:
	 * M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
	 * Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions
	 * on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. 	 */
	typedef mersenne_twister<
		unsigned long, 32, 624, 397, 31,
		0x9908b0dful, 11, 7,
		0x9d2c5680ul, 15,
		0xefc60000ul, 18
	> mt19937;


	/**
	 * @brief The Marsaglia-Zaman generator.
	 * 
	 * This is a model of a Generalized Fibonacci discrete random number
	 * generator, sometimes referred to as the SWC generator.
	 *
	 * A discrete random number generator that produces pseudorandom numbers using
	 * @f$x_{i}\leftarrow(x_{i - s} - x_{i - r} - carry_{i-1}) \bmod m @f$.
	 *
	 * The size of the state is @f$ r @f$
	 * and the maximum period of the generator is @f$ m^r - m^s -1 @f$.
	 *
	 * N1688[4.13] says "the template parameter _IntType shall denote an integral
	 * type large enough to store values up to m."
	 *
	 * @if maint
	 * @var _M_x     The state of te generator.  This is a ring buffer.
	 * @var _M_carry The carry.
	 * @var _M_p     Current index of x(i - r).
	 * @endif
	 */
	template<typename _IntType, _IntType m, int s, int r>
		class subtract_with_carry
		{
			__glibcxx_class_requires(_IntType, _IntegerConcept)

		public:
			/** The type of the generated random value. */
			typedef _IntType result_type;
			
			// parameter values
			static const _IntType modulus   = m;
			static const int      long_lag  = r;
			static const int      short_lag = s;
			
		public:
			/**
			 * Constructs a default-initialized % subtract_with_carry random number
			 * generator.
			 */
			subtract_with_carry()
			{ this->seed(); }

			/**
			 * Constructs an explicitly seeded % subtract_with_carry random number
			 * generator.
			 */
			explicit
			subtract_with_carry(_IntType __value)
			{ this->seed(__value); }

			/**
			 * Constructs a % subtract_with_carry random number generator seeded from
			 * the PAD iterated by [__first, last).
			 */
			template<class Gen>
				subtract_with_carry(Gen& g)
				{ this->seed(g); }

			/**
			 * Seeds the initial state @f$ x_0 @f$ of the random number generator.
			 *
			 * @note This implementation follows the tr1 specification but will
			 * obviously not work correctly on all platforms, since it has hardcoded
			 * values that may overflow ints on some platforms.
			 *
			 * N1688[4.19] modifies this as follows.
			 * If @p __value == 0, sets value to 19780503.  In any case, with a linear
			 * congruential generator lcg(i) having parameters @f$ m_{lcg} =
			 * 2147483563, a_{lcg} = 40014, c_{lcg} = 0, and lcg(0) = value @f$, sets
			 * @f$ x_{-r} \dots x_{-1} @f$ to
			 * @f$ lcg(1) \bmod m \dots lcg(r) \bmod m @f$ respectively.
			 * If @f$ x_{-1} = 0 @f$ set carry to 1, otherwise sets carry to 0.
			 */
			void
			seed(_IntType __value = 19780503);

			/**
			 * Seeds the initial state @f$ x_0 @f$ of the % subtract_with_carry
			 * random number generator.
			 */
			template<class Gen>
				void
				seed(Gen& g)
				{ seed(g, typename is_fundamental<Gen>::type()); }

			/**
			 * Gets the inclusive minimum value of the range of random integers
			 * returned by this generator.
			 */
			result_type
			min() const
			{ return 0; }
			
			/**
			 * Gets the inclusive maximum value of the range of random integers
			 * returned by this generator.
			 */
			result_type
			max() const
			{ return this->modulus - 1; }
			
			/**
			 * Gets the next random number in the sequence.
			 */
			result_type
			operator()();

			/**
			 * Compares two % subtract_with_carry random number generator objects of
			 * the same type for equality.
			 *
			 * @param __lhs A % subtract_with_carry random number generator object.
			 * @param __rhs Another % subtract_with_carry random number generator
			 *              object.
			 *
			 * @returns true if the two objects are equal, false otherwise.
			 */
			friend bool
			operator==(const subtract_with_carry& __lhs,
								 const subtract_with_carry& __rhs)
			{ 
				return ((__lhs._M_x[0] == __rhs._M_x[0])
						    && (__lhs._M_x[r-1] == __rhs._M_x[r-1]));
			}

			/**
			 * Compares two % subtract_with_carry random number generator objects of
			 * the same type for inequality.
			 *
			 * @param __lhs A % subtract_with_carry random number generator object.
			 * @param __rhs Another % subtract_with_carry random number generator
			 *              object.
			 *
			 * @returns true if the two objects are not equal, false otherwise.
			 */
			friend bool
			operator!=(const subtract_with_carry& __lhs,
								 const subtract_with_carry& __rhs)
			{ return !(__lhs == __rhs); }

			/**
			 * Inserts the current state of a % subtract_with_carry random number
			 * genator engine @p x into the output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A % subtract_with_carry random number generator engine.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const subtract_with_carry& __x)
				{
					std::copy(__x._M_x, __x._M_x + r,
					          std::ostream_iterator<_IntType>(__os, " "));
					return __os << __x._M_carry;
				}

			/**
			 * Extracts the current state of a % subtract_with_carry random number
			 * gerator engine @p x from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __x  A % subtract_with_carry random number generator engine.
			 *
			 * @returns The input stream with the state of @p x extracted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 subtract_with_carry& __x)
				{
					for (int __i = 0; __i < r; ++__i)
					{
						__is >> __x._M_x[__i];
					}
					__is >> __x._M_carry;
					return __is;
				}

		private:
			template<class Gen>
				void
				seed(Gen& g, true_type)
				{ return seed(static_cast<unsigned long>(g)); }

			template<class Gen>
				void
				seed(Gen& g, false_type);

		private:
			int         _M_p;
			result_type _M_x[long_lag];
			result_type _M_carry;
		};	


	/**
	 * Produces random numbers from some base engine by discarding blocks of
	 * data.
	 *
	 * 0 <= @p r <= @p p
	 */
	template<class UniformRandomNumberGenerator, int p, int r>
		class discard_block
		{
//			__glibcxx_class_requires(typename base_type::result_type, ArithmeticTypeConcept);

		public:
			/** The type of the underlying generator engine. */
			typedef UniformRandomNumberGenerator    base_type;
			/** The type of the generated random value. */
			typedef typename base_type::result_type result_type;

			// parameter values
			static const int block_size = p;
			static const int used_block = r;

			/**
			 * Constructs a default %discard_block engine.
			 *
			 * The underlying engine is default constrcuted as well.
			 */
			discard_block()
			: _M_n(0)
			{ }

			/**
			 * Copy constructs a %discard_block engine.
			 *
			 * Copies an existing base class random number geenerator.
			 * @param rng An existing (base class) engine object.
			 */
			explicit discard_block(const base_type& rng)
			: _M_b(rng)
			, _M_n(0)
			{ }

			/**
			 * Seed constructs a %discard_block engine.
			 *
			 * Constructs the underlying generator engine seeded with @p s.
			 * @param s A seed value for the base class engine.
			 */
			explicit discard_block(unsigned long s)
			: _M_b(s)
			, _M_n(0)
			{ }

			/**
			 * Generator constructs a %discard_block engine.
			 *
			 * @param g A seed generator function.
			 */
			template<class Gen>
				discard_block(Gen& g)
				: _M_b(g)
				, _M_n(0)
				{ }

			/**
			 * Reseeds the %discard_block object with the default seed for the
			 * underlying base class generator engine.
			 */
			void seed()
			{
				_M_b.seed();
				_M_n = 0;
			}

			/**
			 * Reseeds the %discard_block object with the given seed generator
			 * function.
			 * @param g A seed generator function.
			 */
			template<class Gen>
				void seed(Gen& g)
				{
					_M_b.seed(g);
					_M_n = 0;
				}

			/**
			 * Gets a const reference to the underlying generator engine object.
			 */
			const base_type&
			base() const
			{ return _M_b; }

			/**
			 * Gets the minimum value in the generated random number range.
			 */
			result_type
			min() const
			{ return _M_b.min(); }

			/**
			 * Gets the maximum value in the generated random number range.
			 */
			result_type
			max() const
			{ return _M_b.max(); }

			/**
			 * Gets the next value in the generated random number sequence.
			 */
			result_type
			operator()();

			/**
			 * Compares two %discard_block random number generator objects of
			 * the same type for equality.
			 *
			 * @param __lhs A %discard_block random number generator object.
			 * @param __rhs Another %discard_block random number generator
			 *              object.
			 *
			 * @returns true if the two objects are equal, false otherwise.
			 */
			friend bool
			operator==(const discard_block& __lhs,
								 const discard_block& __rhs)
			{ 
				return ((__lhs._M_b == __rhs._M_b)
						    && (__lhs._M_n == __rhs._M_n));
			}

			/**
			 * Compares two %discard_block random number generator objects of
			 * the same type for inequality.
			 *
			 * @param __lhs A %discard_block random number generator object.
			 * @param __rhs Another %discard_block random number generator
			 *              object.
			 *
			 * @returns true if the two objects are not equal, false otherwise.
			 */
			friend bool
			operator!=(const discard_block& __lhs,
								 const discard_block& __rhs)
			{ return !(__lhs == __rhs); }

			/**
			 * Inserts the current state of a %discard_block random number
			 * genator engine @p x into the output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %discard_block random number generator engine.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const discard_block& __x)
				{ return __os << __x._M_b << " " << __x._M_n; }

			/**
			 * Extracts the current state of a % subtract_with_carry random number
			 * gerator engine @p x from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __x  A %discard_block random number generator engine.
			 *
			 * @returns The input stream with the state of @p x extracted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 discard_block& __x)
				{ return __is >> __x._M_b >> __x._M_n; }
			
		private:
			base_type _M_b;
			int       _M_n;
		};


	/**
	 * James's luxury-level-3 integer adaptation of Luescher's generator.
	 */
	typedef discard_block<
	          subtract_with_carry<int, (1<<24), 10, 24>,
		        223,
						24
					> ranlux3;

	/**
	 * James's luxury-level-4 integer adaptation of Luescher's generator.
	 */
	typedef discard_block<
	          subtract_with_carry<int, (1<<24), 10, 24>,
		        389,
						24
					> ranlux4;


	/**
	 * A random number generator adaptor class that combines two random number
	 * generator engines into a single output sequence.
	 */
	template<class UniformRandomNumberGenerator1, int s1,
	         class UniformRandomNumberGenerator2, int s2>
		class xor_combine
		{
//			__glibcxx_class_requires(typename UniformRandomNumberGenerator1::result_type, ArithmeticTypeConcept);
//			__glibcxx_class_requires(typename UniformRandomNumberGenerator2::result_type, ArithmeticTypeConcept);
		public:
			/** The type of the the first underlying generator engine. */
			typedef UniformRandomNumberGenerator1 base1_type;
			/** The type of the the second underlying generator engine. */
			typedef UniformRandomNumberGenerator2 base2_type;
			/** The type of the generated random value. */
			typedef typename _Private::_Select<
			                             (sizeof(base1_type) > sizeof(base2_type)),
							                     base1_type,
			                             base2_type
			                           >::Type result_type;

			// parameter values
			static const int shift1 = s1;
			static const int shift2 = s2;

			// constructors and member function
			xor_combine()
			{ }

			xor_combine(const base1_type& rng1, const base2_type& rng2)
			: _M_b1(rng1)
			, _M_b2(rng2)
			{ }

			xor_combine(unsigned long s)
			: _M_b1(s)
			, _M_b2(s + 1)
			{ }

			template<class Gen>
				xor_combine(Gen& g)
				: _M_b1(g)
				, _M_b2(g)
				{ }

			void
			seed()
			{
				_M_b1.seed();
				_M_b2.seed();
			}

			template<class Gen> 
				void
				seed(Gen& g)
				{
					_M_b1.seed(g);
					_M_b2.seed(g);
				}

			const base1_type&
			base1() const
			{ return _M_b1; }

			const base2_type&
			base2() const
			{ return _M_b2; }

			result_type
			min() const
			{ return _M_b1.min() ^ _M_b2.min(); }

			result_type
			max() const
			{ return _M_b1.max() | _M_b2.max(); }

			/**
			 * Gets the next random number in the sequence.
			 */
			result_type
			operator()()
			{
				return ((_M_b1() << shift1) ^ (_M_b2() << shift2));
			}

			/**
			 * Compares two %xor_combine random number generator objects of
			 * the same type for equality.
			 *
			 * @param __lhs A %xor_combine random number generator object.
			 * @param __rhs Another %xor_combine random number generator
			 *              object.
			 *
			 * @returns true if the two objects are equal, false otherwise.
			 */
			friend bool
			operator==(const xor_combine& __lhs,
								 const xor_combine& __rhs)
			{ 
				return (__lhs.base1() == __rhs.base1())
						&& (__lhs.base2() == __rhs.base2());
			}

			/**
			 * Compares two %xor_combine random number generator objects of
			 * the same type for inequality.
			 *
			 * @param __lhs A %xor_combine random number generator object.
			 * @param __rhs Another %xor_combine random number generator
			 *              object.
			 *
			 * @returns true if the two objects are not equal, false otherwise.
			 */
			friend bool
			operator!=(const xor_combine& __lhs,
								 const xor_combine& __rhs)
			{ 
				return !(__lhs == __rhs);
			}

			/**
			 * Inserts the current state of a %xor_combine random number
			 * genator engine @p x into the output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %xor_combine random number generator engine.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const xor_combine& __x)
				{ return __os << __x.base1() << " " << __x.base1(); }

			/**
			 * Extracts the current state of a %xor_combine random number
			 * gerator engine @p x from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __x  A %xor_combine random number generator engine.
			 *
			 * @returns The input stream with the state of @p x extracted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 xor_combine& __x)
				{ return __is >> __x._M_b1 >> __x._M_b2; }

		private:
			base1_type _M_b1;
			base2_type _M_b2;
		};


	/**
	 * A standard interface to a platform-specific non-deterministic random number
	 * generator (if any are available).
	 *
	 * @todo The underlying interface is system-specific and needs to be factored
	 * into the generated configury mechs.  For example, the use of "/dev/random"
	 * under a Linux OS would be appropriate.
	 */
	class random_device
	{
		public:
			// types
			typedef unsigned int result_type;

			// constructors, destructors and member functions
			explicit random_device(const std::string& token = "unimplemented");
			result_type min() const;
			result_type max() const;
			double entropy() const;
			result_type operator()();

		private:
			random_device(const random_device &);
			void operator=(const random_device &);
	};

	/* @} */ // group tr1_random_generators

	/**
	 * @addtogroup tr1_random_distributions Random Number Distributions
	 * @ingroup tr1_random
	 * @{
	 */

	/**
	 * @addtogroup tr1_random_distributions_discrete Discrete Distributions
	 * @ingroup tr1_random_distributions
	 * @{
	 */

	/**
	 * @brief Uniform discrete distribution for random numbers.
	 * A discrete random distribution on the range @f$[min, max]@f$ with equal
	 * probability throughout the range.
	 */
	template<typename _IntType = int>
		class uniform_int
		{
			__glibcxx_class_requires(_IntType, _IntegerConcept)

		public:
			/** The type of the parameters of the distribution. */
			typedef _IntType input_type;
			/** The type of the range of the distribution. */
			typedef _IntType result_type;

		public:
			/**
			 * Constructs a uniform distribution object.
			 */
			explicit
			uniform_int(_IntType __min = 0, _IntType __max = 9)
			: _M_min(__min)
			, _M_max(__max)
			{
				_GLIBCXX_DEBUG_ASSERT(_M_min <= _M_max);
			}

			/**
			 * Gets the inclusive lower bound of the distribution range.
			 */
			result_type
			min() const
			{ return _M_min; }

			/**
			 * Gets the inclusive upper bound of the distribution range.
			 */
			result_type
			max() const
			{ return _M_max; }

			/**
			 * Resets the distribution state.
			 *
			 * Does nothing for the uniform integer distribution.
			 */
			void
			reset()
			{ }

			/**
			 * Gets a uniformly distributed random number in the range
			 * @f$(min, max)@f$.
			 */
			template<typename _UniformRandomNumberGenerator>
				result_type
				operator()(_UniformRandomNumberGenerator& __urng)
				{ return (__urng() % (_M_max - _M_min + 1)) + _M_min; }

			/**
			 * Gets a uniform random number in the range @f$[0, n)@f$.
			 *
			 * This function is aimed at use with std::random_shuffle.
			 */
			template<typename _UniformRandomNumberGenerator>
				result_type
				operator()(_UniformRandomNumberGenerator& __urng, result_type __n)
				{
					return __urng() % __n;
				}

			/**
			 * Inserts a %uniform_int random number distribution @p x into the
			 * output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %uniform_int random number distribution.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const uniform_int& __x)
				{ return __os << __x._M_min << " " << __x._M_max; }

			/**
			 * Extracts a %unform_int random number distribution
			 * @p u from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __u  A %uniform_int random number generator engine.
			 *
			 * @returns The input stream with @p u extracted or in an error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 uniform_int& __u)
				{ return __is >> __u._M_min >> __u._M_max; }

		private:
			_IntType _M_min;
			_IntType _M_max;
		};


	/**
	 * @brief A Bernoulli random number distribution.
	 *
	 * Generates a sequence of true and false values with likelihood @f$ p @f$
	 * that true will come up and @f$ (1 - p) @f$ that false will appear.
	 */
	class bernoulli_distribution
	{
	public:
		typedef int  input_type;
		typedef bool result_type;
		
	public:
		/**
		 * Constructs a Bernoulli distribution with likelihood @p p.
		 *
		 * @param p  [IN]  The likelihood of a true result being returned.  Must
		 * be in the interval @f$ [0, 1] @f$.
		 */
		explicit
		bernoulli_distribution(double __p = 0.5)
		: _M_p(__p)
		{ 
			_GLIBCXX_DEBUG_ASSERT((_M_p >= 0.0) && (_M_p <= 1.0));
		}
		
		/**
		 * Gets the @p p parameter of the distribution.
		 */
		double
		p() const
		{ return _M_p; }
		
		/**
		 * Gets the inclusive lower bound of the distribution range.
		 */
		result_type
		min() const
		{ return false; }

		/**
		 * Gets the inclusive upper bound of the distribution range.
		 */
		result_type
		max() const
		{ return true; }

		/**
		 * Resets the distribution state.
		 *
		 * Does nothing for a bernoulli distribution.
		 */
		void
		reset()
		{ }
		
		/**
		 * Gets the next value in the Bernoullian sequence.
		 */
		template<class UniformRandomNumberGenerator>
			result_type
			operator()(UniformRandomNumberGenerator& __urng)
			{
				if (__urng() < _M_p)
				{
					return true;
				}
				return false;
			}

		/**
		 * Inserts a %bernoulli_distribution random number distribution
		 * @p x into the output stream @p os.
		 *
		 * @param __os An output stream.
		 * @param __x  A %bernoulli_distribution random number distribution.
		 *
		 * @returns The output stream with the state of @p x inserted or in an
		 * error state.
		 */
		template<typename _CharT, typename _Traits>
			friend basic_ostream<_CharT, _Traits>&
			operator<<(basic_ostream<_CharT, _Traits>& __os,
								 const bernoulli_distribution& __x)
			{ return __os << __x.p(); }

		/**
		 * Extracts a %bernoulli_distribution random number distribution
		 * @p u from the input stream @p is.
		 *
		 * @param __is An input stream.
		 * @param __u  A %bernoulli_distribution random number generator engine.
		 *
		 * @returns The input stream with @p u extracted or in an error state.
		 */
		template<typename _CharT, typename _Traits>
			friend basic_istream<_CharT, _Traits>&
			operator>>(basic_istream<_CharT, _Traits>& __is,
								 bernoulli_distribution& __u)
			{ return __is >> __u._M_p; }

	protected:
		double _M_p;
	};


	/**
	 * @brief A discrete geometric random number distribution.
	 *
	 * The formula for the geometric probability mass function is 
	 * @f$ p(i) = (1 - p)p^{i-1} @f$ where @f$ p @f$ is the parameter of the
	 * distribution.
	 */
	template<typename _IntType = int, typename _RealType = double>
		class geometric_distribution
		{
		public:
			// types
			typedef _RealType input_type;
			typedef _IntType  result_type;
			
			// constructors and member function

			explicit
			geometric_distribution(const _RealType& __p = _RealType(0.5))
			: _M_p(__p), _M_log_p(std::log(_M_p))
			{
				_GLIBCXX_DEBUG_ASSERT((_M_p >= 0.0) && (_M_p <= 1.0));
			}
			
			/**
			 * Gets the distribution parameter @p.
			 */
			_RealType
			p() const
			{ return _M_p; }
			
			/**
			 * Gets the inclusive lower bound of the distribution range.
			 */
			result_type
			min() const;

			/**
			 * Gets the inclusive upper bound of the distribution range.
			 */
			result_type
			max() const;

			void
			reset()
			{ }
			
			template<class _UniformRandomNumberGenerator>
				result_type
				operator()(_UniformRandomNumberGenerator& __urng)
				{
					return result_type(std::floor(std::log(_RealType(1.0) - __urng())
																				/ _M_log_p))
						     + result_type(1);
				}

			/**
			 * Inserts a %geometric_distribution random number distribution
			 * @p x into the output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %geometric_distribution random number distribution.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const geometric_distribution& __x)
				{ return __os << __x.p(); }

			/**
			 * Extracts a %geometric_distribution random number distribution
			 * @p u from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __u  A %geometric_distribution random number generator engine.
			 *
			 * @returns The input stream with @p u extracted or in an error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 geometric_distribution& __u)
				{
					__is >> __u._M_p;
					__u._M_log_p = std::log(__u._M_p);
				  return __is;
				}

		protected:
			_RealType _M_p;
			_RealType _M_log_p;
		};

	/* @} */ // group tr1_random_distributions_discrete

	/**
	 * @addtogroup tr1_random_distributions_continuous Continuous Distributions
	 * @ingroup tr1_random_distributions
	 * @{
	 */

	/**
	 * @brief Uniform continuous distribution for random numbers.
	 *
	 * A continuous random distribution on the range [min, max) with equal
	 * probability throughout the range.  The URNG should be real-valued and
	 * deliver number in the range [0, 1).
	 */
	template<typename _RealType = double>
		class uniform_real
		{
		public:
			// types
			typedef _RealType input_type;
			typedef _RealType result_type;
			
		public:
			/**
			 * Constructs a uniform_real object.
			 *
			 * @param min [IN]  The lower bound of the distribution.
			 * @param max [IN]  The upper bound of the distribution.
			 */
			explicit
			uniform_real(_RealType min = _RealType(0), _RealType max = _RealType(1));
			
			result_type
			min() const;
			
			result_type
			max() const;
			
			void reset();
			
			template<class _UniformRandomNumberGenerator>
				result_type
				operator()(_UniformRandomNumberGenerator& __urng)
				{ return (__urng() * (max() - min())) + min(); }

			/**
			 * Inserts a %uniform_real random number distribution @p x into the
			 * output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %uniform_real random number distribution.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const uniform_real& __x)
				{ return __os << __x.min() << " " << __x.max(); }

			/**
			 * Extracts a %unform_real random number distribution
			 * @p u from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __u  A %uniform_real random number generator engine.
			 *
			 * @returns The input stream with @p u extracted or in an error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 uniform_real& __u)
				{ return __is >> __u._M_min >> __u._M_max; }

		protected:
			_RealType _M_min;
			_RealType _M_max;
		};


	/**
	 * @brief An exponential continuous distribution for random numbers.
	 *
	 * The formula for the exponential probability mass function is 
	 * @f$ p(x) = \lambda e^{-\lambda x} @f$.
	 *
	 * <table border=1 cellpadding=10 cellspacing=0>
	 * <caption align=top>Distribution Statistics</caption>
	 * <tr><td>Mean</td><td>@f$ \frac{1}{\lambda} @f$</td></tr>
	 * <tr><td>Median</td><td>@f$ \frac{\ln 2}{\lambda} @f$</td></tr>
	 * <tr><td>Mode</td><td>@f$ zero @f$</td></tr>
	 * <tr><td>Range</td><td>@f$[0, \infty]@f$</td></tr>
	 * <tr><td>Standard Deviation</td><td>@f$ \frac{1}{\lambda} @f$</td></tr>
	 * </table>
	 */
	template<typename _RealType = double>
		class exponential_distribution
		{
		public:
			// types
			typedef _RealType input_type;
			typedef _RealType result_type;
			
		public:
			/**
			 * Constructs an exponential distribution with inverse scale parameter
			 * @f$ \lambda @f$.
			 */
			explicit
			exponential_distribution(const result_type& __lambda = result_type(1))
			: _M_lambda(__lambda)
			{ }
			
			/**
			 * Gets the inverse scale parameter of the distribution.
			 */
			_RealType
			lambda() const
			{ return _M_lambda; }
			
			/**
			 * Resets the distribution.
			 *
			 * Has no effect on exponential distributions.
			 */
			void
			reset()
			{ }
			
			template<class _UniformRandomNumberGenerator>
				result_type
				operator()(_UniformRandomNumberGenerator& __urng)
				{
					return std::log(__urng) / _M_lambda;
				}

			/**
			 * Inserts a %exponential_distribution random number distribution
			 * @p x into the output stream @p os.
			 *
			 * @param __os An output stream.
			 * @param __x  A %exponential_distribution random number distribution.
			 *
			 * @returns The output stream with the state of @p x inserted or in an
			 * error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_ostream<_CharT, _Traits>&
				operator<<(basic_ostream<_CharT, _Traits>& __os,
									 const exponential_distribution& __x)
				{ return __os << __x.lambda(); }

			/**
			 * Extracts a %exponential_distribution random number distribution
			 * @p u from the input stream @p is.
			 *
			 * @param __is An input stream.
			 * @param __u  A %exponential_distribution random number generator engine.
			 *
			 * @returns The input stream with @p u extracted or in an error state.
			 */
			template<typename _CharT, typename _Traits>
				friend basic_istream<_CharT, _Traits>&
				operator>>(basic_istream<_CharT, _Traits>& __is,
									 exponential_distribution& __u)
				{ return __is >> __u._M_lambda; }

		private:
			result_type _M_lambda;
		};

	/* @} */ // group tr1_random_distributions_continuous
	/* @} */ // group tr1_random_distributions
	/* @} */ // group tr1_random

} // namespace tr1
} // namespace std

#include "tr1/random.tcc"

#endif // _STD_TR1_RANDOM
