code:
/* boost random/piecewise_linear_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_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED

#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <boost/assert.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/discrete_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>

namespace boost {
namespace random {

/**
 * The class @c piecewise_linear_distribution models a \random_distribution.
 */
template<class RealType = double>
class piecewise_linear_distribution {
public:
    typedef std::size_t input_type;
    typedef RealType result_type;

    class param_type {
    public:

        typedef piecewise_linear_distribution distribution_type;

        /**
         * Constructs a @c param_type object, representing a distribution
         * that produces values uniformly distributed in the range [0, 1).
         */
        param_type()
        {
            _weights.push_back(RealType(1));
            _weights.push_back(RealType(1));
            _intervals.push_back(RealType(0));
            _intervals.push_back(RealType(1));
        }
        /**
         * Constructs a @c param_type object from two iterator ranges
         * containing the interval boundaries and weights at the boundaries.
         * If there are fewer than two boundaries, then this is equivalent to
         * the default constructor and the distribution will produce values
         * uniformly distributed in the range [0, 1).
         *
         * The values of the interval boundaries must be strictly
         * increasing, and the number of weights must be the same as
         * the number of interval boundaries.  If there are extra
         * weights, they are ignored.
         */
        template<class IntervalIter, class WeightIter>
        param_type(IntervalIter intervals_first, IntervalIter intervals_last,
                   WeightIter weight_first)
          : _intervals(intervals_first, intervals_last)
        {
            if(_intervals.size() < 2) {
                _intervals.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            } else {
                _weights.reserve(_intervals.size());
                for(std::size_t i = 0; i < _intervals.size(); ++i) {
                    _weights.push_back(*weight_first++);
                }
            }
        }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
        /**
         * Constructs a @c param_type object from an initializer_list
         * containing the interval boundaries and a unary function
         * specifying the weights at the boundaries.  Each weight is
         * determined by calling the function at the corresponding point.
         *
         * If the initializer_list contains fewer than two elements,
         * this is equivalent to the default constructor and the
         * distribution will produce values uniformly distributed
         * in the range [0, 1).
         */
        template<class T, class F>
        param_type(const std::initializer_list<T>& il, F f)
          : _intervals(il.begin(), il.end())
        {
            if(_intervals.size() < 2) {
                _intervals.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            } else {
                _weights.reserve(_intervals.size());
                for(typename std::vector<RealType>::const_iterator
                    iter = _intervals.begin(), end = _intervals.end();
                    iter != end; ++iter)
                {
                    _weights.push_back(f(*iter));
                }
            }
        }
#endif
        /**
         * Constructs a @c param_type object from Boost.Range ranges holding
         * the interval boundaries and the weights at the boundaries.  If
         * there are fewer than two interval boundaries, this is equivalent
         * to the default constructor and the distribution will produce
         * values uniformly distributed in the range [0, 1).  The
         * number of weights must be equal to the number of
         * interval boundaries.
         */
        template<class IntervalRange, class WeightRange>
        param_type(const IntervalRange& intervals_arg,
                   const WeightRange& weights_arg)
          : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
            _weights(boost::begin(weights_arg), boost::end(weights_arg))
        {
            if(_intervals.size() < 2) {
                _weights.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.clear();
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            }
        }

        /**
         * Constructs the parameters for a distribution that approximates a
         * function.  The range of the distribution is [xmin, xmax).  This
         * range is divided into nw equally sized intervals and the weights
         * are found by calling the unary function f on the boundaries of the
         * intervals.
         */
        template<class F>
        param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
        {
            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) {
                _weights.push_back(f(xmin + k*delta));
                _intervals.push_back(xmin + k*delta);
            }
            _weights.push_back(f(xmax));
            _intervals.push_back(xmax);
        }

        /**  Returns a vector containing the interval boundaries. */
        std::vector<RealType> intervals() const { return _intervals; }

        /**
         * Returns a vector containing the probability densities
         * at all the interval boundaries.
         */
        std::vector<RealType> densities() const
        {
            RealType sum = static_cast<RealType>(0);
            for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
                RealType width = _intervals[i + 1] - _intervals[i];
                sum += (_weights[i] + _weights[i + 1]) * width / 2;
            }
            std::vector<RealType> result;
            result.reserve(_weights.size());
            for(typename std::vector<RealType>::const_iterator
                iter = _weights.begin(), end = _weights.end();
                iter != end; ++iter)
            {
                result.push_back(*iter / sum);
            }
            return result;
        }

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

        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            std::vector<RealType> new_intervals;
            std::vector<RealType> new_weights;
            detail::read_vector(is, new_intervals);
            detail::read_vector(is, new_weights);
            if(is) {
                parm._intervals.swap(new_intervals);
                parm._weights.swap(new_weights);
            }
            return is;
        }

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

    private:
        friend class piecewise_linear_distribution;

        std::vector<RealType> _intervals;
        std::vector<RealType> _weights;
    };

    /**
     * Creates a new @c piecewise_linear_distribution that
     * produces values uniformly distributed in the range [0, 1).
     */
    piecewise_linear_distribution()
    {
        default_init();
    }
    /**
     * Constructs a piecewise_linear_distribution from two iterator ranges
     * containing the interval boundaries and the weights at the boundaries.
     * If there are fewer than two boundaries, then this is equivalent to
     * the default constructor and creates a distribution that
     * produces values uniformly distributed in the range [0, 1).
     *
     * The values of the interval boundaries must be strictly
     * increasing, and the number of weights must be equal to
     * the number of interval boundaries.  If there are extra
     * weights, they are ignored.
     *
     * For example,
     *
     * @code
     * double intervals[] = { 0.0, 1.0, 2.0 };
     * double weights[] = { 0.0, 1.0, 0.0 };
     * piecewise_constant_distribution<> dist(
     *     &intervals[0], &intervals[0] + 3, &weights[0]);
     * @endcode
     *
     * produces a triangle distribution.
     */
    template<class IntervalIter, class WeightIter>
    piecewise_linear_distribution(IntervalIter first_interval,
                                  IntervalIter last_interval,
                                  WeightIter first_weight)
      : _intervals(first_interval, last_interval)
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            _weights.reserve(_intervals.size());
            for(std::size_t i = 0; i < _intervals.size(); ++i) {
                _weights.push_back(*first_weight++);
            }
            init();
        }
    }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    /**
     * Constructs a piecewise_linear_distribution from an
     * initializer_list containing the interval boundaries
     * and a unary function specifying the weights.  Each
     * weight is determined by calling the function at the
     * corresponding interval boundary.
     *
     * If the initializer_list contains fewer than two elements,
     * this is equivalent to the default constructor and the
     * distribution will produce values uniformly distributed
     * in the range [0, 1).
     */
    template<class T, class F>
    piecewise_linear_distribution(std::initializer_list<T> il, F f)
      : _intervals(il.begin(), il.end())
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            _weights.reserve(_intervals.size());
            for(typename std::vector<RealType>::const_iterator
                iter = _intervals.begin(), end = _intervals.end();
                iter != end; ++iter)
            {
                _weights.push_back(f(*iter));
            }
            init();
        }
    }
