code:
/* boost random/detail/signed_unsigned_tools.hpp header file
 *
 * Copyright Jens Maurer 2006
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 */

#ifndef BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS
#define BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS

#include <boost/limits.hpp>
#include <boost/config.hpp>
#include <boost/random/traits.hpp>

namespace boost {
namespace random {
namespace detail {

/*
 * Compute x - y, we know that x >= y, return an unsigned value.
 */

template<class T, bool sgn = std::numeric_limits<T>::is_signed && std::numeric_limits<T>::is_bounded>
struct subtract { };

template<class T>
struct subtract<T, /* signed */ false>
{
  typedef T result_type;
  result_type operator()(T x, T y) { return x - y; }
};

template<class T>
struct subtract<T, /* signed */ true>
{
  typedef typename boost::random::traits::make_unsigned_or_unbounded<T>::type result_type;
  result_type operator()(T x, T y)
  {
    if (y >= 0)   // because x >= y, it follows that x >= 0, too
      return result_type(x) - result_type(y);
    if (x >= 0)   // y < 0
      // avoid the nasty two's complement case for y == min()
      return result_type(x) + result_type(-(y+1)) + 1;
    // both x and y are negative: no signed overflow
    return result_type(x - y);
  }
};

/*
 * Compute x + y, x is unsigned, result fits in type of "y".
 */

template<class T1, class T2, bool sgn = (std::numeric_limits<T2>::is_signed && (std::numeric_limits<T1>::digits >= std::numeric_limits<T2>::digits))>
struct add { };

template<class T1, class T2>
struct add<T1, T2, /* signed or else T2 has more digits than T1 so the cast always works - needed when T2 is a multiprecision type and T1 is a native integer */ false>
{
  typedef T2 result_type;
  result_type operator()(T1 x, T2 y) { return T2(x) + y; }
};

template<class T1, class T2>
struct add<T1, T2, /* signed */ true>
{
  typedef T2 result_type;
  result_type operator()(T1 x, T2 y)
  {
    if (y >= 0)
      return T2(x) + y;
    // y < 0
    if (x > T1(-(y+1)))  // result >= 0 after subtraction
      // avoid the nasty two's complement edge case for y == min()
      return T2(x - T1(-(y+1)) - 1);
    // abs(x) < abs(y), thus T2 able to represent x
    return T2(x) + y;
  }
};

} // namespace detail
} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS


code:
/* boost random/detail/uniform_int_float.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 */

#ifndef BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP
#define BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP

#include <boost/limits.hpp>
#include <boost/config.hpp>
#include <boost/integer.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/generator_bits.hpp>

#include <boost/random/detail/disable_warnings.hpp>

namespace boost {
namespace random {
namespace detail {

template<class URNG>
class uniform_int_float
{
public:
    typedef URNG base_type;
    typedef typename base_type::result_type base_result;

    typedef typename boost::uint_t<
        (std::numeric_limits<boost::uintmax_t>::digits <
            std::numeric_limits<base_result>::digits)?
        std::numeric_limits<boost::uintmax_t>::digits :
        std::numeric_limits<base_result>::digits
    >::fast result_type;

    uniform_int_float(base_type& rng)
      : _rng(rng) {}

    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return 0; }
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    {
        std::size_t digits = std::numeric_limits<result_type>::digits;
        if(detail::generator_bits<URNG>::value() < digits) {
            digits = detail::generator_bits<URNG>::value();
        }
        return (result_type(2) << (digits - 1)) - 1;
    }
    base_type& base() { return _rng; }
    const base_type& base() const { return _rng; }

    result_type operator()()
    {
        base_result range = static_cast<base_result>((max)())+1;
        return static_cast<result_type>(_rng() * range);
    }

private:
    base_type& _rng;
};

} // namespace detail
} // namespace random
} // namespace boost

#include <boost/random/detail/enable_warnings.hpp>

#endif // BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP

license: mit


code:
/* boost random/vector_io.hpp header file
 *
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_DETAIL_VECTOR_IO_HPP
#define BOOST_RANDOM_DETAIL_VECTOR_IO_HPP

#include <vector>
#include <iosfwd>
#include <istream>

namespace boost {
namespace random {
namespace detail {

template<class CharT, class Traits, class T>
void print_vector(std::basic_ostream<CharT, Traits>& os,
                  const std::vector<T>& vec)
{
    typename std::vector<T>::const_iterator
        iter = vec.begin(),
        end =  vec.end();
    os << os.widen('[');
    if(iter != end) {
        os << *iter;
        ++iter;
        for(; iter != end; ++iter)
        {
            os << os.widen(' ') << *iter;
        }
    }
    os << os.widen(']');
}

template<class CharT, class Traits, class T>
void read_vector(std::basic_istream<CharT, Traits>& is, std::vector<T>& vec)
{
    CharT ch;
    if(!(is >> ch)) {
        return;
    }
    if(ch != is.widen('[')) {
        is.putback(ch);
        is.setstate(std::ios_base::failbit);
        return;
    }
    T val;
    while(is >> std::ws >> val) {
        vec.push_back(val);
    }
    if(is.fail()) {
        is.clear();
        if(!(is >> ch)) {
            return;
        }
        if(ch != is.widen(']')) {
            is.putback(ch);
            is.setstate(std::ios_base::failbit);
        }
    }
}

}
}
}

#endif // BOOST_RANDOM_DETAIL_VECTOR_IO_HPP


code:
/* boost random/discard_block.hpp header file
 *
 * Copyright Jens Maurer 2002
 * Copyright Steven Watanabe 2010
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 * Revision history
 *  2001-03-02  created
 */

#ifndef BOOST_RANDOM_DISCARD_BLOCK_HPP
#define BOOST_RANDOM_DISCARD_BLOCK_HPP

#include <iostream>
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/seed_impl.hpp>

namespace boost {
namespace random {

/**
 * The class template \discard_block_engine is a model of
 * \pseudo_random_number_generator.  It modifies
 * another generator by discarding parts of its output.
 * Out of every block of @c p results, the first @c r
 * will be returned and the rest discarded.
 *
 * Requires: 0 < p <= r
 */
template<class UniformRandomNumberGenerator, std::size_t p, std::size_t r>
class discard_block_engine
{
    typedef typename detail::seed_type<
        typename UniformRandomNumberGenerator::result_type>::type seed_type;
public:
    typedef UniformRandomNumberGenerator base_type;
    typedef typename base_type::result_type result_type;

    BOOST_STATIC_CONSTANT(std::size_t, block_size = p);
    BOOST_STATIC_CONSTANT(std::size_t, used_block = r);

    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    BOOST_STATIC_CONSTANT(std::size_t, total_block = p);
    BOOST_STATIC_CONSTANT(std::size_t, returned_block = r);

    BOOST_STATIC_ASSERT(total_block >= returned_block);

    /** Uses the default seed for the base generator. */
    discard_block_engine() : _rng(), _n(0) { }
    /** Constructs a new \discard_block_engine with a copy of rng. */
    explicit discard_block_engine(const base_type & rng) : _rng(rng), _n(0) { }

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
    /** Constructs a new \discard_block_engine with rng. */
    explicit discard_block_engine(base_type && rng) : _rng(rng), _n(0) { }
#endif

    /**
     * Creates a new \discard_block_engine and seeds the underlying
     * generator with @c value
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(discard_block_engine,
                                               seed_type, value)
    { _rng.seed(value); _n = 0; }

    /**
     * Creates a new \discard_block_engine and seeds the underlying
     * generator with @c seq
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(discard_block_engine, SeedSeq, seq)
    { _rng.seed(seq); _n = 0; }

    /**
     * Creates a new \discard_block_engine and seeds the underlying
     * generator with first and last.
     */
    template<class It> discard_block_engine(It& first, It last)
      : _rng(first, last), _n(0) { }

    /** default seeds the underlying generator. */
    void seed() { _rng.seed(); _n = 0; }
    /** Seeds the underlying generator with s. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(discard_block_engine, seed_type, s)
    { _rng.seed(s); _n = 0; }
    /** Seeds the underlying generator with seq. */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(discard_block_engine, SeedSeq, seq)
    { _rng.seed(seq); _n = 0; }
    /** Seeds the underlying generator with first and last. */
    template<class It> void seed(It& first, It last)
    { _rng.seed(first, last); _n = 0; }

    /** Returns the underlying engine. */
    const base_type& base() const { return _rng; }

    /** Returns the next value of the generator. */
    result_type operator()()
    {
        if(_n >= returned_block) {
            // discard values of random number generator
            // Don't use discard, since we still need to
            // be somewhat compatible with TR1.
            // _rng.discard(total_block - _n);
            for(std::size_t i = 0; i < total_block - _n; ++i) {
                _rng();
            }
            _n = 0;
        }
        ++_n;
        return _rng();
    }

    void discard(boost::uintmax_t z)
    {
        for(boost::uintmax_t j = 0; j < z; ++j) {
            (*this)();
        }
    }

    template<class It>
    void generate(It first, It last)
    { detail::generate(*this, first, last); }

    /**
     * Returns the smallest value that the generator can produce.
     * This is the same as the minimum of the underlying generator.
     */
    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (base_type::min)(); }
    /**
     * Returns the largest value that the generator can produce.
     * This is the same as the maximum of the underlying generator.
     */
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (base_type::max)(); }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
    /** Writes a \discard_block_engine to a @c std::ostream. */
    template<class CharT, class Traits>
    friend std::basic_ostream<CharT,Traits>&
    operator<<(std::basic_ostream<CharT,Traits>& os,
               const discard_block_engine& s)
    {
        os << s._rng << ' ' << s._n;
        return os;
    }

    /** Reads a \discard_block_engine from a @c std::istream. */
    template<class CharT, class Traits>
    friend std::basic_istream<CharT,Traits>&
    operator>>(std::basic_istream<CharT,Traits>& is, discard_block_engine& s)
    {
        is >> s._rng >> std::ws >> s._n;
        return is;
    }
#endif

    /** Returns true if the two generators will produce identical sequences. */
    friend bool operator==(const discard_block_engine& x,
                           const discard_block_engine& y)
    { return x._rng == y._rng && x._n == y._n; }
    /** Returns true if the two generators will produce different sequences. */
    friend bool operator!=(const discard_block_engine& x,
                           const discard_block_engine& y)
    { return !(x == y); }

private:
    base_type _rng;
    std::size_t _n;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class URNG, std::size_t p, std::size_t r>
const bool discard_block_engine<URNG, p, r>::has_fixed_range;
template<class URNG, std::size_t p, std::size_t r>
const std::size_t discard_block_engine<URNG, p, r>::total_block;
template<class URNG, std::size_t p, std::size_t r>
const std::size_t discard_block_engine<URNG, p, r>::returned_block;
template<class URNG, std::size_t p, std::size_t r>
const std::size_t discard_block_engine<URNG, p, r>::block_size;
template<class URNG, std::size_t p, std::size_t r>
const std::size_t discard_block_engine<URNG, p, r>::used_block;
#endif

/// \cond \show_deprecated

template<class URNG, int p, int r>
class discard_block : public discard_block_engine<URNG, p, r>
{
    typedef discard_block_engine<URNG, p, r> base_t;
public:
    typedef typename base_t::result_type result_type;
    discard_block() {}
    template<class T>
    discard_block(T& arg) : base_t(arg) {}
    template<class T>
    discard_block(const T& arg) : base_t(arg) {}
    template<class It>
    discard_block(It& first, It last) : base_t(first, last) {}
    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (this->base().min)(); }
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (this->base().max)(); }
};

/// \endcond

namespace detail {

    template<class Engine>
    struct generator_bits;

    template<class URNG, std::size_t p, std::size_t r>
    struct generator_bits<discard_block_engine<URNG, p, r> > {
        static std::size_t value() { return generator_bits<URNG>::value(); }
    };

    template<class URNG, int p, int r>
    struct generator_bits<discard_block<URNG, p, r> > {
        static std::size_t value() { return generator_bits<URNG>::value(); }
    };

}

} // namespace random

} // namespace boost

#endif // BOOST_RANDOM_DISCARD_BLOCK_HPP