#endif
    /**
     * Constructs a piecewise_linear_distribution from Boost.Range
     * ranges holding the interval boundaries and the weights.  If
     * there are fewer than two interval boundaries, this is equivalent
     * to the default constructor and the distribution will produce
     * values uniformly distributed in the range [0, 1).  The
     * number of weights must be equal to the number of
     * interval boundaries.
     */
    template<class IntervalsRange, class WeightsRange>
    piecewise_linear_distribution(const IntervalsRange& intervals_arg,
                                  const WeightsRange& weights_arg)
      : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
        _weights(boost::begin(weights_arg), boost::end(weights_arg))
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            init();
        }
    }
    /**
     * Constructs a piecewise_linear_distribution that approximates a
     * function.  The range of the distribution is [xmin, xmax).  This
     * range is divided into nw equally sized intervals and the weights
     * are found by calling the unary function f on the interval boundaries.
     */
    template<class F>
    piecewise_linear_distribution(std::size_t nw,
                                  RealType xmin,
                                  RealType xmax,
                                  F f)
    {
        if(nw == 0) { nw = 1; }
        RealType delta = (xmax - xmin) / nw;
        _intervals.reserve(nw + 1);
        for(std::size_t i = 0; i < nw; ++i) {
            RealType x = xmin + i * delta;
            _intervals.push_back(x);
            _weights.push_back(f(x));
        }
        _intervals.push_back(xmax);
        _weights.push_back(f(xmax));
        init();
    }
    /**
     * Constructs a piecewise_linear_distribution from its parameters.
     */
    explicit piecewise_linear_distribution(const param_type& parm)
      : _intervals(parm._intervals),
        _weights(parm._weights)
    {
        init();
    }

    /**
     * Returns a value distributed according to the parameters of the
     * piecewise_linear_distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng) const
    {
        std::size_t i = _bins(urng);
        bool is_in_rectangle = (i % 2 == 0);
        i = i / 2;
        uniform_real<RealType> dist(_intervals[i], _intervals[i+1]);
        if(is_in_rectangle) {
            return dist(urng);
        } else if(_weights[i] < _weights[i+1]) {
            return (std::max)(dist(urng), dist(urng));
        } else {
            return (std::min)(dist(urng), dist(urng));
        }
    }

    /**
     * Returns a value distributed according to the parameters
     * specified by param.
     */
    template<class URNG>
    RealType operator()(URNG& urng, const param_type& parm) const
    {
        return piecewise_linear_distribution(parm)(urng);
    }

    /** Returns the smallest value that the distribution can produce. */
    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return _intervals.front(); }
    /** Returns the largest value that the distribution can produce. */
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return _intervals.back(); }

    /**
     * Returns a vector containing the probability densities
     * at the interval boundaries.
     */
    std::vector<RealType> densities() const
    {
        RealType sum = static_cast<RealType>(0);
        for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
            RealType width = _intervals[i + 1] - _intervals[i];
            sum += (_weights[i] + _weights[i + 1]) * width / 2;
        }
        std::vector<RealType> result;
        result.reserve(_weights.size());
        for(typename std::vector<RealType>::const_iterator
            iter = _weights.begin(), end = _weights.end();
            iter != end; ++iter)
        {
            result.push_back(*iter / sum);
        }
        return result;
    }
    /**  Returns a vector containing the interval boundaries. */
    std::vector<RealType> intervals() const { return _intervals; }

    /** Returns the parameters of the distribution. */
    param_type param() const
    {
        return param_type(_intervals, _weights);
    }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        std::vector<RealType> new_intervals(parm._intervals);
        std::vector<RealType> new_weights(parm._weights);
        init(new_intervals, new_weights);
        _intervals.swap(new_intervals);
        _weights.swap(new_weights);
    }

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

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

    /** Reads a distribution from a @c std::istream */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
        is, piecewise_linear_distribution, pld)
    {
        param_type parm;
        if(is >> parm) {
            pld.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(
        piecewise_linear_distribution, lhs,  rhs)
    {
        return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights;
    }
    /**
     * Returns true if the two distributions may return different
     * sequences of values, when passed equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution)

private:

    /// @cond \show_private

    void init(const std::vector<RealType>& intervals_arg,
              const std::vector<RealType>& weights_arg)
    {
        using std::abs;
        std::vector<RealType> bin_weights;
        bin_weights.reserve((intervals_arg.size() - 1) * 2);
        for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) {
            RealType width = intervals_arg[i + 1] - intervals_arg[i];
            RealType w1 = weights_arg[i];
            RealType w2 = weights_arg[i + 1];
            bin_weights.push_back((std::min)(w1, w2) * width);
            bin_weights.push_back(abs(w1 - w2) * width / 2);
        }
        typedef discrete_distribution<std::size_t, RealType> bins_type;
        typename bins_type::param_type bins_param(bin_weights);
        _bins.param(bins_param);
    }

    void init()
    {
        init(_intervals, _weights);
    }

    void default_init()
    {
        _intervals.clear();
        _intervals.push_back(RealType(0));
        _intervals.push_back(RealType(1));
        _weights.clear();
        _weights.push_back(RealType(1));
        _weights.push_back(RealType(1));
        init();
    }

    discrete_distribution<std::size_t, RealType> _bins;
    std::vector<RealType> _intervals;
    std::vector<RealType> _weights;

    /// @endcond
};

}
}

#endif


code:
/* boost random/poisson_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_POISSON_DISTRIBUTION_HPP
#define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP

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

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

namespace boost {
namespace random {

namespace detail {

template<class RealType>
struct poisson_table {
    static RealType value[10];
};

template<class RealType>
RealType poisson_table<RealType>::value[10] = {
    0.0,
    0.0,
    0.69314718055994529,
    1.7917594692280550,
    3.1780538303479458,
    4.7874917427820458,
    6.5792512120101012,
    8.5251613610654147,
    10.604602902745251,
    12.801827480081469
};

}

/**
 * An instantiation of the class template @c poisson_distribution is a
 * model of \random_distribution.  The poisson distribution has
 * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$
 *
 * This implementation is based on the PTRD algorithm described
 *
 *  @blockquote
 *  "The transformed rejection method for generating Poisson random variables",
 *  Wolfgang Hormann, Insurance: Mathematics and Economics
 *  Volume 12, Issue 1, February 1993, Pages 39-45
 *  @endblockquote
 */
template<class IntType = int, class RealType = double>
class poisson_distribution {
public:
    typedef IntType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef poisson_distribution distribution_type;
        /**
         * Construct a param_type object with the parameter "mean"
         *
         * Requires: mean > 0
         */
        explicit param_type(RealType mean_arg = RealType(1))
          : _mean(mean_arg)
        {
            BOOST_ASSERT(_mean > 0);
        }
        /* Returns the "mean" parameter of the distribution. */
        RealType mean() const { return _mean; }
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
        /** Writes the parameters of the 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 param_type& parm)
        {
            os << parm._mean;
            return os;
        }

        /** Reads the parameters of the distribution 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._mean;
            return is;
        }
#endif
        /** Returns true if the parameters have the same values. */
        friend bool operator==(const param_type& lhs, const param_type& rhs)
        {
            return lhs._mean == rhs._mean;
        }
        /** Returns true if the parameters have different values. */
        friend bool operator!=(const param_type& lhs, const param_type& rhs)
        {
            return !(lhs == rhs);
        }
    private:
        RealType _mean;
    };

    /**
     * Constructs a @c poisson_distribution with the parameter @c mean.
     *
     * Requires: mean > 0
     */
    explicit poisson_distribution(RealType mean_arg = RealType(1))
      : _mean(mean_arg)
    {
        BOOST_ASSERT(_mean > 0);
        init();
    }

    /**
     * Construct an @c poisson_distribution object from the
     * parameters.
     */
    explicit poisson_distribution(const param_type& parm)
      : _mean(parm.mean())
    {
        init();
    }

    /**
     * Returns a random variate distributed according to the
     * poisson distribution.
     */
    template<class URNG>
    IntType operator()(URNG& urng) const
    {
        if(use_inversion()) {
            return invert(urng);
        } else {
            return generate(urng);
        }
    }

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

    /** Returns the "mean" parameter of the distribution. */
    RealType mean() const { return _mean; }

    /** Returns the smallest value that the distribution can produce. */
    IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 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(_mean); }
    /** Sets parameters of the distribution. */
    void param(const param_type& parm)
    {
        _mean = parm.mean();
        init();
    }

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

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
    /** Writes the parameters of the 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 poisson_distribution& pd)
    {
        os << pd.param();
        return os;
    }

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

    /** Returns true if the two distributions will produce the same
        sequence of values, given equal generators. */
    friend bool operator==(const poisson_distribution& lhs,
                           const poisson_distribution& rhs)
    {
        return lhs._mean == rhs._mean;
    }
    /** Returns true if the two distributions could produce different
        sequences of values, given equal generators. */
    friend bool operator!=(const poisson_distribution& lhs,
                           const poisson_distribution& rhs)
    {
        return !(lhs == rhs);
    }

private:

    /// @cond show_private

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

    bool use_inversion() const
    {
        return _mean < 10;
    }

    static RealType log_factorial(IntType k)
    {
        BOOST_ASSERT(k >= 0);
        BOOST_ASSERT(k < 10);
        return detail::poisson_table<RealType>::value[k];
    }

    void init()
    {
        using std::sqrt;
        using std::exp;

        if(use_inversion()) {
            _exp_mean = exp(-_mean);
        } else {
            _ptrd.smu = sqrt(_mean);
            _ptrd.b = 0.931 + 2.53 * _ptrd.smu;
            _ptrd.a = -0.059 + 0.02483 * _ptrd.b;
            _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4);
            _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2);
        }
    }

    template<class URNG>
    IntType generate(URNG& urng) const
    {
        using std::floor;
        using std::abs;
        using std::log;

        while(true) {
            RealType u;
            RealType v = uniform_01<RealType>()(urng);
            if(v <= 0.86 * _ptrd.v_r) {
                u = v / _ptrd.v_r - 0.43;
                return static_cast<IntType>(floor(
                    (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445));
            }

            if(v >= _ptrd.v_r) {
                u = uniform_01<RealType>()(urng) - 0.5;
            } else {
                u = v/_ptrd.v_r - 0.93;
                u = ((u < 0)? -0.5 : 0.5) - u;
                v = uniform_01<RealType>()(urng) * _ptrd.v_r;
            }

            RealType us = 0.5 - abs(u);
            if(us < 0.013 && v > us) {
                continue;
            }

            RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445);
            v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b);

            RealType log_sqrt_2pi = 0.91893853320467267;

            if(k >= 10) {
                if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k)
                               - _mean
                               - log_sqrt_2pi
                               + k
                               - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) {
                    return static_cast<IntType>(k);
                }
            } else if(k >= 0) {
                if(log(v) <= k*log(_mean)
                           - _mean
                           - log_factorial(static_cast<IntType>(k))) {
                    return static_cast<IntType>(k);
                }
            }
        }
    }

    template<class URNG>
    IntType invert(URNG& urng) const
    {
        RealType p = _exp_mean;
        IntType x = 0;
        RealType u = uniform_01<RealType>()(urng);
        while(u > p) {
            u = u - p;
            ++x;
            p = _mean * p / x;
        }
        return x;
    }

    RealType _mean;

    union {
        // for ptrd
        struct {
            RealType v_r;
            RealType a;
            RealType b;
            RealType smu;
            RealType inv_alpha;
        } _ptrd;
        // for inversion
        RealType _exp_mean;
    };

    /// @endcond
};

} // namespace random