code:
/* boost random/discrete_distribution.hpp header file
 *
 * Copyright Steven Watanabe 2009-2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED

#include <vector>
#include <limits>
#include <numeric>
#include <utility>
#include <iterator>
#include <boost/assert.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/vector_io.hpp>

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
#include <initializer_list>
#endif

#include <boost/range/begin.hpp>
#include <boost/range/end.hpp>

#include <boost/random/detail/disable_warnings.hpp>

namespace boost {
namespace random {
namespace detail {

template<class IntType, class WeightType>
struct integer_alias_table {
    WeightType get_weight(IntType bin) const {
        WeightType result = _average;
        if(bin < _excess) ++result;
        return result;
    }
    template<class Iter>
    WeightType init_average(Iter begin, Iter end) {
        WeightType weight_average = 0;
        IntType excess = 0;
        IntType n = 0;
        // weight_average * n + excess == current partial sum
        // This is a bit messy, but it's guaranteed not to overflow
        for(Iter iter = begin; iter != end; ++iter) {
            ++n;
            if(*iter < weight_average) {
                WeightType diff = weight_average - *iter;
                weight_average -= diff / n;
                if(diff % n > excess) {
                    --weight_average;
                    excess += n - diff % n;
                } else {
                    excess -= diff % n;
                }
            } else {
                WeightType diff = *iter - weight_average;
                weight_average += diff / n;
                if(diff % n < n - excess) {
                    excess += diff % n;
                } else {
                    ++weight_average;
                    excess -= n - diff % n;
                }
            }
        }
        _alias_table.resize(static_cast<std::size_t>(n));
        _average = weight_average;
        _excess = excess;
        return weight_average;
    }
    void init_empty()
    {
        _alias_table.clear();
        _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
                                              static_cast<IntType>(0)));
        _average = static_cast<WeightType>(1);
        _excess = static_cast<IntType>(0);
    }
    bool operator==(const integer_alias_table& other) const
    {
        return _alias_table == other._alias_table &&
            _average == other._average && _excess == other._excess;
    }
    static WeightType normalize(WeightType val, WeightType average)
    {
        return val;
    }
    static void normalize(std::vector<WeightType>&) {}
    template<class URNG>
    WeightType test(URNG &urng) const
    {
        return uniform_int_distribution<WeightType>(0, _average)(urng);
    }
    bool accept(IntType result, WeightType val) const
    {
        return result < _excess || val < _average;
    }
    static WeightType try_get_sum(const std::vector<WeightType>& weights)
    {
        WeightType result = static_cast<WeightType>(0);
        for(typename std::vector<WeightType>::const_iterator
                iter = weights.begin(), end = weights.end();
            iter != end; ++iter)
        {
            if((std::numeric_limits<WeightType>::max)() - result > *iter) {
                return static_cast<WeightType>(0);
            }
            result += *iter;
        }
        return result;
    }
    template<class URNG>
    static WeightType generate_in_range(URNG &urng, WeightType max)
    {
        return uniform_int_distribution<WeightType>(
            static_cast<WeightType>(0), max-1)(urng);
    }
    typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
    alias_table_t _alias_table;
    WeightType _average;
    IntType _excess;
};

template<class IntType, class WeightType>
struct real_alias_table {
    WeightType get_weight(IntType) const
    {
        return WeightType(1.0);
    }
    template<class Iter>
    WeightType init_average(Iter first, Iter last)
    {
        std::size_t size = std::distance(first, last);
        WeightType weight_sum =
            std::accumulate(first, last, static_cast<WeightType>(0));
        _alias_table.resize(size);
        return weight_sum / size;
    }
    void init_empty()
    {
        _alias_table.clear();
        _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
                                              static_cast<IntType>(0)));
    }
    bool operator==(const real_alias_table& other) const
    {
        return _alias_table == other._alias_table;
    }
    static WeightType normalize(WeightType val, WeightType average)
    {
        return val / average;
    }
    static void normalize(std::vector<WeightType>& weights)
    {
        WeightType sum =
            std::accumulate(weights.begin(), weights.end(),
                            static_cast<WeightType>(0));
        for(typename std::vector<WeightType>::iterator
                iter = weights.begin(),
                end = weights.end();
            iter != end; ++iter)
        {
            *iter /= sum;
        }
    }
    template<class URNG>
    WeightType test(URNG &urng) const
    {
        return uniform_01<WeightType>()(urng);
    }
    bool accept(IntType, WeightType) const
    {
        return true;
    }
    static WeightType try_get_sum(const std::vector<WeightType>& weights)
    {
        return static_cast<WeightType>(1);
    }
    template<class URNG>
    static WeightType generate_in_range(URNG &urng, WeightType)
    {
        return uniform_01<WeightType>()(urng);
    }
    typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
    alias_table_t _alias_table;
};

template<bool IsIntegral>
struct select_alias_table;

template<>
struct select_alias_table<true> {
    template<class IntType, class WeightType>
    struct apply {
        typedef integer_alias_table<IntType, WeightType> type;
    };
};

template<>
struct select_alias_table<false> {
    template<class IntType, class WeightType>
    struct apply {
        typedef real_alias_table<IntType, WeightType> type;
    };
};

}

/**
 * The class @c discrete_distribution models a \random_distribution.
 * It produces integers in the range [0, n) with the probability
 * of producing each value is specified by the parameters of the
 * distribution.
 */
template<class IntType = int, class WeightType = double>
class discrete_distribution {
public:
    typedef WeightType input_type;
    typedef IntType result_type;

    class param_type {
    public:

        typedef discrete_distribution distribution_type;

        /**
         * Constructs a @c param_type object, representing a distribution
         * with \f$p(0) = 1\f$ and \f$p(k|k>0) = 0\f$.
         */
        param_type() : _probabilities(1, static_cast<WeightType>(1)) {}
        /**
         * If @c first == @c last, equivalent to the default constructor.
         * Otherwise, the values of the range represent weights for the
         * possible values of the distribution.
         */
        template<class Iter>
        param_type(Iter first, Iter last) : _probabilities(first, last)
        {
            normalize();
        }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
        /**
         * If wl.size() == 0, equivalent to the default constructor.
         * Otherwise, the values of the @c initializer_list represent
         * weights for the possible values of the distribution.
         */
        param_type(const std::initializer_list<WeightType>& wl)
          : _probabilities(wl)
        {
            normalize();
        }
#endif
        /**
         * If the range is empty, equivalent to the default constructor.
         * Otherwise, the elements of the range represent
         * weights for the possible values of the distribution.
         */
        template<class Range>
        explicit param_type(const Range& range)
          : _probabilities(boost::begin(range), boost::end(range))
        {
            normalize();
        }

        /**
         * If nw is zero, equivalent to the default constructor.
         * Otherwise, the range of the distribution is [0, nw),
         * and the weights are found by  calling fw with values
         * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
         * \f$\mbox{xmax} - \delta/2\f$, where
         * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
         */
        template<class Func>
        param_type(std::size_t nw, double xmin, double xmax, Func fw)
        {
            std::size_t n = (nw == 0) ? 1 : nw;
            double delta = (xmax - xmin) / n;
            BOOST_ASSERT(delta > 0);
            for(std::size_t k = 0; k < n; ++k) {
                _probabilities.push_back(fw(xmin + k*delta + delta/2));
            }
            normalize();
        }

        /**
         * Returns a vector containing the probabilities of each possible
         * value of the distribution.
         */
        std::vector<WeightType> probabilities() const
        {
            return _probabilities;
        }