using random::poisson_distribution;

} // namespace boost

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

#endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP

license: mit


code:
/* boost random/random_number_generator.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * 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_RANDOM_NUMBER_GENERATOR_HPP
#define BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP

#include <boost/assert.hpp>
#include <boost/random/uniform_int_distribution.hpp>

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

namespace boost {
namespace random {

/**
 * Instantiations of class template random_number_generator model a
 * RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle]). On
 * each invocation, it returns a uniformly distributed integer in
 * the range [0..n).
 *
 * The template parameter IntType shall denote some integer-like value type.
 */
template<class URNG, class IntType = long>
class random_number_generator
{
public:
    typedef URNG base_type;
    typedef IntType argument_type;
    typedef IntType result_type;
    /**
     * Constructs a random_number_generator functor with the given
     * \uniform_random_number_generator as the underlying source of
     * random numbers.
     */
    random_number_generator(base_type& rng) : _rng(rng) {}

    // compiler-generated copy ctor is fine
    // assignment is disallowed because there is a reference member

    /**
     * Returns a value in the range [0, n)
     */
    result_type operator()(argument_type n)
    {
        BOOST_ASSERT(n > 0);
        return uniform_int_distribution<IntType>(0, n-1)(_rng);
    }

private:
    base_type& _rng;
};

} // namespace random

using random::random_number_generator;

} // namespace boost

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

#endif // BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP


code:
/* boost random/ranlux.hpp header file
 *
 * Copyright Jens Maurer 2002
 * 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  created
 */

#ifndef BOOST_RANDOM_RANLUX_HPP
#define BOOST_RANDOM_RANLUX_HPP

#include <boost/config.hpp>
#include <boost/random/subtract_with_carry.hpp>
#include <boost/random/discard_block.hpp>

namespace boost {

namespace random {

namespace detail {
/**
 * The ranlux family of generators are described in
 *
 *  @blockquote
 *  "A portable high-quality random number generator for lattice field theory
 *  calculations", M. Luescher, Computer Physics Communications, 79 (1994)
 *  pp 100-110.
 *  @endblockquote
 *
 * The levels are given in
 *
 *  @blockquote
 *  "RANLUX: A Fortran implementation of the high-quality
 *  pseudorandom number generator of Luescher", F. James,
 *  Computer Physics Communications 79 (1994) 111-114
 *  @endblockquote
 */
class ranlux_documentation {};
}

typedef subtract_with_carry_engine<uint32_t, 24, 10, 24> ranlux_base;
typedef subtract_with_carry_01_engine<float, 24, 10, 24> ranlux_base_01;
typedef subtract_with_carry_01_engine<double, 48, 10, 24> ranlux64_base_01;

/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux_base, 223, 24> ranlux3;
/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux_base, 389, 24> ranlux4;

/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux_base_01, 223, 24> ranlux3_01;
/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux_base_01, 389, 24> ranlux4_01;

/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux64_base_01, 223, 24> ranlux64_3_01;
/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux64_base_01, 389, 24> ranlux64_4_01;

#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
typedef subtract_with_carry_engine<uint64_t, 48, 10, 24> ranlux64_base;
/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux64_base, 223, 24> ranlux64_3;
/** @copydoc boost::random::detail::ranlux_documentation */
typedef discard_block_engine<ranlux64_base, 389, 24> ranlux64_4;
#endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */

typedef subtract_with_carry_engine<uint32_t, 24, 10, 24> ranlux24_base;
typedef subtract_with_carry_engine<uint64_t, 48, 5, 12> ranlux48_base;

typedef discard_block_engine<ranlux24_base, 223, 23> ranlux24;
#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
typedef discard_block_engine<ranlux48_base, 389, 11> ranlux48;
#endif
}

using random::ranlux3;
using random::ranlux4;
using random::ranlux3_01;
using random::ranlux4_01;
using random::ranlux64_3_01;
using random::ranlux64_4_01;
#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
using random::ranlux64_3;
using random::ranlux64_4;
#endif

} // namespace boost

#endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP


code:
/* boost random/seed_seq.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_SEED_SEQ_HPP
#define BOOST_RANDOM_SEED_SEQ_HPP

#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/range/begin.hpp>
#include <boost/range/end.hpp>
#include <cstddef>
#include <vector>
#include <algorithm>
#include <iterator>

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
#include <initializer_list>
#endif

namespace boost {
namespace random {

/**
 * The class @c seed_seq stores a sequence of 32-bit words
 * for seeding a \pseudo_random_number_generator.  These
 * words will be combined to fill the entire state of the
 * generator.
 */
class seed_seq {
public:
    typedef boost::uint_least32_t result_type;

    /** Initializes a seed_seq to hold an empty sequence. */
    seed_seq() {}
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    /** Initializes the sequence from an initializer_list. */
    template<class T>
    seed_seq(const std::initializer_list<T>& il) : v(il.begin(), il.end()) {}
#endif
    /** Initializes the sequence from an iterator range. */
    template<class Iter>
    seed_seq(Iter first, Iter last) : v(first, last) {}
    /** Initializes the sequence from Boost.Range range. */
    template<class Range>
    explicit seed_seq(const Range& range)
      : v(boost::begin(range), boost::end(range)) {}

    /**
     * Fills a range with 32-bit values based on the stored sequence.
     *
     * Requires: Iter must be a Random Access Iterator whose value type
     * is an unsigned integral type at least 32 bits wide.
     */
    template<class Iter>
    void generate(Iter first, Iter last) const
    {
        typedef typename std::iterator_traits<Iter>::value_type value_type;
        std::fill(first, last, static_cast<value_type>(0x8b8b8b8bu));
        std::size_t s = v.size();
        std::size_t n = last - first;
        std::size_t t =
            (n >= 623) ? 11 :
            (n >=  68) ?  7 :
            (n >=  39) ?  5 :
            (n >=   7) ?  3 :
            (n - 1)/2;
        std::size_t p = (n - t) / 2;
        std::size_t q = p + t;
        std::size_t m = (std::max)(s+1, n);
        value_type mask = 0xffffffffu;
        for(std::size_t k = 0; k < m; ++k) {
            value_type r1 = static_cast<value_type>
                (*(first + k%n) ^ *(first + (k+p)%n) ^ *(first + (k+n-1)%n));
            r1 = r1 ^ (r1 >> 27);
            r1 = (r1 * 1664525u) & mask;
            value_type r2 = static_cast<value_type>(r1 +
                ((k == 0) ? s :
                 (k <= s) ? k % n + v[k - 1] :
                 (k % n)));
            *(first + (k+p)%n) = (*(first + (k+p)%n) + r1) & mask;
            *(first + (k+q)%n) = (*(first + (k+q)%n) + r2) & mask;
            *(first + k%n) = r2;
        }
        for(std::size_t k = m; k < m + n; ++k) {
            value_type r3 = static_cast<value_type>
                ((*(first + k%n) + *(first + (k+p)%n) + *(first + (k+n-1)%n))
                & mask);
            r3 = r3 ^ (r3 >> 27);
            r3 = (r3 * 1566083941u) & mask;
            value_type r4 = static_cast<value_type>(r3 - k%m);
            *(first + (k+p)%n) ^= r3;
            *(first + (k+q)%n) ^= r4;
            *(first + k%n) = r4;
        }
    }
    /** Returns the size of the sequence. */
    std::size_t size() const { return v.size(); }
    /** Writes the stored sequence to iter. */
    template<class Iter>
    void param(Iter out) { std::copy(v.begin(), v.end(), out); }
private:
    std::vector<result_type> v;
};

}
}

#endif


code:
/* boost random/shuffle_order.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * 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_SHUFFLE_ORDER_HPP
#define BOOST_RANDOM_SHUFFLE_ORDER_HPP

#include <iostream>
#include <algorithm>     // std::copy
#include <cassert>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/cstdint.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/random/linear_congruential.hpp>

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

namespace boost {
namespace random {

/**
 * Instatiations of class template @c shuffle_order_engine model a
 * \pseudo_random_number_generator. It mixes the output
 * of some (usually \linear_congruential_engine)
 * \uniform_random_number_generator to get better statistical properties.
 * The algorithm is described in
 *
 *  @blockquote
 *  "Improving a poor random number generator", Carter Bays
 *  and S.D. Durham, ACM Transactions on Mathematical Software,
 *  Vol 2, No. 1, March 1976, pp. 59-64.
 *  http://doi.acm.org/10.1145/355666.355670
 *  @endblockquote
 *
 * The output of the base generator is buffered in an array of
 * length k. Every output X(n) has a second role: It gives an
 * index into the array where X(n+1) will be retrieved. Used
 * array elements are replaced with fresh output from the base
 * generator.
 *
 * Template parameters are the base generator and the array
 * length k, which should be around 100.
 */
template<class UniformRandomNumberGenerator, std::size_t k>
class shuffle_order_engine
{
public:
    typedef UniformRandomNumberGenerator base_type;
    typedef typename base_type::result_type result_type;

    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    BOOST_STATIC_CONSTANT(std::size_t, buffer_size = k);
    BOOST_STATIC_CONSTANT(std::size_t, table_size = k);

    BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
    BOOST_STATIC_ASSERT(k > 0);

    /**
     * Constructs a @c shuffle_order_engine by invoking the
     * default constructor of the base generator.
     *
     * Complexity: Exactly k+1 invocations of the base generator.
     */
    shuffle_order_engine() : _rng() { init(); }
    /**
     * Constructs a @c shuffle_output_engine by invoking the one-argument
     * constructor of the base generator with the parameter seed.
     *
     * Complexity: Exactly k+1 invocations of the base generator.
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(shuffle_order_engine,
                                               result_type, s)
    { _rng.seed(s); init(); }
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(shuffle_order_engine, SeedSeq, seq)
    { _rng.seed(seq); init(); }
    /**
     * Constructs a @c shuffle_output_engine by using a copy
     * of the provided generator.
     *
     * Precondition: The template argument UniformRandomNumberGenerator
     * shall denote a CopyConstructible type.
     *
     * Complexity: Exactly k+1 invocations of the base generator.
     */
    explicit shuffle_order_engine(const base_type & rng) : _rng(rng) { init(); }

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
    explicit shuffle_order_engine(base_type&& rng) : _rng(rng) { init(); }
#endif

    template<class It> shuffle_order_engine(It& first, It last)
      : _rng(first, last) { init(); }
    void seed() { _rng.seed(); init(); }
    /**
     * Invokes the one-argument seed method of the base generator
     * with the parameter seed and re-initializes the internal buffer array.
     *
     * Complexity: Exactly k+1 invocations of the base generator.
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(shuffle_order_engine,
        result_type, seed_arg)
    { _rng.seed(seed_arg); init(); }
    /**
     * Invokes the one-argument seed method of the base generator
     * with the parameter seq and re-initializes the internal buffer array.
     *
     * Complexity: Exactly k+1 invocations of the base generator.
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(shuffle_order_engine, SeedSeq, seq)
    { _rng.seed(seq); init(); }
    template<class It> void seed(It& first, It last)
    { _rng.seed(first, last); init(); }

    const base_type& base() const { return _rng; }

    result_type operator()() {
        // calculating the range every time may seem wasteful.  However, this
        // makes the information locally available for the optimizer.
        typedef typename boost::random::traits::make_unsigned<result_type>::type base_unsigned;
        const base_unsigned brange =
            detail::subtract<result_type>()((max)(), (min)());
        const base_unsigned off =
            detail::subtract<result_type>()(y, (min)());

        base_unsigned j;
        if(k == 1) {
            j = 0;
        } else if(brange < (std::numeric_limits<base_unsigned>::max)() / k) {
            // try to do it in the native type if we know that it won't
            // overflow
            j = k * off / (brange + 1);
        } else if(brange < (std::numeric_limits<uintmax_t>::max)() / k) {
            // Otherwise try to use uint64_t
            j = static_cast<base_unsigned>(
                static_cast<uintmax_t>(off) * k /
                (static_cast<uintmax_t>(brange) + 1));
        } else {
            boost::uintmax_t divisor =
                static_cast<boost::uintmax_t>(brange) + 1;
            j = static_cast<base_unsigned>(detail::muldiv(off, k, divisor));
        }
        // assert(0 <= j && j < k);
        y = v[j];
        v[j] = _rng();
        return y;
    }

    /** Advances the generator by z steps. */
    void discard(boost::uintmax_t z)
    {
        for(boost::uintmax_t j = 0; j < z; ++j) {
            (*this)();
        }
    }

    /** Fills a range with pseudo-random values. */
    template<class Iter>
    void generate(Iter first, Iter last)
    { detail::generate_from_int(*this, first, last); }

    /** Returns the smallest value that the generator can produce. */
    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (base_type::min)(); }
    /** Returns the largest value that the generator can produce. */
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return (base_type::max)(); }

    /** Writes a @c shuffle_order_engine to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, shuffle_order_engine, s)
    {
        os << s._rng;
        for(std::size_t i = 0; i < k; ++i)
            os << ' ' << s.v[i];
        os << ' ' << s.y;
        return os;
    }

    /** Reads a @c shuffle_order_engine from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, shuffle_order_engine, s)
    {
        is >> s._rng;
        for(std::size_t i = 0; i < k; ++i)
            is >> std::ws >> s.v[i];
        is >> std::ws >> s.y;
        return is;
    }

    /** Returns true if the two generators will produce identical sequences. */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(shuffle_order_engine, x, y)
    { return x._rng == y._rng && x.y == y.y && std::equal(x.v, x.v+k, y.v); }
    /** Returns true if the two generators will produce different sequences. */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(shuffle_order_engine)

private:

    /// \cond show_private

    void init()
    {
        // we cannot use std::generate, because it uses pass-by-value for _rng
        for(result_type * p = v; p != v+k; ++p)
            *p = _rng();
        y = _rng();
    }

    /// \endcond

    base_type _rng;
    result_type v[k];
    result_type y;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class URNG, std::size_t k>
const bool shuffle_order_engine<URNG, k>::has_fixed_range;
template<class URNG, std::size_t k>
const std::size_t shuffle_order_engine<URNG, k>::table_size;
template<class URNG, std::size_t k>
const std::size_t shuffle_order_engine<URNG, k>::buffer_size;
#endif

/**
 * According to Harry Erwin (private e-mail), the specialization
 * @c kreutzer1986 was suggested in:
 *
 * @blockquote
 * "System Simulation: Programming Styles and Languages (International
 * Computer Science Series)", Wolfgang Kreutzer, Addison-Wesley, December 1986.
 * @endblockquote
 */
typedef shuffle_order_engine<
    linear_congruential_engine<uint32_t, 1366, 150889, 714025>,
    97> kreutzer1986;

/**
 * The specialization @c knuth_b is specified by the C++ standard.
 * It is described in
 *
 * @blockquote
 * "The Art of Computer Programming, Second Edition, Volume 2,
 * Seminumerical Algorithms", Donald Knuth, Addison-Wesley, 1981.
 * @endblockquote
 */
typedef shuffle_order_engine<minstd_rand0, 256> knuth_b;

} // namespace random

using random::kreutzer1986;

} // namespace boost

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

#endif // BOOST_RANDOM_SHUFFLE_OUTPUT_HPP


code:
/* boost random/shuffle_output.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * 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_SHUFFLE_OUTPUT_HPP
#define BOOST_RANDOM_SHUFFLE_OUTPUT_HPP

#include <boost/random/shuffle_order.hpp>

namespace boost {
namespace random {

/// \cond

template<typename URNG, int k,
         typename URNG::result_type val = 0>
class shuffle_output : public shuffle_order_engine<URNG, k>
{
    typedef shuffle_order_engine<URNG, k> base_t;
public:
    typedef typename base_t::result_type result_type;
    shuffle_output() {}
    template<class T>
    explicit shuffle_output(T& arg) : base_t(arg) {}
    template<class T>
    explicit shuffle_output(const T& arg) : base_t(arg) {}
    template<class It>
    shuffle_output(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

}
}

#endif // BOOST_RANDOM_SHUFFLE_OUTPUT_HPP


code:
/* boost random/student_t_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_STUDENT_T_DISTRIBUTION_HPP
#define BOOST_RANDOM_STUDENT_T_DISTRIBUTION_HPP

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

namespace boost {
namespace random {

/**
 * The Student t distribution is a real valued distribution with one
 * parameter n, the number of degrees of freedom.
 *
 * It has \f$\displaystyle p(x) =
 *   \frac{1}{\sqrt{n\pi}}
 *   \frac{\Gamma((n+1)/2)}{\Gamma(n/2)}
 *   \left(1+\frac{x^2}{n}\right)^{-(n+1)/2}
 * \f$.
 */