        /** Writes the parameters to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            detail::print_vector(os, parm._probabilities);
            return os;
        }

        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            std::vector<WeightType> temp;
            detail::read_vector(is, temp);
            if(is) {
                parm._probabilities.swap(temp);
            }
            return is;
        }

        /** Returns true if the two sets of parameters are the same. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        {
            return lhs._probabilities == rhs._probabilities;
        }
        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
    private:
        /// @cond show_private
        friend class discrete_distribution;
        explicit param_type(const discrete_distribution& dist)
          : _probabilities(dist.probabilities())
        {}
        void normalize()
        {
            impl_type::normalize(_probabilities);
        }
        std::vector<WeightType> _probabilities;
        /// @endcond
    };

    /**
     * Creates a new @c discrete_distribution object that has
     * \f$p(0) = 1\f$ and \f$p(i|i>0) = 0\f$.
     */
    discrete_distribution()
    {
        _impl.init_empty();
    }
    /**
     * Constructs a discrete_distribution from an iterator range.
     * If @c first == @c last, equivalent to the default constructor.
     * Otherwise, the values of the range represent weights for the
     * possible values of the distribution.
     */
    template<class Iter>
    discrete_distribution(Iter first, Iter last)
    {
        init(first, last);
    }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    /**
     * Constructs a @c discrete_distribution from a @c std::initializer_list.
     * If the @c initializer_list is empty, equivalent to the default
     * constructor.  Otherwise, the values of the @c initializer_list
     * represent weights for the possible values of the distribution.
     * For example, given the distribution
     *
     * @code
     * discrete_distribution<> dist{1, 4, 5};
     * @endcode
     *
     * The probability of a 0 is 1/10, the probability of a 1 is 2/5,
     * the probability of a 2 is 1/2, and no other values are possible.
     */
    discrete_distribution(std::initializer_list<WeightType> wl)
    {
        init(wl.begin(), wl.end());
    }
#endif
    /**
     * Constructs a discrete_distribution from a Boost.Range range.
     * If the range is empty, equivalent to the default constructor.
     * Otherwise, the values of the range represent weights for the
     * possible values of the distribution.
     */
    template<class Range>
    explicit discrete_distribution(const Range& range)
    {
        init(boost::begin(range), boost::end(range));
    }
    /**
     * Constructs a discrete_distribution that approximates a function.
     * If nw is zero, equivalent to the default constructor.
     * Otherwise, the range of the distribution is [0, nw),
     * and the weights are found by  calling fw with values
     * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
     * \f$\mbox{xmax} - \delta/2\f$, where
     * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
     */
    template<class Func>
    discrete_distribution(std::size_t nw, double xmin, double xmax, Func fw)
    {
        std::size_t n = (nw == 0) ? 1 : nw;
        double delta = (xmax - xmin) / n;
        BOOST_ASSERT(delta > 0);
        std::vector<WeightType> weights;
        for(std::size_t k = 0; k < n; ++k) {
            weights.push_back(fw(xmin + k*delta + delta/2));
        }
        init(weights.begin(), weights.end());
    }
    /**
     * Constructs a discrete_distribution from its parameters.
     */
    explicit discrete_distribution(const param_type& parm)
    {
        param(parm);
    }

    /**
     * Returns a value distributed according to the parameters of the
     * discrete_distribution.
     */
    template<class URNG>
    IntType operator()(URNG& urng) const
    {
        BOOST_ASSERT(!_impl._alias_table.empty());
        IntType result;
        WeightType test;
        do {
            result = uniform_int_distribution<IntType>((min)(), (max)())(urng);
            test = _impl.test(urng);
        } while(!_impl.accept(result, test));
        if(test < _impl._alias_table[static_cast<std::size_t>(result)].first) {
            return result;
        } else {
            return(_impl._alias_table[static_cast<std::size_t>(result)].second);
        }
    }

    /**
     * Returns a value distributed according to the parameters
     * specified by param.
     */
    template<class URNG>
    IntType operator()(URNG& urng, const param_type& parm) const
    {
        if(WeightType limit = impl_type::try_get_sum(parm._probabilities)) {
            WeightType val = impl_type::generate_in_range(urng, limit);
            WeightType sum = 0;
            std::size_t result = 0;
            for(typename std::vector<WeightType>::const_iterator
                    iter = parm._probabilities.begin(),
                    end = parm._probabilities.end();
                iter != end; ++iter, ++result)
            {
                sum += *iter;
                if(sum > val) {
                    return result;
                }
            }
            // This shouldn't be reachable, but round-off error
            // can prevent any match from being found when val is
            // very close to 1.
            return static_cast<IntType>(parm._probabilities.size() - 1);
        } else {
            // WeightType is integral and sum(parm._probabilities)
            // would overflow.  Just use the easy solution.
            return discrete_distribution(parm)(urng);
        }
    }

    /** Returns the smallest value that the distribution can produce. */
    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
    /** Returns the largest value that the distribution can produce. */
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return static_cast<result_type>(_impl._alias_table.size() - 1); }

    /**
     * Returns a vector containing the probabilities of each
     * value of the distribution.  For example, given
     *
     * @code
     * discrete_distribution<> dist = { 1, 4, 5 };
     * std::vector<double> p = dist.param();
     * @endcode
     *
     * the vector, p will contain {0.1, 0.4, 0.5}.
     *
     * If @c WeightType is integral, then the weights
     * will be returned unchanged.
     */
    std::vector<WeightType> probabilities() const
    {
        std::vector<WeightType> result(_impl._alias_table.size(), static_cast<WeightType>(0));
        std::size_t i = 0;
        for(typename impl_type::alias_table_t::const_iterator
                iter = _impl._alias_table.begin(),
                end = _impl._alias_table.end();
                iter != end; ++iter, ++i)
        {
            WeightType val = iter->first;
            result[i] += val;
            result[static_cast<std::size_t>(iter->second)] += _impl.get_weight(i) - val;
        }
        impl_type::normalize(result);
        return(result);
    }

    /** Returns the parameters of the distribution. */
    param_type param() const
    {
        return param_type(*this);
    }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        init(parm._probabilities.begin(), parm._probabilities.end());
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() {}