template<class RealType = double>
class student_t_distribution {
public:
    typedef RealType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef student_t_distribution distribution_type;

        /**
         * Constructs a @c param_type with "n" degrees of freedom.
         *
         * Requires: n > 0
         */
        explicit param_type(RealType n_arg = RealType(1.0))
          : _n(n_arg)
        {}

        /** Returns the number of degrees of freedom 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._n; return os; }

        /** Reads a @c param_type from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        { is >> 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._n == rhs._n; }

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

    private:
        RealType _n;
    };

    /**
     * Constructs an @c student_t_distribution with "n" degrees of freedom.
     *
     * Requires: n > 0
     */
    explicit student_t_distribution(RealType n_arg = RealType(1.0))
      : _normal(), _chi_squared(n_arg)
    {}
    /** Constructs an @c student_t_distribution from its parameters. */
    explicit student_t_distribution(const param_type& parm)
      : _normal(), _chi_squared(parm.n())
    {}

    /**
     * Returns a random variate distributed according to the
     * Student t distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng)
    {
        using std::sqrt;
        return _normal(urng) / sqrt(_chi_squared(urng) / n());
    }

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

    /** Returns the number of degrees of freedom. */
    RealType n() const { return _chi_squared.n(); }

    /** 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(n()); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        typedef chi_squared_distribution<RealType> chi_squared_type;
        typename chi_squared_type::param_type chi_squared_param(parm.n());
        _chi_squared.param(chi_squared_param);
    }

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

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

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

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

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

private:
    normal_distribution<RealType> _normal;
    chi_squared_distribution<RealType> _chi_squared;
};

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_STUDENT_T_DISTRIBUTION_HPP


code:
/* boost random/subtract_with_carry.hpp header file
 *
 * Copyright Jens Maurer 2002
 * 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
 *  2002-03-02  created
 */

#ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
#define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP

#include <boost/config/no_tr1/cmath.hpp>         // std::pow
#include <iostream>
#include <algorithm>     // std::equal
#include <stdexcept>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/cstdint.hpp>
#include <boost/static_assert.hpp>
#include <boost/integer/static_log2.hpp>
#include <boost/integer/integer_mask.hpp>
#include <boost/detail/workaround.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/seed_impl.hpp>
#include <boost/random/detail/generator_seed_seq.hpp>
#include <boost/random/linear_congruential.hpp>

namespace boost {
namespace random {

namespace detail {

struct subtract_with_carry_discard
{
    template<class Engine>
    static void apply(Engine& eng, boost::uintmax_t z)
    {
        typedef typename Engine::result_type IntType;
        const std::size_t short_lag = Engine::short_lag;
        const std::size_t long_lag = Engine::long_lag;
        std::size_t k = eng.k;
        IntType carry = eng.carry;
        if(k != 0) {
            // increment k until it becomes 0.
            if(k < short_lag) {
                std::size_t limit = (short_lag - k) < z?
                    short_lag : (k + static_cast<std::size_t>(z));
                for(std::size_t j = k; j < limit; ++j) {
                    carry = eng.do_update(j, j + long_lag - short_lag, carry);
                }
            }
            std::size_t limit = (long_lag - k) < z?
                long_lag : (k + static_cast<std::size_t>(z));
            std::size_t start = (k < short_lag ? short_lag : k);
            for(std::size_t j = start; j < limit; ++j) {
                carry = eng.do_update(j, j - short_lag, carry);
            }
        }

        k = ((z % long_lag) + k) % long_lag;

        if(k < z) {
            // main loop: update full blocks from k = 0 to long_lag
            for(std::size_t i = 0; i < (z - k) / long_lag; ++i) {
                for(std::size_t j = 0; j < short_lag; ++j) {
                    carry = eng.do_update(j, j + long_lag - short_lag, carry);
                }
                for(std::size_t j = short_lag; j < long_lag; ++j) {
                    carry = eng.do_update(j, j - short_lag, carry);
                }
            }

            // Update the last partial block
            std::size_t limit = short_lag < k? short_lag : k;
            for(std::size_t j = 0; j < limit; ++j) {
                carry = eng.do_update(j, j + long_lag - short_lag, carry);
            }
            for(std::size_t j = short_lag; j < k; ++j) {
                carry = eng.do_update(j, j - short_lag, carry);
            }
        }
        eng.carry = carry;
        eng.k = k;
    }
};

}

/**
 * Instantiations of @c subtract_with_carry_engine model a
 * \pseudo_random_number_generator.  The algorithm is
 * described in
 *
 *  @blockquote
 *  "A New Class of Random Number Generators", George
 *  Marsaglia and Arif Zaman, Annals of Applied Probability,
 *  Volume 1, Number 3 (1991), 462-480.
 *  @endblockquote
 */
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
class subtract_with_carry_engine
{
public:
    typedef IntType result_type;
    BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
    BOOST_STATIC_CONSTANT(std::size_t, long_lag = r);
    BOOST_STATIC_CONSTANT(std::size_t, short_lag = s);
    BOOST_STATIC_CONSTANT(uint32_t, default_seed = 19780503u);

    // Required by the old Boost.Random concepts
    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    // Backwards compatibility
    BOOST_STATIC_CONSTANT(result_type, modulus = (result_type(1) << w));

    BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);

    /**
     * Constructs a new @c subtract_with_carry_engine and seeds
     * it with the default seed.
     */
    subtract_with_carry_engine() { seed(); }
    /**
     * Constructs a new @c subtract_with_carry_engine and seeds
     * it with @c value.
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_engine,
                                               IntType, value)
    { seed(value); }
    /**
     * Constructs a new @c subtract_with_carry_engine and seeds
     * it with values produced by @c seq.generate().
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(subtract_with_carry_engine,
                                             SeedSeq, seq)
    { seed(seq); }
    /**
     * Constructs a new @c subtract_with_carry_engine and seeds
     * it with values from a range.  first is updated to point
     * one past the last value consumed.  If there are not
     * enough elements in the range to fill the entire state of
     * the generator, throws @c std::invalid_argument.
     */
    template<class It> subtract_with_carry_engine(It& first, It last)
    { seed(first,last); }

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

    /** Seeds the generator with the default seed. */
    void seed() { seed(default_seed); }
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_engine,
                                        IntType, value)
    {
        typedef linear_congruential_engine<uint32_t,40014,0,2147483563> gen_t;
        gen_t intgen(static_cast<boost::uint32_t>(value == 0 ? default_seed : value));
        detail::generator_seed_seq<gen_t> gen(intgen);
        seed(gen);
    }

    /** Seeds the generator with values produced by @c seq.generate(). */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(subtract_with_carry, SeedSeq, seq)
    {
        detail::seed_array_int<w>(seq, x);
        carry = (x[long_lag-1] == 0);
        k = 0;
    }

    /**
     * Seeds the generator with values from a range.  Updates @c first to
     * point one past the last consumed value.  If the range does not
     * contain enough elements to fill the entire state of the generator,
     * throws @c std::invalid_argument.
     */
    template<class It>
    void seed(It& first, It last)
    {
        detail::fill_array_int<w>(first, last, x);
        carry = (x[long_lag-1] == 0);
        k = 0;
    }