    /** Writes a distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, discrete_distribution, dd)
    {
        os << dd.param();
        return os;
    }

    /** Reads a distribution from a @c std::istream */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, discrete_distribution, dd)
    {
        param_type parm;
        if(is >> parm) {
            dd.param(parm);
        }
        return is;
    }

    /**
     * Returns true if the two distributions will return the
     * same sequence of values, when passed equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution, lhs, rhs)
    {
        return lhs._impl == rhs._impl;
    }
    /**
     * Returns true if the two distributions may return different
     * sequences of values, when passed equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(discrete_distribution)

private:

    /// @cond show_private

    template<class Iter>
    void init(Iter first, Iter last, std::input_iterator_tag)
    {
        std::vector<WeightType> temp(first, last);
        init(temp.begin(), temp.end());
    }
    template<class Iter>
    void init(Iter first, Iter last, std::forward_iterator_tag)
    {
        std::vector<std::pair<WeightType, IntType> > below_average;
        std::vector<std::pair<WeightType, IntType> > above_average;
        WeightType weight_average = _impl.init_average(first, last);
        WeightType normalized_average = _impl.get_weight(0);
        std::size_t i = 0;
        for(; first != last; ++first, ++i) {
            WeightType val = impl_type::normalize(*first, weight_average);
            std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i));
            if(val < normalized_average) {
                below_average.push_back(elem);
            } else {
                above_average.push_back(elem);
            }
        }

        typename impl_type::alias_table_t::iterator
            b_iter = below_average.begin(),
            b_end = below_average.end(),
            a_iter = above_average.begin(),
            a_end = above_average.end()
            ;
        while(b_iter != b_end && a_iter != a_end) {
            _impl._alias_table[static_cast<std::size_t>(b_iter->second)] =
                std::make_pair(b_iter->first, a_iter->second);
            a_iter->first -= (_impl.get_weight(b_iter->second) - b_iter->first);
            if(a_iter->first < normalized_average) {
                *b_iter = *a_iter++;
            } else {
                ++b_iter;
            }
        }
        for(; b_iter != b_end; ++b_iter) {
            _impl._alias_table[static_cast<std::size_t>(b_iter->second)].first =
                _impl.get_weight(b_iter->second);
        }
        for(; a_iter != a_end; ++a_iter) {
            _impl._alias_table[static_cast<std::size_t>(a_iter->second)].first =
                _impl.get_weight(a_iter->second);
        }
    }
    template<class Iter>
    void init(Iter first, Iter last)
    {
        if(first == last) {
            _impl.init_empty();
        } else {
            typename std::iterator_traits<Iter>::iterator_category category;
            init(first, last, category);
        }
    }
    typedef typename detail::select_alias_table<
        (::boost::is_integral<WeightType>::value)
    >::template apply<IntType, WeightType>::type impl_type;
    impl_type _impl;
    /// @endcond
};

}
}

#include <boost/random/detail/enable_warnings.hpp>

#endif


code:
/* boost random/exponential_distribution.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <iosfwd>
#include <boost/assert.hpp>
#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/uniform_01.hpp>

namespace boost {
namespace random {

/**
 * The exponential distribution is a model of \random_distribution with
 * a single parameter lambda.
 *
 * It has \f$\displaystyle p(x) = \lambda e^{-\lambda x}\f$
 */
template<class RealType = double>
class exponential_distribution
{
public:
    typedef RealType input_type;
    typedef RealType result_type;

    class param_type
    {
    public:

        typedef exponential_distribution distribution_type;

        /**
         * Constructs parameters with a given lambda.
         *
         * Requires: lambda > 0
         */
        param_type(RealType lambda_arg = RealType(1.0))
          : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }

        /** Returns the lambda parameter of the distribution. */
        RealType lambda() const { return _lambda; }

        /** Writes the parameters to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            os << parm._lambda;
            return os;
        }

        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            is >> parm._lambda;
            return is;
        }

        /** Returns true if the two sets of parameters are equal. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._lambda == rhs._lambda; }

        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _lambda;
    };

    /**
     * Constructs an exponential_distribution with a given lambda.
     *
     * Requires: lambda > 0
     */
    explicit exponential_distribution(RealType lambda_arg = RealType(1.0))
      : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }

    /**
     * Constructs an exponential_distribution from its parameters
     */
    explicit exponential_distribution(const param_type& parm)
      : _lambda(parm.lambda()) {}

    // compiler-generated copy ctor and assignment operator are fine

    /** Returns the lambda parameter of the distribution. */
    RealType lambda() const { return _lambda; }

    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return RealType(0); }
    /** Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return (std::numeric_limits<RealType>::infinity)(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_lambda); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm) { _lambda = parm.lambda(); }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { }

    /**
     * Returns a random variate distributed according to the
     * exponential distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng) const
    {
        using std::log;
        return -result_type(1) /
            _lambda * log(result_type(1)-uniform_01<RealType>()(eng));
    }

    /**
     * Returns a random variate distributed according to the exponential
     * distribution with parameters specified by param.
     */
    template<class Engine>
    result_type operator()(Engine& eng, const param_type& parm) const
    {
        return exponential_distribution(parm)(eng);
    }

    /** Writes the distribution to a std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, exponential_distribution, ed)
    {
        os << ed._lambda;
        return os;
    }

    /** Reads the distribution from a std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, exponential_distribution, ed)
    {
        is >> ed._lambda;
        return is;
    }

    /**
     * Returns true iff the two distributions will produce identical
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(exponential_distribution, lhs, rhs)
    { return lhs._lambda == rhs._lambda; }

    /**
     * Returns true iff the two distributions will produce different
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(exponential_distribution)

private:
    result_type _lambda;
};

} // namespace random

using random::exponential_distribution;

} // namespace boost

#endif // BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP


code:
/* boost random/extreme_value_distribution.hpp header file
 *
 * Copyright Steven Watanabe 2010
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP
#define BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <iosfwd>
#include <istream>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/uniform_01.hpp>

namespace boost {
namespace random {

/**
 * The extreme value distribution is a real valued distribution with two
 * parameters a and b.
 *
 * It has \f$\displaystyle p(x) = \frac{1}{b}e^{\frac{a-x}{b} - e^\frac{a-x}{b}}\f$.
 */
template<class RealType = double>
class extreme_value_distribution {
public:
    typedef RealType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef extreme_value_distribution distribution_type;

        /**
         * Constructs a @c param_type from the "a" and "b" parameters
         * of the distribution.
         *
         * Requires: b > 0
         */
        explicit param_type(RealType a_arg = 1.0, RealType b_arg = 1.0)
          : _a(a_arg), _b(b_arg)
        {}

        /** Returns the "a" parameter of the distribtuion. */
        RealType a() const { return _a; }
        /** Returns the "b" parameter of the distribution. */
        RealType b() const { return _b; }

        /** Writes a @c param_type to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        { os << parm._a << ' ' << parm._b; return os; }

        /** Reads a @c param_type from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        { is >> parm._a >> std::ws >> parm._b; return is; }

        /** Returns true if the two sets of parameters are the same. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._a == rhs._a && lhs._b == rhs._b; }

        /** Returns true if the two sets of parameters are the different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _a;
        RealType _b;
    };

    /**
     * Constructs an @c extreme_value_distribution from its "a" and "b" parameters.
     *
     * Requires: b > 0
     */
    explicit extreme_value_distribution(RealType a_arg = 1.0, RealType b_arg = 1.0)
      : _a(a_arg), _b(b_arg)
    {}
    /** Constructs an @c extreme_value_distribution from its parameters. */
    explicit extreme_value_distribution(const param_type& parm)
      : _a(parm.a()), _b(parm.b())
    {}

    /**
     * Returns a random variate distributed according to the
     * @c extreme_value_distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng) const
    {
        using std::log;
        return _a - log(-log(uniform_01<RealType>()(urng))) * _b;
    }

    /**
     * Returns a random variate distributed accordint to the extreme
     * value distribution with parameters specified by @c param.
     */
    template<class URNG>
    RealType operator()(URNG& urng, const param_type& parm) const
    {
        return extreme_value_distribution(parm)(urng);
    }

    /** Returns the "a" parameter of the distribution. */
    RealType a() const { return _a; }
    /** Returns the "b" parameter of the distribution. */
    RealType b() const { return _b; }

    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return -std::numeric_limits<RealType>::infinity(); }
    /** Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return std::numeric_limits<RealType>::infinity(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_a, _b); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _a = parm.a();
        _b = parm.b();
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { }

    /** Writes an @c extreme_value_distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, extreme_value_distribution, wd)
    {
        os << wd.param();
        return os;
    }

    /** Reads an @c extreme_value_distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, extreme_value_distribution, wd)
    {
        param_type parm;
        if(is >> parm) {
            wd.param(parm);
        }
        return is;
    }

    /**
     * Returns true if the two instances of @c extreme_value_distribution will
     * return identical sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(extreme_value_distribution, lhs, rhs)
    { return lhs._a == rhs._a && lhs._b == rhs._b; }

    /**
     * Returns true if the two instances of @c extreme_value_distribution will
     * return different sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(extreme_value_distribution)

private:
    RealType _a;
    RealType _b;
};

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP


code:
/* boost random/fisher_f_distribution.hpp header file
 *
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_FISHER_F_DISTRIBUTION_HPP
#define BOOST_RANDOM_FISHER_F_DISTRIBUTION_HPP

#include <iosfwd>
#include <istream>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/chi_squared_distribution.hpp>

namespace boost {
namespace random {

/**
 * The Fisher F distribution is a real valued distribution with two
 * parameters m and n.
 *
 * It has \f$\displaystyle p(x) =
 *   \frac{\Gamma((m+n)/2)}{\Gamma(m/2)\Gamma(n/2)}
 *   \left(\frac{m}{n}\right)^{m/2}
 *   x^{(m/2)-1} \left(1+\frac{mx}{n}\right)^{-(m+n)/2}
 * \f$.
 */
template<class RealType = double>
class fisher_f_distribution {
public:
    typedef RealType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef fisher_f_distribution distribution_type;

        /**
         * Constructs a @c param_type from the "m" and "n" parameters
         * of the distribution.
         *
         * Requires: m > 0 and n > 0
         */
        explicit param_type(RealType m_arg = RealType(1.0),
                            RealType n_arg = RealType(1.0))
          : _m(m_arg), _n(n_arg)
        {}

        /** Returns the "m" parameter of the distribtuion. */
        RealType m() const { return _m; }
        /** Returns the "n" parameter of the distribution. */
        RealType n() const { return _n; }

        /** Writes a @c param_type to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        { os << parm._m << ' ' << parm._n; return os; }

        /** Reads a @c param_type from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        { is >> parm._m >> std::ws >> parm._n; return is; }

        /** Returns true if the two sets of parameters are the same. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._m == rhs._m && lhs._n == rhs._n; }

        /** Returns true if the two sets of parameters are the different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _m;
        RealType _n;
    };

    /**
     * Constructs a @c fisher_f_distribution from its "m" and "n" parameters.
     *
     * Requires: m > 0 and n > 0
     */
    explicit fisher_f_distribution(RealType m_arg = RealType(1.0),
                                   RealType n_arg = RealType(1.0))
      : _impl_m(m_arg), _impl_n(n_arg)
    {}
    /** Constructs an @c fisher_f_distribution from its parameters. */
    explicit fisher_f_distribution(const param_type& parm)
      : _impl_m(parm.m()), _impl_n(parm.n())
    {}

    /**
     * Returns a random variate distributed according to the
     * F distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng)
    {
        return (_impl_m(urng) * n()) / (_impl_n(urng) * m());
    }

    /**
     * Returns a random variate distributed according to the
     * F distribution with parameters specified by @c param.
     */
    template<class URNG>
    RealType operator()(URNG& urng, const param_type& parm) const
    {
        return fisher_f_distribution(parm)(urng);
    }

    /** Returns the "m" parameter of the distribution. */
    RealType m() const { return _impl_m.n(); }
    /** Returns the "n" parameter of the distribution. */
    RealType n() const { return _impl_n.n(); }

    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
    /** Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return std::numeric_limits<RealType>::infinity(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(m(), n()); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        typedef chi_squared_distribution<RealType> impl_type;
        typename impl_type::param_type m_param(parm.m());
        _impl_m.param(m_param);
        typename impl_type::param_type n_param(parm.n());
        _impl_n.param(n_param);
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { }

    /** Writes an @c fisher_f_distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, fisher_f_distribution, fd)
    {
        os << fd.param();
        return os;
    }

    /** Reads an @c fisher_f_distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, fisher_f_distribution, fd)
    {
        param_type parm;
        if(is >> parm) {
            fd.param(parm);
        }
        return is;
    }

    /**
     * Returns true if the two instances of @c fisher_f_distribution will
     * return identical sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(fisher_f_distribution, lhs, rhs)
    { return lhs._impl_m == rhs._impl_m && lhs._impl_n == rhs._impl_n; }

    /**
     * Returns true if the two instances of @c fisher_f_distribution will
     * return different sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(fisher_f_distribution)

private:
    chi_squared_distribution<RealType> _impl_m;
    chi_squared_distribution<RealType> _impl_n;
};

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP


code:
/* boost random/gamma_distribution.hpp header file
 *
 * Copyright Jens Maurer 2002
 * Copyright Steven Watanabe 2010
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 */

#ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
#define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <istream>
#include <iosfwd>
#include <boost/assert.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/exponential_distribution.hpp>

namespace boost {
namespace random {

// The algorithm is taken from Knuth

/**
 * The gamma distribution is a continuous distribution with two
 * parameters alpha and beta.  It produces values > 0.
 *
 * It has
 * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
 */
template<class RealType = double>
class gamma_distribution
{
public:
    typedef RealType input_type;
    typedef RealType result_type;

    class param_type
    {
    public:
        typedef gamma_distribution distribution_type;