    /** Returns the smallest value that the generator can produce. */
    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return 0; }
    /** Returns the largest value that the generator can produce. */
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return boost::low_bits_mask_t<w>::sig_bits; }

    /** Returns the next value of the generator. */
    result_type operator()()
    {
        std::size_t short_index =
            (k < short_lag)?
                (k + long_lag - short_lag) :
                (k - short_lag);
        carry = do_update(k, short_index, carry);
        IntType result = x[k];
        ++k;
        if(k >= long_lag)
            k = 0;
        return result;
    }

    /** Advances the state of the generator by @c z. */
    void discard(boost::uintmax_t z)
    {
        detail::subtract_with_carry_discard::apply(*this, z);
    }

    /** Fills a range with random values. */
    template<class It>
    void generate(It first, It last)
    { detail::generate_from_int(*this, first, last); }

    /** Writes a @c subtract_with_carry_engine to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, subtract_with_carry_engine, f)
    {
        for(unsigned int j = 0; j < f.long_lag; ++j)
            os << f.compute(j) << ' ';
        os << f.carry;
        return os;
    }

    /** Reads a @c subtract_with_carry_engine from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, subtract_with_carry_engine, f)
    {
        for(unsigned int j = 0; j < f.long_lag; ++j)
            is >> f.x[j] >> std::ws;
        is >> f.carry;
        f.k = 0;
        return is;
    }

    /**
     * Returns true if the two generators will produce identical
     * sequences of values.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(subtract_with_carry_engine, x, y)
    {
        for(unsigned int j = 0; j < r; ++j)
            if(x.compute(j) != y.compute(j))
                return false;
        return true;
    }

    /**
     * Returns true if the two generators will produce different
     * sequences of values.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(subtract_with_carry_engine)

private:
    /// \cond show_private
    // returns x(i-r+index), where index is in 0..r-1
    IntType compute(unsigned int index) const
    {
        return x[(k+index) % long_lag];
    }

    friend struct detail::subtract_with_carry_discard;

    IntType do_update(std::size_t current, std::size_t short_index, IntType carry)
    {
        IntType delta;
        IntType temp = x[current] + carry;
        if (x[short_index] >= temp) {
            // x(n) >= 0
            delta =  x[short_index] - temp;
            carry = 0;
        } else {
            // x(n) < 0
            delta = modulus - temp + x[short_index];
            carry = 1;
        }
        x[current] = delta;
        return carry;
    }
    /// \endcond

    // state representation; next output (state) is x(i)
    //   x[0]  ... x[k] x[k+1] ... x[long_lag-1]     represents
    //  x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1)
    // speed: base: 20-25 nsec
    // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec
    // This state representation makes operator== and save/restore more
    // difficult, because we've already computed "too much" and thus
    // have to undo some steps to get at x(i-r) etc.

    // state representation: next output (state) is x(i)
    //   x[0]  ... x[k] x[k+1]          ... x[long_lag-1]     represents
    //  x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1)
    // speed: base 28 nsec
    // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec
    IntType x[long_lag];
    std::size_t k;
    IntType carry;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const bool subtract_with_carry_engine<IntType, w, s, r>::has_fixed_range;
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const IntType subtract_with_carry_engine<IntType, w, s, r>::modulus;
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_engine<IntType, w, s, r>::word_size;
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_engine<IntType, w, s, r>::long_lag;
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_engine<IntType, w, s, r>::short_lag;
template<class IntType, std::size_t w, std::size_t s, std::size_t r>
const uint32_t subtract_with_carry_engine<IntType, w, s, r>::default_seed;
#endif

// use a floating-point representation to produce values in [0..1)
/**
 * Instantiations of \subtract_with_carry_01_engine model a
 * \pseudo_random_number_generator.  The algorithm is
 * described in
 *
 *  @blockquote
 *  "A New Class of Random Number Generators", George
 *  Marsaglia and Arif Zaman, Annals of Applied Probability,
 *  Volume 1, Number 3 (1991), 462-480.
 *  @endblockquote
 */
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
class subtract_with_carry_01_engine
{
public:
    typedef RealType result_type;
    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
    BOOST_STATIC_CONSTANT(std::size_t, long_lag = r);
    BOOST_STATIC_CONSTANT(std::size_t, short_lag = s);
    BOOST_STATIC_CONSTANT(boost::uint32_t, default_seed = 19780503u);

    BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);

    /** Creates a new \subtract_with_carry_01_engine using the default seed. */
    subtract_with_carry_01_engine() { init_modulus(); seed(); }
    /** Creates a new subtract_with_carry_01_engine and seeds it with value. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_01_engine,
                                               boost::uint32_t, value)
    { init_modulus(); seed(value); }
    /**
     * Creates a new \subtract_with_carry_01_engine and seeds with values
     * produced by seq.generate().
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(subtract_with_carry_01_engine,
                                             SeedSeq, seq)
    { init_modulus(); seed(seq); }
    /**
     * Creates a new \subtract_with_carry_01_engine and seeds it with values
     * from a range.  Advances first to point one past the last consumed
     * value.  If the range does not contain enough elements to fill the
     * entire state, throws @c std::invalid_argument.
     */
    template<class It> subtract_with_carry_01_engine(It& first, It last)
    { init_modulus(); seed(first,last); }

private:
    /// \cond show_private
    void init_modulus()
    {
#ifndef BOOST_NO_STDC_NAMESPACE
        // allow for Koenig lookup
        using std::pow;
#endif
        _modulus = pow(RealType(2), RealType(word_size));
    }
    /// \endcond

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

    /** Seeds the generator with the default seed. */
    void seed() { seed(default_seed); }

    /** Seeds the generator with @c value. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_01_engine,
                                        boost::uint32_t, value)
    {
        typedef linear_congruential_engine<uint32_t, 40014, 0, 2147483563> gen_t;
        gen_t intgen(value == 0 ? default_seed : value);
        detail::generator_seed_seq<gen_t> gen(intgen);
        seed(gen);
    }

    /** Seeds the generator with values produced by @c seq.generate(). */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(subtract_with_carry_01_engine,
                                      SeedSeq, seq)
    {
        detail::seed_array_real<w>(seq, x);
        carry = (x[long_lag-1] ? result_type(0) : result_type(1 / _modulus));
        k = 0;
    }

    /**
     * Seeds the generator with values from a range.  Updates first to
     * point one past the last consumed element.  If there are not
     * enough elements in the range to fill the entire state, throws
     * @c std::invalid_argument.
     */
    template<class It>
    void seed(It& first, It last)
    {
        detail::fill_array_real<w>(first, last, x);
        carry = (x[long_lag-1] ? result_type(0) : result_type(1 / _modulus));
        k = 0;
    }

    /** Returns the smallest value that the generator can produce. */
    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return result_type(0); }
    /** Returns the largest value that the generator can produce. */
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
    { return result_type(1); }

    /** Returns the next value of the generator. */
    result_type operator()()
    {
        std::size_t short_index =
            (k < short_lag) ?
                (k + long_lag - short_lag) :
                (k - short_lag);
        carry = do_update(k, short_index, carry);
        RealType result = x[k];
        ++k;
        if(k >= long_lag)
            k = 0;
        return result;
    }

    /** Advances the state of the generator by @c z. */
    void discard(boost::uintmax_t z)
    { detail::subtract_with_carry_discard::apply(*this, z); }

    /** Fills a range with random values. */
    template<class Iter>
    void generate(Iter first, Iter last)
    { detail::generate_from_real(*this, first, last); }

    /** Writes a \subtract_with_carry_01_engine to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, subtract_with_carry_01_engine, f)
    {
        std::ios_base::fmtflags oldflags =
            os.flags(os.dec | os.fixed | os.left);
        for(unsigned int j = 0; j < f.long_lag; ++j)
            os << (f.compute(j) * f._modulus) << ' ';
        os << (f.carry * f._modulus);
        os.flags(oldflags);
        return os;
    }

    /** Reads a \subtract_with_carry_01_engine from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, subtract_with_carry_01_engine, f)
    {
        RealType value;
        for(unsigned int j = 0; j < long_lag; ++j) {
            is >> value >> std::ws;
            f.x[j] = value / f._modulus;
        }
        is >> value;
        f.carry = value / f._modulus;
        f.k = 0;
        return is;
    }

    /** Returns true if the two generators will produce identical sequences. */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(subtract_with_carry_01_engine, x, y)
    {
        for(unsigned int j = 0; j < r; ++j)
            if(x.compute(j) != y.compute(j))
                return false;
        return true;
    }

    /** Returns true if the two generators will produce different sequences. */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(subtract_with_carry_01_engine)

private:
    /// \cond show_private
    RealType compute(unsigned int index) const
    {
        return x[(k+index) % long_lag];
    }

    friend struct detail::subtract_with_carry_discard;

    RealType do_update(std::size_t current, std::size_t short_index, RealType carry)
    {
        RealType delta = x[short_index] - x[current] - carry;
        if(delta < 0) {
            delta += RealType(1);
            carry = RealType(1)/_modulus;
        } else {
            carry = 0;
        }
        x[current] = delta;
        return carry;
    }
    /// \endcond
    std::size_t k;
    RealType carry;
    RealType x[long_lag];
    RealType _modulus;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
const bool subtract_with_carry_01_engine<RealType, w, s, r>::has_fixed_range;
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::word_size;
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::long_lag;
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::short_lag;
template<class RealType, std::size_t w, std::size_t s, std::size_t r>
const uint32_t subtract_with_carry_01_engine<RealType, w, s, r>::default_seed;
#endif

/// \cond show_deprecated