        /**
         * Constructs a @c param_type object from the "alpha" and "beta"
         * parameters.
         *
         * Requires: alpha > 0 && beta > 0
         */
        param_type(const RealType& alpha_arg = RealType(1.0),
                   const RealType& beta_arg = RealType(1.0))
          : _alpha(alpha_arg), _beta(beta_arg)
        {
        }

        /** Returns the "alpha" parameter of the distribution. */
        RealType alpha() const { return _alpha; }
        /** Returns the "beta" parameter of the distribution. */
        RealType beta() const { return _beta; }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
        /** Writes the parameters to a @c std::ostream. */
        template<class CharT, class Traits>
        friend std::basic_ostream<CharT, Traits>&
        operator<<(std::basic_ostream<CharT, Traits>& os,
                   const param_type& parm)
        {
            os << parm._alpha << ' ' << parm._beta;
            return os;
        }

        /** Reads the parameters from a @c std::istream. */
        template<class CharT, class Traits>
        friend std::basic_istream<CharT, Traits>&
        operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
        {
            is >> parm._alpha >> std::ws >> parm._beta;
            return is;
        }
#endif

        /** Returns true if the two sets of parameters are the same. */
        friend bool operator==(const param_type& lhs, const param_type& rhs)
        {
            return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
        }
        /** Returns true if the two sets fo parameters are different. */
        friend bool operator!=(const param_type& lhs, const param_type& rhs)
        {
            return !(lhs == rhs);
        }
    private:
        RealType _alpha;
        RealType _beta;
    };

#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
    BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
#endif

    /**
     * Creates a new gamma_distribution with parameters "alpha" and "beta".
     *
     * Requires: alpha > 0 && beta > 0
     */
    explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
                                const result_type& beta_arg = result_type(1.0))
      : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
    {
        BOOST_ASSERT(_alpha > result_type(0));
        BOOST_ASSERT(_beta > result_type(0));
        init();
    }

    /** Constructs a @c gamma_distribution from its parameters. */
    explicit gamma_distribution(const param_type& parm)
      : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
    {
        init();
    }

    // compiler-generated copy ctor and assignment operator are fine

    /** Returns the "alpha" paramter of the distribution. */
    RealType alpha() const { return _alpha; }
    /** Returns the "beta" parameter of the distribution. */
    RealType beta() const { return _beta; }
    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
    /* Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return (std::numeric_limits<RealType>::infinity)(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_alpha, _beta); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _alpha = parm.alpha();
        _beta = parm.beta();
        init();
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { _exp.reset(); }

    /**
     * Returns a random variate distributed according to
     * the gamma distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng)
    {
#ifndef BOOST_NO_STDC_NAMESPACE
        // allow for Koenig lookup
        using std::tan; using std::sqrt; using std::exp; using std::log;
        using std::pow;
#endif
        if(_alpha == result_type(1)) {
            return _exp(eng) * _beta;
        } else if(_alpha > result_type(1)) {
            // Can we have a boost::mathconst please?
            const result_type pi = result_type(3.14159265358979323846);
            for(;;) {
                result_type y = tan(pi * uniform_01<RealType>()(eng));
                result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
                    + _alpha-result_type(1);
                if(x <= result_type(0))
                    continue;
                if(uniform_01<RealType>()(eng) >
                    (result_type(1)+y*y) * exp((_alpha-result_type(1))
                                               *log(x/(_alpha-result_type(1)))
                                               - sqrt(result_type(2)*_alpha
                                                      -result_type(1))*y))
                    continue;
                return x * _beta;
            }
        } else /* alpha < 1.0 */ {
            for(;;) {
                result_type u = uniform_01<RealType>()(eng);
                result_type y = _exp(eng);
                result_type x, q;
                if(u < _p) {
                    x = exp(-y/_alpha);
                    q = _p*exp(-x);
                } else {
                    x = result_type(1)+y;
                    q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
                }
                if(u >= q)
                    continue;
                return x * _beta;
            }
        }
    }

    template<class URNG>
    RealType operator()(URNG& urng, const param_type& parm) const
    {
        return gamma_distribution(parm)(urng);
    }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
    /** Writes a @c gamma_distribution to a @c std::ostream. */
    template<class CharT, class Traits>
    friend std::basic_ostream<CharT,Traits>&
    operator<<(std::basic_ostream<CharT,Traits>& os,
               const gamma_distribution& gd)
    {
        os << gd.param();
        return os;
    }

    /** Reads a @c gamma_distribution from a @c std::istream. */
    template<class CharT, class Traits>
    friend std::basic_istream<CharT,Traits>&
    operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
    {
        gd.read(is);
        return is;
    }
#endif

    /**
     * Returns true if the two distributions will produce identical
     * sequences of random variates given equal generators.
     */
    friend bool operator==(const gamma_distribution& lhs,
                           const gamma_distribution& rhs)
    {
        return lhs._alpha == rhs._alpha
            && lhs._beta == rhs._beta
            && lhs._exp == rhs._exp;
    }

    /**
     * Returns true if the two distributions can produce different
     * sequences of random variates, given equal generators.
     */
    friend bool operator!=(const gamma_distribution& lhs,
                           const gamma_distribution& rhs)
    {
        return !(lhs == rhs);
    }

private:
    /// \cond hide_private_members

    template<class CharT, class Traits>
    void read(std::basic_istream<CharT, Traits>& is)
    {
        param_type parm;
        if(is >> parm) {
            param(parm);
        }
    }

    void init()
    {
#ifndef BOOST_NO_STDC_NAMESPACE
        // allow for Koenig lookup
        using std::exp;
#endif
        _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
    }
    /// \endcond

    exponential_distribution<RealType> _exp;
    result_type _alpha;
    result_type _beta;
    // some data precomputed from the parameters
    result_type _p;
};

} // namespace random

using random::gamma_distribution;

} // namespace boost

#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP


code:
/* boost random/generate_canonical.hpp header file
 *
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 */

#ifndef BOOST_RANDOM_GENERATE_CANONICAL_HPP
#define BOOST_RANDOM_GENERATE_CANONICAL_HPP

#include <algorithm>
#include <boost/assert.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/limits.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/random/detail/generator_bits.hpp>