template<class IntType, IntType m, unsigned s, unsigned r, IntType v>
class subtract_with_carry :
    public subtract_with_carry_engine<IntType,
        boost::static_log2<m>::value, s, r>
{
    typedef subtract_with_carry_engine<IntType,
        boost::static_log2<m>::value, s, r> base_type;
public:
    subtract_with_carry() {}
    BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Gen, gen)
    { seed(gen); }
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry,
                                               IntType, val)
    { seed(val); }
    template<class It>
    subtract_with_carry(It& first, It last) : base_type(first, last) {}
    void seed() { base_type::seed(); }
    BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry, Gen, gen)
    {
        detail::generator_seed_seq<Gen> seq(gen);
        base_type::seed(seq);
    }
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, IntType, val)
    { base_type::seed(val); }
    template<class It>
    void seed(It& first, It last) { base_type::seed(first, last); }
};

template<class RealType, int w, unsigned s, unsigned r, int v = 0>
class subtract_with_carry_01 :
    public subtract_with_carry_01_engine<RealType, w, s, r>
{
    typedef subtract_with_carry_01_engine<RealType, w, s, r> base_type;
public:
    subtract_with_carry_01() {}
    BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry_01, Gen, gen)
    { seed(gen); }
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_01,
                                               uint32_t, val)
    { seed(val); }
    template<class It>
    subtract_with_carry_01(It& first, It last) : base_type(first, last) {}
    void seed() { base_type::seed(); }
    BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry_01, Gen, gen)
    {
        detail::generator_seed_seq<Gen> seq(gen);
        base_type::seed(seq);
    }
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_01, uint32_t, val)
    { base_type::seed(val); }
    template<class It>
    void seed(It& first, It last) { base_type::seed(first, last); }
};

/// \endcond

namespace detail {

template<class Engine>
struct generator_bits;

template<class RealType, std::size_t w, std::size_t s, std::size_t r>
struct generator_bits<subtract_with_carry_01_engine<RealType, w, s, r> > {
    static std::size_t value() { return w; }
};

template<class RealType, int w, unsigned s, unsigned r, int v>
struct generator_bits<subtract_with_carry_01<RealType, w, s, r, v> > {
    static std::size_t value() { return w; }
};

}

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP


code:
/* boost random/taus88.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/libs/random for documentation.
 *
 * $Id$
 *
 */

#ifndef BOOST_RANDOM_TAUS88_HPP
#define BOOST_RANDOM_TAUS88_HPP

#include <boost/random/linear_feedback_shift.hpp>
#include <boost/random/xor_combine.hpp>

namespace boost {
namespace random {

/**
 * The specialization taus88 was suggested in
 *
 *  @blockquote
 *  "Maximally Equidistributed Combined Tausworthe Generators",
 *  Pierre L'Ecuyer, Mathematics of Computation, Volume 65,
 *  Number 213, January 1996, Pages 203-213
 *  @endblockquote
 */
typedef xor_combine_engine<
    xor_combine_engine<
        linear_feedback_shift_engine<uint32_t, 32, 31, 13, 12>, 0,
        linear_feedback_shift_engine<uint32_t, 32, 29, 2, 4>, 0>, 0,
    linear_feedback_shift_engine<uint32_t, 32, 28, 3, 17>, 0> taus88;

} // namespace random

using random::taus88;

} // namespace boost

#endif // BOOST_RANDOM_TAUS88_HPP


code:
/* boost random/traits.hpp header file
 *
 * Copyright John Maddock 2015
 * 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.
 *
 * These traits classes serve two purposes: they are designed to mostly
 * work out of the box for multiprecision types (ie number types that are
 * C++ class types and not integers or floats from type-traits point of view),
 * they are also a potential point of specialization for user-defined
 * number types.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_TRAITS_HPP
#define BOOST_RANDOM_TRAITS_HPP

#include <boost/type_traits/is_signed.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/type_traits/make_unsigned.hpp>
#include <boost/mpl/bool.hpp>
#include <limits>

namespace boost {
namespace random {
namespace traits {
   // \cond show_private
   template <class T, bool intrinsic>
   struct make_unsigned_imp
   {
      typedef typename boost::make_unsigned<T>::type type;
   };
   template <class T>
   struct make_unsigned_imp<T, false>
   {
      BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
      BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_signed == false);
      BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_integer == true);
      typedef T type;
   };
   // \endcond
   /** \brief Converts the argument type T to an unsigned type.
   *
   * This trait has a single member `type` which is the unsigned type corresponding to T.
   * Note that
   * if T is signed, then member `type` *should define a type with one more bit precision than T*.  For built-in
   * types this trait defaults to `boost::make_unsigned<T>::type`.  For user defined types it simply asserts that
   * the argument type T is an unsigned integer (using std::numeric_limits).
   * User defined specializations may be provided for other cases.
   */
   template <class T>
   struct make_unsigned
   // \cond show_private
      : public make_unsigned_imp < T, boost::is_integral<T>::value >
      // \endcond
   {};
   // \cond show_private
   template <class T, bool intrinsic>
   struct make_unsigned_or_unbounded_imp
   {
      typedef typename boost::make_unsigned<T>::type type;
   };
   template <class T>
   struct make_unsigned_or_unbounded_imp<T, false>
   {
      BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
      BOOST_STATIC_ASSERT((std::numeric_limits<T>::is_signed == false) || (std::numeric_limits<T>::is_bounded == false));
      BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_integer == true);
      typedef T type;
   };
   // \endcond
   /** \brief Converts the argument type T to either an unsigned type or an unbounded integer type.
   *
   * This trait has a single member `type` which is either the unsigned type corresponding to T or an unbounded
   * integer type.  This trait is used to generate types suitable for the calculation of a range: as a result
   * if T is signed, then member `type` *should define a type with one more bit precision than T*.  For built-in
   * types this trait defaults to `boost::make_unsigned<T>::type`.  For user defined types it simply asserts that
   * the argument type T is either an unbounded integer, or an unsigned one (using std::numeric_limits).
   * User defined specializations may be provided for other cases.
   */
   template <class T>
   struct make_unsigned_or_unbounded
      // \cond show_private
      : public make_unsigned_or_unbounded_imp < T, boost::is_integral<T>::value >
      // \endcond
   {};
   /** \brief Traits class that indicates whether type T is an integer
   */
   template <class T>
   struct is_integral
      : public mpl::bool_<boost::is_integral<T>::value || (std::numeric_limits<T>::is_integer)>
   {};
   /** \brief Traits class that indicates whether type T is a signed integer
   */
   template <class T> struct is_signed
      : public mpl::bool_ < boost::is_signed<T>::value || (std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_signed)>
   {};

}
}
}

#endif


code:
/* boost random/triangle_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_TRIANGLE_DISTRIBUTION_HPP
#define BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <iosfwd>
#include <ios>
#include <istream>
#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 {

/**
 * Instantiations of @c triangle_distribution model a \random_distribution.
 * A @c triangle_distribution has three parameters, @c a, @c b, and @c c,
 * which are the smallest, the most probable and the largest values of
 * the distribution respectively.
 */
template<class RealType = double>
class triangle_distribution
{
public:
    typedef RealType input_type;
    typedef RealType result_type;

    class param_type
    {
    public:

        typedef triangle_distribution distribution_type;

        /** Constructs the parameters of a @c triangle_distribution. */
        explicit param_type(RealType a_arg = RealType(0.0),
                            RealType b_arg = RealType(0.5),
                            RealType c_arg = RealType(1.0))
          : _a(a_arg), _b(b_arg), _c(c_arg)
        {
            BOOST_ASSERT(_a <= _b && _b <= _c);
        }

        /** Returns the minimum value of the distribution. */
        RealType a() const { return _a; }
        /** Returns the mode of the distribution. */
        RealType b() const { return _b; }
        /** Returns the maximum value of the distribution. */
        RealType c() const { return _c; }

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

        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            double a_in, b_in, c_in;
            if(is >> a_in >> std::ws >> b_in >> std::ws >> c_in) {
                if(a_in <= b_in && b_in <= c_in) {
                    parm._a = a_in;
                    parm._b = b_in;
                    parm._c = c_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._a == rhs._a && lhs._b == rhs._b && lhs._c == rhs._c; }

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

    private:
        RealType _a;
        RealType _b;
        RealType _c;
    };

    /**
     * Constructs a @c triangle_distribution with the parameters
     * @c a, @c b, and @c c.
     *
     * Preconditions: a <= b <= c.
     */
    explicit triangle_distribution(RealType a_arg = RealType(0.0),
                                   RealType b_arg = RealType(0.5),
                                   RealType c_arg = RealType(1.0))
      : _a(a_arg), _b(b_arg), _c(c_arg)
    {
        BOOST_ASSERT(_a <= _b && _b <= _c);
        init();
    }

    /** Constructs a @c triangle_distribution from its parameters. */
    explicit triangle_distribution(const param_type& parm)
      : _a(parm.a()), _b(parm.b()), _c(parm.c())
    {
        init();
    }

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

    /** Returns the @c a parameter of the distribution */
    result_type a() const { return _a; }
    /** Returns the @c b parameter of the distribution */
    result_type b() const { return _b; }
    /** Returns the @c c parameter of the distribution */
    result_type c() const { return _c; }

    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _a; }
    /** Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _c; }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_a, _b, _c); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _a = parm.a();
        _b = parm.b();
        _c = parm.c();
        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
     * triangle distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng)
    {
        using std::sqrt;
        result_type u = uniform_01<result_type>()(eng);
        if( u <= q1 )
            return _a + p1*sqrt(u);
        else
        return _c - d3*sqrt(d2*u-d1);
    }

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

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

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

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

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

private:
    /// \cond show_private
    void init()
    {
        using std::sqrt;
        d1 = _b - _a;
        d2 = _c - _a;
        d3 = sqrt(_c - _b);
        q1 = d1 / d2;
        p1 = sqrt(d1 * d2);
    }
    /// \endcond

    RealType _a, _b, _c;
    RealType d1, d2, d3, q1, p1;
};

} // namespace random

using random::triangle_distribution;

} // namespace boost

#endif // BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP


code:
/* boost random/uniform_01.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * 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_UNIFORM_01_HPP
#define BOOST_RANDOM_UNIFORM_01_HPP

#include <iostream>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/ptr_helper.hpp>

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

namespace boost {
namespace random {

#ifdef BOOST_RANDOM_DOXYGEN

/**
 * The distribution function uniform_01 models a \random_distribution.
 * On each invocation, it returns a random floating-point value
 * uniformly distributed in the range [0..1).
 *
 * The template parameter RealType shall denote a float-like value type
 * with support for binary operators +, -, and /.
 *
 * Note: The current implementation is buggy, because it may not fill
 * all of the mantissa with random bits. I'm unsure how to fill a
 * (to-be-invented) @c boost::bigfloat class with random bits efficiently.
 * It's probably time for a traits class.
 */
template<class RealType = double>
class uniform_01
{
public:
  typedef RealType input_type;
  typedef RealType result_type;
  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const;
  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const;
  void reset();

  template<class Engine>
  result_type operator()(Engine& eng);

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const new_uniform_01&)
  {
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT,Traits>&
  operator>>(std::basic_istream<CharT,Traits>& is, new_uniform_01&)
  {
    return is;
  }
#endif
};

#else

namespace detail {

template<class RealType>
class new_uniform_01
{
public:
  typedef RealType input_type;
  typedef RealType result_type;
  // compiler-generated copy ctor and copy assignment are fine
  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
  void reset() { }

  template<class Engine>
  result_type operator()(Engine& eng) {
    for (;;) {
      typedef typename Engine::result_type base_result;
      result_type factor = result_type(1) /
              (result_type(base_result((eng.max)()-(eng.min)())) +
               result_type(std::numeric_limits<base_result>::is_integer ? 1 : 0));
      result_type result = result_type(base_result(eng() - (eng.min)())) * factor;
      if (result < result_type(1))
        return result;
    }
  }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const new_uniform_01&)
  {
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT,Traits>&
  operator>>(std::basic_istream<CharT,Traits>& is, new_uniform_01&)
  {
    return is;
  }
#endif
};

template<class UniformRandomNumberGenerator, class RealType>
class backward_compatible_uniform_01
{
  typedef boost::random::detail::ptr_helper<UniformRandomNumberGenerator> traits;
public:
  typedef UniformRandomNumberGenerator base_type;
  typedef RealType result_type;

  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

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

  explicit backward_compatible_uniform_01(typename traits::rvalue_type rng)
    : _rng(rng),
      _factor(result_type(1) /
              (result_type((base().max)()-(base().min)()) +
               result_type(std::numeric_limits<base_result>::is_integer ? 1 : 0)))
  {
  }
  // compiler-generated copy ctor and copy assignment are fine

  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
  typename traits::value_type& base() { return traits::ref(_rng); }
  const typename traits::value_type& base() const { return traits::ref(_rng); }
  void reset() { }

  result_type operator()() {
    for (;;) {
      result_type result = result_type(base()() - (base().min)()) * _factor;
      if (result < result_type(1))
        return result;
    }
  }

#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const backward_compatible_uniform_01& u)
  {
    os << u._rng;
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT,Traits>&
  operator>>(std::basic_istream<CharT,Traits>& is, backward_compatible_uniform_01& u)
  {
    is >> u._rng;
    return is;
  }
#endif

private:
  typedef typename traits::value_type::result_type base_result;
  UniformRandomNumberGenerator _rng;
  result_type _factor;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, class RealType>
const bool backward_compatible_uniform_01<UniformRandomNumberGenerator, RealType>::has_fixed_range;
#endif

template<class UniformRandomNumberGenerator, bool is_number = std::numeric_limits<UniformRandomNumberGenerator>::is_specialized>
struct select_uniform_01
{
  template<class RealType>
  struct apply
  {
    typedef backward_compatible_uniform_01<UniformRandomNumberGenerator, RealType> type;
  };
};

template<class Num>
struct select_uniform_01<Num, true>
{
  template<class RealType>
  struct apply
  {
    typedef new_uniform_01<Num> type;
  };
};

}

// Because it is so commonly used: uniform distribution on the real [0..1)
// range.  This allows for specializations to avoid a costly int -> float
// conversion plus float multiplication
template<class UniformRandomNumberGenerator = double, class RealType = double>
class uniform_01
  : public detail::select_uniform_01<UniformRandomNumberGenerator>::BOOST_NESTED_TEMPLATE apply<RealType>::type
{
  typedef typename detail::select_uniform_01<UniformRandomNumberGenerator>::BOOST_NESTED_TEMPLATE apply<RealType>::type impl_type;
  typedef boost::random::detail::ptr_helper<UniformRandomNumberGenerator> traits;
public:

  uniform_01() {}

  explicit uniform_01(typename traits::rvalue_type rng)
    : impl_type(rng)
  {
  }

#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_01& u)
  {
    os << static_cast<const impl_type&>(u);
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT,Traits>&
  operator>>(std::basic_istream<CharT,Traits>& is, uniform_01& u)
  {
    is >> static_cast<impl_type&>(u);
    return is;
  }
#endif
};

#endif

} // namespace random

using random::uniform_01;

} // namespace boost

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

#endif // BOOST_RANDOM_UNIFORM_01_HPP


code:
/* boost random/uniform_int.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * 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-04-08  added min<max assertion (N. Becker)
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_UNIFORM_INT_HPP
#define BOOST_RANDOM_UNIFORM_INT_HPP

#include <boost/assert.hpp>
#include <boost/random/uniform_int_distribution.hpp>

namespace boost {

/**
 * The distribution function uniform_int models a \random_distribution.
 * On each invocation, it returns a random integer value uniformly
 * distributed in the set of integer numbers {min, min+1, min+2, ..., max}.
 *
 * The template parameter IntType shall denote an integer-like value type.
 *
 * This class is deprecated.  Please use @c uniform_int_distribution in
 * new code.
 */
template<class IntType = int>
class uniform_int : public random::uniform_int_distribution<IntType>
{
    typedef random::uniform_int_distribution<IntType> base_type;
public:

    class param_type : public base_type::param_type
    {
    public:
        typedef uniform_int distribution_type;
        /**
         * Constructs the parameters of a uniform_int distribution.
         *
         * Requires: min <= max
         */
        explicit param_type(IntType min_arg = 0, IntType max_arg = 9)
          : base_type::param_type(min_arg, max_arg)
        {}
    };

    /**
     * Constructs a uniform_int object. @c min and @c max are
     * the parameters of the distribution.
     *
     * Requires: min <= max
     */
    explicit uniform_int(IntType min_arg = 0, IntType max_arg = 9)
      : base_type(min_arg, max_arg)
    {}

    /** Constructs a uniform_int distribution from its parameters. */
    explicit uniform_int(const param_type& parm)
      : base_type(parm)
    {}

    /** Returns the parameters of the distribution */
    param_type param() const { return param_type(this->a(), this->b()); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm) { this->base_type::param(parm); }

    // Codergear seems to have trouble with a using declaration here

    template<class Engine>
    IntType operator()(Engine& eng) const
    {
        return static_cast<const base_type&>(*this)(eng);
    }

    template<class Engine>
    IntType operator()(Engine& eng, const param_type& parm) const
    {
        return static_cast<const base_type&>(*this)(eng, parm);
    }

    template<class Engine>
    IntType operator()(Engine& eng, IntType n) const
    {
        BOOST_ASSERT(n > 0);
        return static_cast<const base_type&>(*this)(eng, param_type(0, n - 1));
    }
};

} // namespace boost

#endif // BOOST_RANDOM_UNIFORM_INT_HPP