namespace boost {
namespace random {

namespace detail {

template<class RealType, std::size_t bits, class URNG>
RealType generate_canonical_impl(URNG& g, boost::mpl::true_ /*is_integral*/)
{
    using std::pow;
    typedef typename URNG::result_type base_result;
    std::size_t digits = std::numeric_limits<RealType>::digits;
    RealType R = RealType((g.max)()) - RealType((g.min)()) + 1;
    RealType mult = R;
    RealType limit =
        pow(RealType(2),
            RealType((std::min)(static_cast<std::size_t>(bits), digits)));
    RealType S = RealType(detail::subtract<base_result>()(g(), (g.min)()));
    while(mult < limit) {
        RealType inc = RealType(detail::subtract<base_result>()(g(), (g.min)()));
        S += inc * mult;
        mult *= R;
    }
    return S / mult;
}

template<class RealType, std::size_t bits, class URNG>
RealType generate_canonical_impl(URNG& g, boost::mpl::false_ /*is_integral*/)
{
    using std::pow;
    using std::floor;
    BOOST_ASSERT((g.min)() == 0);
    BOOST_ASSERT((g.max)() == 1);
    std::size_t digits = std::numeric_limits<RealType>::digits;
    std::size_t engine_bits = detail::generator_bits<URNG>::value();
    std::size_t b = (std::min)(bits, digits);
    RealType R = pow(RealType(2), RealType(engine_bits));
    RealType mult = R;
    RealType limit = pow(RealType(2), RealType(b));
    RealType S = RealType(g() - (g.min)());
    while(mult < limit) {
        RealType inc(floor((RealType(g()) - RealType((g.min)())) * R));
        S += inc * mult;
        mult *= R;
    }
    return S / mult;
}

}

/**
 * Returns a value uniformly distributed in the range [0, 1)
 * with at least @c bits random bits.
 */
template<class RealType, std::size_t bits, class URNG>
RealType generate_canonical(URNG& g)
{
    RealType result = detail::generate_canonical_impl<RealType, bits>(
        g, boost::random::traits::is_integral<typename URNG::result_type>());
    BOOST_ASSERT(result >= 0);
    BOOST_ASSERT(result <= 1);
    if(result == 1) {
        result -= std::numeric_limits<RealType>::epsilon() / 2;
        BOOST_ASSERT(result != 1);
    }
    return result;
}

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_GENERATE_CANONICAL_HPP


code:
/* boost random/geometric_distribution.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP
#define BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>          // std::log
#include <iosfwd>
#include <ios>
#include <boost/assert.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/uniform_01.hpp>

namespace boost {
namespace random {

/**
 * An instantiation of the class template @c geometric_distribution models
 * a \random_distribution.  The distribution produces positive
 * integers which are the number of bernoulli trials
 * with probability @c p required to get one that fails.
 *
 * For the geometric distribution, \f$p(i) = p(1-p)^{i}\f$.
 *
 * @xmlwarning
 * This distribution has been updated to match the C++ standard.
 * Its behavior has changed from the original
 * boost::geometric_distribution.  A backwards compatible
 * wrapper is provided in namespace boost.
 * @endxmlwarning
 */
template<class IntType = int, class RealType = double>
class geometric_distribution
{
public:
    typedef RealType input_type;
    typedef IntType result_type;

    class param_type
    {
    public:

        typedef geometric_distribution distribution_type;

        /** Constructs the parameters with p. */
        explicit param_type(RealType p_arg = RealType(0.5))
          : _p(p_arg)
        {
            BOOST_ASSERT(RealType(0) < _p && _p < RealType(1));
        }

        /** Returns the p parameter of the distribution. */
        RealType p() const { return _p; }

        /** Writes the parameters to a std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            os << parm._p;
            return os;
        }

        /** Reads the parameters from a std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            double p_in;
            if(is >> p_in) {
                if(p_in > RealType(0) && p_in < RealType(1)) {
                    parm._p = p_in;
                } else {
                    is.setstate(std::ios_base::failbit);
                }
            }
            return is;
        }

        /** Returns true if the two sets of parameters are equal. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._p == rhs._p; }

        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _p;
    };

    /**
     * Contructs a new geometric_distribution with the paramter @c p.
     *
     * Requires: 0 < p < 1
     */
    explicit geometric_distribution(const RealType& p_arg = RealType(0.5))
      : _p(p_arg)
    {
        BOOST_ASSERT(RealType(0) < _p && _p < RealType(1));
        init();
    }

    /** Constructs a new geometric_distribution from its parameters. */
    explicit geometric_distribution(const param_type& parm)
      : _p(parm.p())
    {
        init();
    }

    // compiler-generated copy ctor and assignment operator are fine

    /** Returns: the distribution parameter @c p  */
    RealType p() const { return _p; }

    /** Returns the smallest value that the distribution can produce. */
    IntType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return IntType(0); }

    /** Returns the largest value that the distribution can produce. */
    IntType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return (std::numeric_limits<IntType>::max)(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_p); }

    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _p = parm.p();
        init();
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { }

    /**
     * Returns a random variate distributed according to the
     * geometric_distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng) const
    {
        using std::log;
        using std::floor;
        RealType x = RealType(1) - boost::uniform_01<RealType>()(eng);
        return IntType(floor(log(x) / _log_1mp));
    }

    /**
     * Returns a random variate distributed according to the
     * geometric distribution with parameters specified by param.
     */
    template<class Engine>
    result_type operator()(Engine& eng, const param_type& parm) const
    { return geometric_distribution(parm)(eng); }

    /** Writes the distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, geometric_distribution, gd)
    {
        os << gd._p;
        return os;
    }

    /** Reads the distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, geometric_distribution, gd)
    {
        param_type parm;
        if(is >> parm) {
            gd.param(parm);
        }
        return is;
    }

    /**
     * Returns true if the two distributions will produce identical
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(geometric_distribution, lhs, rhs)
    { return lhs._p == rhs._p; }

    /**
     * Returns true if the two distributions may produce different
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(geometric_distribution)

private:

    /// \cond show_private

    void init()
    {
        using std::log;
        _log_1mp = log(1 - _p);
    }

    RealType _p;
    RealType _log_1mp;

    /// \endcond
};

} // namespace random

/// \cond show_deprecated

/**
 * Provided for backwards compatibility.  This class is
 * deprecated.  It provides the old behavior of geometric_distribution
 * with \f$p(i) = (1-p) p^{i-1}\f$.
 */
template<class IntType = int, class RealType = double>
class geometric_distribution
{
public:
    typedef RealType input_type;
    typedef IntType result_type;

    explicit geometric_distribution(RealType p_arg = RealType(0.5))
      : _impl(1 - p_arg) {}

    RealType p() const { return 1 - _impl.p(); }

    void reset() {}

    template<class Engine>
    IntType operator()(Engine& eng) const { return _impl(eng) + IntType(1); }

    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, geometric_distribution, gd)
    {
        os << gd.p();
        return os;
    }

    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, geometric_distribution, gd)
    {
        RealType val;
        if(is >> val) {
            typename impl_type::param_type impl_param(1 - val);
            gd._impl.param(impl_param);
        }
        return is;
    }

private:
    typedef random::geometric_distribution<IntType, RealType> impl_type;
    impl_type _impl;
};

/// \endcond

} // namespace boost

#endif // BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP