code:
/* boost random/linear_feedback_shift.hpp header file
 *
 * Copyright Jens Maurer 2002
 * 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_LINEAR_FEEDBACK_SHIFT_HPP
#define BOOST_RANDOM_LINEAR_FEEDBACK_SHIFT_HPP

#include <iosfwd>
#include <stdexcept>
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/static_assert.hpp>
#include <boost/integer/integer_mask.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>

namespace boost {
namespace random {

/**
 * Instatiations of @c linear_feedback_shift model a
 * \pseudo_random_number_generator.  It was originally
 * proposed in
 *
 *  @blockquote
 *  "Random numbers generated by linear recurrence modulo two.",
 *  Tausworthe, R. C.(1965), Mathematics of Computation 19, 201-209.
 *  @endblockquote
 */
template<class UIntType, int w, int k, int q, int s>
class linear_feedback_shift_engine
{
public:
    typedef UIntType result_type;
    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    BOOST_STATIC_CONSTANT(int, word_size = w);
    BOOST_STATIC_CONSTANT(int, exponent1 = k);
    BOOST_STATIC_CONSTANT(int, exponent2 = q);
    BOOST_STATIC_CONSTANT(int, step_size = s);
    BOOST_STATIC_CONSTANT(UIntType, default_seed = 341);

    /** 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 wordmask(); }

    BOOST_STATIC_ASSERT(w > 0);
    BOOST_STATIC_ASSERT(q > 0);
    BOOST_STATIC_ASSERT(k < w);
    BOOST_STATIC_ASSERT(0 < 2*q && 2*q < k);
    BOOST_STATIC_ASSERT(0 < s && s <= k-q);

    /** Constructs a @c linear_feedback_shift_engine, using the default seed. */
    linear_feedback_shift_engine() { seed(); }

    /** Constructs a @c linear_feedback_shift_engine, seeding it with s0. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_feedback_shift_engine,
        UIntType, s0)
    { seed(s0); }

    /** Constructs a @c linear_feedback_shift_engine, seeding it with seq. */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_feedback_shift_engine,
        SeedSeq, seq)
    { seed(seq); }

    /**
     * Constructs a @c linear_feedback_shift_engine, seeding it with
     * values from the range [first, last).
     */
    template<class It> linear_feedback_shift_engine(It& first, It last)
    { seed(first, last); }

    /** Seeds a @c linear_feedback_shift_engine with the default seed. */
    void seed() {  seed(default_seed); }

    /** Seeds a @c linear_feedback_shift_engine with @c s0. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(linear_feedback_shift_engine,
        UIntType, s0)
    {
        value = s0 & wordmask();
        if(value < (1 << (w-k))) {
            value += 1 << (w-k);
        }
    }

    /**
     * Seeds a @c linear_feedback_shift_engine with values
     * produced by @c seq.generate().
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(linear_feedback_shift_engine,
        SeedSeq, seq)
    { seed(detail::seed_one_int<UIntType, (UIntType(2) << (w - 1))>(seq)); }

    /**
     * Seeds a @c linear_feedback_shift_engine with values
     * from the range [first, last).
     */
    template<class It> void seed(It& first, It last)
    {
        seed(detail::get_one_int<UIntType, (UIntType(2) << (w - 1))>(first, last));
    }

    /** Returns the next value of the generator. */
    result_type operator()()
    {
        const UIntType b = (((value << q) ^ value) & wordmask()) >> (k-s);
        const UIntType mask = (wordmask() << (w-k)) & wordmask();
        value = ((value & mask) << s) ^ b;
        return value;
    }

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

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

    /**
     * Writes the textual representation of the generator to a @c std::ostream.
     */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, linear_feedback_shift_engine, x)
    {
        os << x.value;
        return os;
    }

    /**
     * Reads the textual representation of the generator from a @c std::istream.
     */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, linear_feedback_shift_engine, x)
    {
        is >> x.value;
        return is;
    }

    /**
     * Returns true if the two generators will produce identical
     * sequences of outputs.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(linear_feedback_shift_engine, x, y)
    { return x.value == y.value; }

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

private:
    /// \cond show_private
    static UIntType wordmask() { return boost::low_bits_mask_t<w>::sig_bits; }
    /// \endcond
    UIntType value;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class UIntType, int w, int k, int q, int s>
const bool linear_feedback_shift_engine<UIntType, w, k, q, s>::has_fixed_range;
template<class UIntType, int w, int k, int q, int s>
const int linear_feedback_shift_engine<UIntType, w, k, q, s>::word_size;
template<class UIntType, int w, int k, int q, int s>
const int linear_feedback_shift_engine<UIntType, w, k, q, s>::exponent1;
template<class UIntType, int w, int k, int q, int s>
const int linear_feedback_shift_engine<UIntType, w, k, q, s>::exponent2;
template<class UIntType, int w, int k, int q, int s>
const int linear_feedback_shift_engine<UIntType, w, k, q, s>::step_size;
template<class UIntType, int w, int k, int q, int s>
const UIntType linear_feedback_shift_engine<UIntType, w, k, q, s>::default_seed;
#endif

/// \cond show_deprecated

/** Provided for backwards compatibility. */
template<class UIntType, int w, int k, int q, int s, UIntType v = 0>
class linear_feedback_shift :
    public linear_feedback_shift_engine<UIntType, w, k, q, s>
{
    typedef linear_feedback_shift_engine<UIntType, w, k, q, s> base_type;
public:
    linear_feedback_shift() {}
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_feedback_shift,
        SeedSeq, seq)
    { seed(seq); }
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_feedback_shift,
        UIntType, val)
    { seed(val); }
    template<class It>
    linear_feedback_shift(It& first, It last) : base_type(first, last) {}
};

/// \endcond

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_LINEAR_FEEDBACK_SHIFT_HPP


code:
/* boost random/lognormal_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_LOGNORMAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>      // std::exp, std::sqrt
#include <cassert>
#include <iosfwd>
#include <istream>
#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/normal_distribution.hpp>

namespace boost {
namespace random {

/**
 * Instantiations of class template lognormal_distribution model a
 * \random_distribution. Such a distribution produces random numbers
 * with \f$\displaystyle p(x) = \frac{1}{x s \sqrt{2\pi}} e^{\frac{-\left(\log(x)-m\right)^2}{2s^2}}\f$
 * for x > 0.
 *
 * @xmlwarning
 * This distribution has been updated to match the C++ standard.
 * Its behavior has changed from the original
 * boost::lognormal_distribution.  A backwards compatible
 * version is provided in namespace boost.
 * @endxmlwarning
 */
template<class RealType = double>
class lognormal_distribution
{
public:
    typedef typename normal_distribution<RealType>::input_type input_type;
    typedef RealType result_type;

    class param_type
    {
    public:

        typedef lognormal_distribution distribution_type;

        /** Constructs the parameters of a lognormal_distribution. */
        explicit param_type(RealType m_arg = RealType(0.0),
                            RealType s_arg = RealType(1.0))
          : _m(m_arg), _s(s_arg) {}

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

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

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

        /** Reads the parameters from a std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            is >> parm._m >> std::ws >> parm._s;
            return is;
        }

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

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

    private:
        RealType _m;
        RealType _s;
    };

    /**
     * Constructs a lognormal_distribution. @c m and @c s are the
     * parameters of the distribution.
     */
    explicit lognormal_distribution(RealType m_arg = RealType(0.0),
                                    RealType s_arg = RealType(1.0))
      : _normal(m_arg, s_arg) {}

    /**
     * Constructs a lognormal_distribution from its parameters.
     */
    explicit lognormal_distribution(const param_type& parm)
      : _normal(parm.m(), parm.s()) {}

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

    /** Returns the m parameter of the distribution. */
    RealType m() const { return _normal.mean(); }
    /** Returns the s parameter of the distribution. */
    RealType s() const { return _normal.sigma(); }

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

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(m(), s()); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        typedef normal_distribution<RealType> normal_type;
        typename normal_type::param_type normal_param(parm.m(), parm.s());
        _normal.param(normal_param);
    }

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

    /**
     * Returns a random variate distributed according to the
     * lognormal distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng)
    {
        using std::exp;
        return exp(_normal(eng));
    }

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

    /** Writes the distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lognormal_distribution, ld)
    {
        os << ld._normal;
        return os;
    }

    /** Reads the distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lognormal_distribution, ld)
    {
        is >> ld._normal;
        return is;
    }

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

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

private:
    normal_distribution<result_type> _normal;
};

} // namespace random

/// \cond show_deprecated

/**
 * Provided for backwards compatibility.  This class is
 * deprecated.  It provides the old behavior of lognormal_distribution with
 * \f$\displaystyle p(x) = \frac{1}{x \sigma_N \sqrt{2\pi}} e^{\frac{-\left(\log(x)-\mu_N\right)^2}{2\sigma_N^2}}\f$
 * for x > 0, where \f$\displaystyle \mu_N = \log\left(\frac{\mu^2}{\sqrt{\sigma^2 + \mu^2}}\right)\f$ and
 * \f$\displaystyle \sigma_N = \sqrt{\log\left(1 + \frac{\sigma^2}{\mu^2}\right)}\f$.
 */
template<class RealType = double>
class lognormal_distribution
{
public:
    typedef typename normal_distribution<RealType>::input_type input_type;
    typedef RealType result_type;

    lognormal_distribution(RealType mean_arg = RealType(1.0),
                           RealType sigma_arg = RealType(1.0))
      : _mean(mean_arg), _sigma(sigma_arg)
    {
        init();
    }
    RealType mean() const { return _mean; }
    RealType sigma() const { return _sigma; }
    void reset() { _normal.reset(); }
    template<class Engine>
    RealType operator()(Engine& eng)
    {
        using std::exp;
        return exp(_normal(eng) * _nsigma + _nmean);
    }
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lognormal_distribution, ld)
    {
        os << ld._normal << " " << ld._mean << " " << ld._sigma;
        return os;
    }
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lognormal_distribution, ld)
    {
        is >> ld._normal >> std::ws >> ld._mean >> std::ws >> ld._sigma;
        ld.init();
        return is;
    }
private:
    /// \cond show_private
    void init()
    {
        using std::log;
        using std::sqrt;
        _nmean = log(_mean*_mean/sqrt(_sigma*_sigma + _mean*_mean));
        _nsigma = sqrt(log(_sigma*_sigma/_mean/_mean+result_type(1)));
    }
    RealType _mean;
    RealType _sigma;
    RealType _nmean;
    RealType _nsigma;
    normal_distribution<RealType> _normal;
    /// \endcond
};

/// \endcond

} // namespace boost

#endif // BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP

license: mit


code:
/* boost random/mersenne_twister.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$
 *
 * Revision history
 *  2013-10-14  fixed some warnings with Wshadow (mgaunard)
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
#define BOOST_RANDOM_MERSENNE_TWISTER_HPP

#include <iosfwd>
#include <istream>
#include <stdexcept>
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/integer/integer_mask.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/ptr_helper.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/seed_impl.hpp>
#include <boost/random/detail/generator_seed_seq.hpp>
#include <boost/random/detail/polynomial.hpp>

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

namespace boost {
namespace random {

/**
 * Instantiations of class template mersenne_twister_engine model a
 * \pseudo_random_number_generator. It uses the algorithm described in
 *
 *  @blockquote
 *  "Mersenne Twister: A 623-dimensionally equidistributed uniform
 *  pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura,
 *  ACM Transactions on Modeling and Computer Simulation: Special Issue on
 *  Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
 *  @endblockquote
 *
 * @xmlnote
 * The boost variant has been implemented from scratch and does not
 * derive from or use mt19937.c provided on the above WWW site. However, it
 * was verified that both produce identical output.
 * @endxmlnote
 *
 * The seeding from an integer was changed in April 2005 to address a
 * weakness.
 *
 * The quality of the generator crucially depends on the choice of the
 * parameters.  User code should employ one of the sensibly parameterized
 * generators such as \mt19937 instead.
 *
 * The generator requires considerable amounts of memory for the storage of
 * its state array. For example, \mt11213b requires about 1408 bytes and
 * \mt19937 requires about 2496 bytes.
 */
template<class UIntType,
         std::size_t w, std::size_t n, std::size_t m, std::size_t r,
         UIntType a, std::size_t u, UIntType d, std::size_t s,
         UIntType b, std::size_t t,
         UIntType c, std::size_t l, UIntType f>
class mersenne_twister_engine
{
public:
    typedef UIntType result_type;
    BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
    BOOST_STATIC_CONSTANT(std::size_t, state_size = n);
    BOOST_STATIC_CONSTANT(std::size_t, shift_size = m);
    BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r);
    BOOST_STATIC_CONSTANT(UIntType, xor_mask = a);
    BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u);
    BOOST_STATIC_CONSTANT(UIntType, tempering_d = d);
    BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s);
    BOOST_STATIC_CONSTANT(UIntType, tempering_b = b);
    BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t);
    BOOST_STATIC_CONSTANT(UIntType, tempering_c = c);
    BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l);
    BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f);
    BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u);

    // backwards compatibility
    BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
    BOOST_STATIC_CONSTANT(std::size_t, output_u = u);
    BOOST_STATIC_CONSTANT(std::size_t, output_s = s);
    BOOST_STATIC_CONSTANT(UIntType, output_b = b);
    BOOST_STATIC_CONSTANT(std::size_t, output_t = t);
    BOOST_STATIC_CONSTANT(UIntType, output_c = c);
    BOOST_STATIC_CONSTANT(std::size_t, output_l = l);

    // old Boost.Random concept requirements
    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

    /**
     * Constructs a @c mersenne_twister_engine and calls @c seed().
     */
    mersenne_twister_engine() { seed(); }

    /**
     * Constructs a @c mersenne_twister_engine and calls @c seed(value).
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine,
                                               UIntType, value)
    { seed(value); }
    template<class It> mersenne_twister_engine(It& first, It last)
    { seed(first,last); }

    /**
     * Constructs a mersenne_twister_engine and calls @c seed(gen).
     *
     * @xmlnote
     * The copy constructor will always be preferred over
     * the templated constructor.
     * @endxmlnote
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine,
                                             SeedSeq, seq)
    { seed(seq); }

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

    /** Calls @c seed(default_seed). */
    void seed() { seed(default_seed); }

    /**
     * Sets the state x(0) to v mod 2w. Then, iteratively,
     * sets x(i) to
     * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w
     * for i = 1 .. n-1. x(n) is the first value to be returned by operator().
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value)
    {
        // New seeding algorithm from
        // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
        // In the previous versions, MSBs of the seed affected only MSBs of the
        // state x[].
        const UIntType mask = (max)();
        x[0] = value & mask;
        for (i = 1; i < n; i++) {
            // See Knuth "The Art of Computer Programming"
            // Vol. 2, 3rd ed., page 106
            x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
        }

        normalize_state();
    }

    /**
     * Seeds a mersenne_twister_engine using values produced by seq.generate().
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq)
    {
        detail::seed_array_int<w>(seq, x);
        i = n;

        normalize_state();
    }

    /** Sets the state of the generator using values from an iterator range. */
    template<class It>
    void seed(It& first, It last)
    {
        detail::fill_array_int<w>(first, last, x);
        i = n;

        normalize_state();
    }

    /** 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; }

    /** Produces the next value of the generator. */
    result_type operator()();

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

    /**
     * Advances the state of the generator by @c z steps.  Equivalent to
     *
     * @code
     * for(unsigned long long i = 0; i < z; ++i) {
     *     gen();
     * }
     * @endcode
     */
    void discard(boost::uintmax_t z)
    {
#ifndef BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD
#define BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD 10000000
#endif
        if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
            discard_many(z);
        } else {
            for(boost::uintmax_t j = 0; j < z; ++j) {
                (*this)();
            }
        }
    }

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

    /** Reads a mersenne_twister_engine from a @c std::istream */
    template<class CharT, class Traits>
    friend std::basic_istream<CharT,Traits>&
    operator>>(std::basic_istream<CharT,Traits>& is,
               mersenne_twister_engine& mt)
    {
        for(std::size_t j = 0; j < mt.state_size; ++j)
            is >> mt.x[j] >> std::ws;
        // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template
        // value parameter "n" available from the class template scope, so use
        // the static constant with the same value
        mt.i = mt.state_size;
        return is;
    }
#endif

    /**
     * Returns true if the two generators are in the same state,
     * and will thus produce identical sequences.
     */
    friend bool operator==(const mersenne_twister_engine& x_,
                           const mersenne_twister_engine& y_)
    {
        if(x_.i < y_.i) return x_.equal_imp(y_);
        else return y_.equal_imp(x_);
    }

    /**
     * Returns true if the two generators are in different states.
     */
    friend bool operator!=(const mersenne_twister_engine& x_,
                           const mersenne_twister_engine& y_)
    { return !(x_ == y_); }

private:
    /// \cond show_private

    void twist();

    /**
     * Does the work of operator==.  This is in a member function
     * for portability.  Some compilers, such as msvc 7.1 and
     * Sun CC 5.10 can't access template parameters or static
     * members of the class from inline friend functions.
     *
     * requires i <= other.i
     */
    bool equal_imp(const mersenne_twister_engine& other) const
    {
        UIntType back[n];
        std::size_t offset = other.i - i;
        for(std::size_t j = 0; j + offset < n; ++j)
            if(x[j] != other.x[j+offset])
                return false;
        rewind(&back[n-1], offset);
        for(std::size_t j = 0; j < offset; ++j)
            if(back[j + n - offset] != other.x[j])
                return false;
        return true;
    }

    /**
     * Does the work of operator<<.  This is in a member function
     * for portability.
     */
    template<class CharT, class Traits>
    void print(std::basic_ostream<CharT, Traits>& os) const
    {
        UIntType data[n];
        for(std::size_t j = 0; j < i; ++j) {
            data[j + n - i] = x[j];
        }
        if(i != n) {
            rewind(&data[n - i - 1], n - i);
        }
        os << data[0];
        for(std::size_t j = 1; j < n; ++j) {
            os << ' ' << data[j];
        }
    }

    /**
     * Copies z elements of the state preceding x[0] into
     * the array whose last element is last.
     */
    void rewind(UIntType* last, std::size_t z) const
    {
        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
        const UIntType lower_mask = ~upper_mask;
        UIntType y0 = x[m-1] ^ x[n-1];
        if(y0 & (static_cast<UIntType>(1) << (w-1))) {
            y0 = ((y0 ^ a) << 1) | 1;
        } else {
            y0 = y0 << 1;
        }
        for(std::size_t sz = 0; sz < z; ++sz) {
            UIntType y1 =
                rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1);
            if(y1 & (static_cast<UIntType>(1) << (w-1))) {
                y1 = ((y1 ^ a) << 1) | 1;
            } else {
                y1 = y1 << 1;
            }
            *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask);
            y0 = y1;
        }
    }

    /**
     * Converts an arbitrary array into a valid generator state.
     * First we normalize x[0], so that it contains the same
     * value we would get by running the generator forwards
     * and then in reverse.  (The low order r bits are redundant).
     * Then, if the state consists of all zeros, we set the
     * high order bit of x[0] to 1.  This function only needs to
     * be called by seed, since the state transform preserves
     * this relationship.
     */
    void normalize_state()
    {
        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
        const UIntType lower_mask = ~upper_mask;
        UIntType y0 = x[m-1] ^ x[n-1];
        if(y0 & (static_cast<UIntType>(1) << (w-1))) {
            y0 = ((y0 ^ a) << 1) | 1;
        } else {
            y0 = y0 << 1;
        }
        x[0] = (x[0] & upper_mask) | (y0 & lower_mask);

        // fix up the state if it's all zeroes.
        for(std::size_t j = 0; j < n; ++j) {
            if(x[j] != 0) return;
        }
        x[0] = static_cast<UIntType>(1) << (w-1);
    }

    /**
     * Given a pointer to the last element of the rewind array,
     * and the current size of the rewind array, finds an element
     * relative to the next available slot in the rewind array.
     */
    UIntType
    rewind_find(UIntType* last, std::size_t size, std::size_t j) const
    {
        std::size_t index = (j + n - size + n - 1) % n;
        if(index < n - size) {
            return x[index];
        } else {
            return *(last - (n - 1 - index));
        }
    }

    /**
     * Optimized algorithm for large jumps.
     *
     * Hiroshi Haramoto, Makoto Matsumoto, and Pierre L'Ecuyer. 2008.
     * A Fast Jump Ahead Algorithm for Linear Recurrences in a Polynomial
     * Space. In Proceedings of the 5th international conference on
     * Sequences and Their Applications (SETA '08).
     * DOI=10.1007/978-3-540-85912-3_26
     */
    void discard_many(boost::uintmax_t z)
    {
        // Compute the minimal polynomial, phi(t)
        // This depends only on the transition function,
        // which is constant.  The characteristic
        // polynomial is the same as the minimal
        // polynomial for a maximum period generator
        // (which should be all specializations of
        // mersenne_twister.)  Even if it weren't,
        // the characteristic polynomial is guaranteed
        // to be a multiple of the minimal polynomial,
        // which is good enough.
        detail::polynomial phi = get_characteristic_polynomial();

        // calculate g(t) = t^z % phi(t)
        detail::polynomial g = mod_pow_x(z, phi);

        // h(s_0, t) = \sum_{i=0}^{2k-1}o(s_i)t^{2k-i-1}
        detail::polynomial h;
        const std::size_t num_bits = w*n - r;
        for(std::size_t j = 0; j < num_bits * 2; ++j) {
            // Yes, we're advancing the generator state
            // here, but it doesn't matter because
            // we're going to overwrite it completely
            // in reconstruct_state.
            if(i >= n) twist();
            h[2*num_bits - j - 1] = x[i++] & UIntType(1);
        }
        // g(t)h(s_0, t)
        detail::polynomial gh = g * h;
        detail::polynomial result;
        for(std::size_t j = 0; j <= num_bits; ++j) {
            result[j] = gh[2*num_bits - j - 1];
        }
        reconstruct_state(result);
    }
    static detail::polynomial get_characteristic_polynomial()
    {
        const std::size_t num_bits = w*n - r;
        detail::polynomial helper;
        helper[num_bits - 1] = 1;
        mersenne_twister_engine tmp;
        tmp.reconstruct_state(helper);
        // Skip the first num_bits elements, since we
        // already know what they are.
        for(std::size_t j = 0; j < num_bits; ++j) {
            if(tmp.i >= n) tmp.twist();
            if(j == num_bits - 1)
                assert((tmp.x[tmp.i] & 1) == 1);
            else
                assert((tmp.x[tmp.i] & 1) == 0);
            ++tmp.i;
        }
        detail::polynomial phi;
        phi[num_bits] = 1;
        detail::polynomial next_bits = tmp.as_polynomial(num_bits);
        for(std::size_t j = 0; j < num_bits; ++j) {
            int val = next_bits[j] ^ phi[num_bits-j-1];
            phi[num_bits-j-1] = val;
            if(val) {
                for(std::size_t k = j + 1; k < num_bits; ++k) {
                    phi[num_bits-k-1] ^= next_bits[k-j-1];
                }
            }
        }
        return phi;
    }
    detail::polynomial as_polynomial(std::size_t size) {
        detail::polynomial result;
        for(std::size_t j = 0; j < size; ++j) {
            if(i >= n) twist();
            result[j] = x[i++] & UIntType(1);
        }
        return result;
    }
    void reconstruct_state(const detail::polynomial& p)
    {
        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
        const UIntType lower_mask = ~upper_mask;
        const std::size_t num_bits = w*n - r;
        for(std::size_t j = num_bits - n + 1; j <= num_bits; ++j)
            x[j % n] = p[j];

        UIntType y0 = 0;
        for(std::size_t j = num_bits + 1; j >= n - 1; --j) {
            UIntType y1 = x[j % n] ^ x[(j + m) % n];
            if(p[j - n + 1])
                y1 = (y1 ^ a) << UIntType(1) | UIntType(1);
            else
                y1 = y1 << UIntType(1);
            x[(j + 1) % n] = (y0 & upper_mask) | (y1 & lower_mask);
            y0 = y1;
        }
        i = 0;
    }

    /// \endcond

    // state representation: next output is o(x(i))
    //   x[0]  ... x[k] x[k+1] ... x[n-1]   represents
    //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1)

    UIntType x[n];
    std::size_t i;
};

/// \cond show_private

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
#define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name)                         \
template<class UIntType, std::size_t w, std::size_t n, std::size_t m,       \
    std::size_t r, UIntType a, std::size_t u, UIntType d, std::size_t s,    \
    UIntType b, std::size_t t, UIntType c, std::size_t l, UIntType f>       \
const type mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::name
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u );
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t);
BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c);
BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l);
BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range);
#undef BOOST_RANDOM_MT_DEFINE_CONSTANT
#endif

template<class UIntType,
         std::size_t w, std::size_t n, std::size_t m, std::size_t r,
         UIntType a, std::size_t u, UIntType d, std::size_t s,
         UIntType b, std::size_t t,
         UIntType c, std::size_t l, UIntType f>
void
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist()
{
    const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
    const UIntType lower_mask = ~upper_mask;

    const std::size_t unroll_factor = 6;
    const std::size_t unroll_extra1 = (n-m) % unroll_factor;
    const std::size_t unroll_extra2 = (m-1) % unroll_factor;

    // split loop to avoid costly modulo operations
    {  // extra scope for MSVC brokenness w.r.t. for scope
        for(std::size_t j = 0; j < n-m-unroll_extra1; j++) {
            UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
            x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
        }
    }
    {
        for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) {
            UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
            x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
        }
    }
    {
        for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) {
            UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
            x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
        }
    }
    {
        for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) {
            UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
            x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
        }
    }
    // last iteration
    UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
    x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a);
    i = 0;
}
/// \endcond

template<class UIntType,
         std::size_t w, std::size_t n, std::size_t m, std::size_t r,
         UIntType a, std::size_t u, UIntType d, std::size_t s,
         UIntType b, std::size_t t,
         UIntType c, std::size_t l, UIntType f>
inline typename
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::operator()()
{
    if(i == n)
        twist();
    // Step 4
    UIntType z = x[i];
    ++i;
    z ^= ((z >> u) & d);
    z ^= ((z << s) & b);
    z ^= ((z << t) & c);
    z ^= (z >> l);
    return z;
}

/**
 * The specializations \mt11213b and \mt19937 are from
 *
 *  @blockquote
 *  "Mersenne Twister: A 623-dimensionally equidistributed
 *  uniform pseudo-random number generator", Makoto Matsumoto
 *  and Takuji Nishimura, ACM Transactions on Modeling and
 *  Computer Simulation: Special Issue on Uniform Random Number
 *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
 *  @endblockquote
 */
typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7,
    11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b;

/**
 * The specializations \mt11213b and \mt19937 are from
 *
 *  @blockquote
 *  "Mersenne Twister: A 623-dimensionally equidistributed
 *  uniform pseudo-random number generator", Makoto Matsumoto
 *  and Takuji Nishimura, ACM Transactions on Modeling and
 *  Computer Simulation: Special Issue on Uniform Random Number
 *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
 *  @endblockquote
 */
typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df,
    11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937;

#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
typedef mersenne_twister_engine<uint64_t,64,312,156,31,
    UINT64_C(0xb5026f5aa96619e9),29,UINT64_C(0x5555555555555555),17,
    UINT64_C(0x71d67fffeda60000),37,UINT64_C(0xfff7eee000000000),43,
    UINT64_C(6364136223846793005)> mt19937_64;
#endif

/// \cond show_deprecated

template<class UIntType,
         int w, int n, int m, int r,
         UIntType a, int u, std::size_t s,
         UIntType b, int t,
         UIntType c, int l, UIntType v>
class mersenne_twister :
    public mersenne_twister_engine<UIntType,
        w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253>
{
    typedef mersenne_twister_engine<UIntType,
        w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> base_type;
public:
    mersenne_twister() {}
    BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen)
    { seed(gen); }
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val)
    { seed(val); }
    template<class It>
    mersenne_twister(It& first, It last) : base_type(first, last) {}
    void seed() { base_type::seed(); }
    BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen)
    {
        detail::generator_seed_seq<Gen> seq(gen);
        base_type::seed(seq);
    }
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val)
    { base_type::seed(val); }
    template<class It>
    void seed(It& first, It last) { base_type::seed(first, last); }
};

/// \endcond

} // namespace random

using random::mt11213b;
using random::mt19937;
using random::mt19937_64;

} // namespace boost

BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)

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

#endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP


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

#ifndef BOOST_RANDOM_NEGATIVE_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_NEGATIVE_BINOMIAL_DISTRIBUTION_HPP_INCLUDED

#include <iosfwd>

#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/poisson_distribution.hpp>

namespace boost {
namespace random {

/**
 * The negative binomial distribution is an integer valued
 * distribution with two parameters, @c k and @c p.  The
 * distribution produces non-negative values.
 *
 * The distribution function is
 * \f$\displaystyle P(i) = {k+i-1\choose i}p^k(1-p)^i\f$.
 *
 * This implementation uses a gamma-poisson mixture.
 */
template<class IntType = int, class RealType = double>
class negative_binomial_distribution {
public:
    typedef IntType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef negative_binomial_distribution distribution_type;
        /**
         * Construct a param_type object.  @c k and @c p
         * are the parameters of the distribution.
         *
         * Requires: k >=0 && 0 <= p <= 1
         */
        explicit param_type(IntType k_arg = 1, RealType p_arg = RealType (0.5))
          : _k(k_arg), _p(p_arg)
        {}
        /** Returns the @c k parameter of the distribution. */
        IntType k() const { return _k; }
        /** Returns the @c p parameter of the distribution. */
        RealType p() const { return _p; }
#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._p << " " << parm._k;
            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._p >> std::ws >> parm._k;
            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._k == rhs._k && lhs._p == rhs._p;
        }
        /** Returns true if the parameters have different values. */
        friend bool operator!=(const param_type& lhs, const param_type& rhs)
        {
            return !(lhs == rhs);
        }
    private:
        IntType _k;
        RealType _p;
    };

    /**
     * Construct a @c negative_binomial_distribution object. @c k and @c p
     * are the parameters of the distribution.
     *
     * Requires: k >=0 && 0 <= p <= 1
     */
    explicit negative_binomial_distribution(IntType k_arg = 1,
                                            RealType p_arg = RealType(0.5))
      : _k(k_arg), _p(p_arg)
    {}

    /**
     * Construct an @c negative_binomial_distribution object from the
     * parameters.
     */
    explicit negative_binomial_distribution(const param_type& parm)
      : _k(parm.k()), _p(parm.p())
    {}

    /**
     * Returns a random variate distributed according to the
     * negative binomial distribution.
     */
    template<class URNG>
    IntType operator()(URNG& urng) const
    {
        gamma_distribution<RealType> gamma(_k, (1-_p)/_p);
        poisson_distribution<IntType, RealType> poisson(gamma(urng));
        return poisson(urng);
    }

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

    /** Returns the @c k parameter of the distribution. */
    IntType k() const { return _k; }
    /** Returns the @c p parameter of the distribution. */
    RealType p() const { return _p; }

    /** 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(_k, _p); }
    /** Sets parameters of the distribution. */
    void param(const param_type& parm)
    {
        _k = parm.k();
        _p = parm.p();
    }

    /**
     * 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 negative_binomial_distribution& bd)
    {
        os << bd.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,
               negative_binomial_distribution& bd)
    {
        bd.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 negative_binomial_distribution& lhs,
                           const negative_binomial_distribution& rhs)
    {
        return lhs._k == rhs._k && lhs._p == rhs._p;
    }
    /** Returns true if the two distributions could produce different
        sequences of values, given equal generators. */
    friend bool operator!=(const negative_binomial_distribution& lhs,
                           const negative_binomial_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);
        }
    }

    // parameters
    IntType _k;
    RealType _p;

    /// @endcond
};

}

}

#endif


code:
/* boost random/non_central_chi_squared_distribution.hpp header file
 *
 * Copyright Thijs van den Berg 2014
 *
 * 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_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP
#define BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP

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

namespace boost {
namespace random {

/**
 * The noncentral chi-squared distribution is a real valued distribution with
 * two parameter, @c k and @c lambda.  The distribution produces values > 0.
 *
 * This is the distribution of the sum of squares of k Normal distributed
 * variates each with variance one and \f$\lambda\f$ the sum of squares of the
 * normal means.
 *
 * The distribution function is
 * \f$\displaystyle P(x) = \frac{1}{2} e^{-(x+\lambda)/2} \left( \frac{x}{\lambda} \right)^{k/4-1/2} I_{k/2-1}( \sqrt{\lambda x} )\f$.
 *  where  \f$\displaystyle I_\nu(z)\f$ is a modified Bessel function of the
 * first kind.
 *
 * The algorithm is taken from
 *
 *  @blockquote
 *  "Monte Carlo Methods in Financial Engineering", Paul Glasserman,
 *  2003, XIII, 596 p, Stochastic Modelling and Applied Probability, Vol. 53,
 *  ISBN 978-0-387-21617-1, p 124, Fig. 3.5.
 *  @endblockquote
 */
template <typename RealType = double>
class non_central_chi_squared_distribution {
public:
    typedef RealType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef non_central_chi_squared_distribution distribution_type;

        /**
         * Constructs the parameters of a non_central_chi_squared_distribution.
         * @c k and @c lambda are the parameter of the distribution.
         *
         * Requires: k > 0 && lambda > 0
         */
        explicit
        param_type(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
        : _k(k_arg), _lambda(lambda_arg)
        {
            BOOST_ASSERT(k_arg > RealType(0));
            BOOST_ASSERT(lambda_arg > RealType(0));
        }

        /** Returns the @c k parameter of the distribution */
        RealType k() const { return _k; }

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

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

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

        /** Returns true if the parameters have the same values. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._k == rhs._k && lhs._lambda == rhs._lambda; }

        /** Returns true if the parameters have different values. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _k;
        RealType _lambda;
    };

    /**
     * Construct a @c non_central_chi_squared_distribution object. @c k and
     * @c lambda are the parameter of the distribution.
     *
     * Requires: k > 0 && lambda > 0
     */
    explicit
    non_central_chi_squared_distribution(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
      : _param(k_arg, lambda_arg)
    {
        BOOST_ASSERT(k_arg > RealType(0));
        BOOST_ASSERT(lambda_arg > RealType(0));
    }

    /**
     * Construct a @c non_central_chi_squared_distribution object from the parameter.
     */
    explicit
    non_central_chi_squared_distribution(const param_type& parm)
      : _param( parm )
    { }

    /**
     * Returns a random variate distributed according to the
     * non central chi squared distribution specified by @c param.
     */
    template<typename URNG>
    RealType operator()(URNG& eng, const param_type& parm) const
    { return non_central_chi_squared_distribution(parm)(eng); }

    /**
     * Returns a random variate distributed according to the
     * non central chi squared distribution.
     */
    template<typename URNG>
    RealType operator()(URNG& eng)
    {
        using std::sqrt;
        if (_param.k() > 1) {
            boost::random::normal_distribution<RealType> n_dist;
            boost::random::chi_squared_distribution<RealType> c_dist(_param.k() - RealType(1));
            RealType _z = n_dist(eng);
            RealType _x = c_dist(eng);
            RealType term1 = _z + sqrt(_param.lambda());
            return term1*term1 + _x;
        }
        else {
            boost::random::poisson_distribution<> p_dist(_param.lambda()/RealType(2));
            boost::random::poisson_distribution<>::result_type _p = p_dist(eng);
            boost::random::chi_squared_distribution<RealType> c_dist(_param.k() + RealType(2)*_p);
            return c_dist(eng);
        }
    }

    /** Returns the @c k parameter of the distribution. */
    RealType k() const { return _param.k(); }

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

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

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

    /** Resets the distribution, so that subsequent uses does not depend on values already produced by it.*/
    void reset() {}

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

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

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

    /** reads the parameters of the distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, non_central_chi_squared_distribution, dist)
    {
        param_type parm;
        if(is >> parm) {
            dist.param(parm);
        }
        return is;
    }

    /** Returns true if two distributions have the same parameters and produce
        the same sequence of random numbers given equal generators.*/
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(non_central_chi_squared_distribution, lhs, rhs)
    { return lhs.param() == rhs.param(); }

    /** Returns true if two distributions have different parameters and/or can produce
       different sequences of random numbers given equal generators.*/
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(non_central_chi_squared_distribution)

private:

    /// @cond show_private
    param_type  _param;
    /// @endcond
};

} // namespace random
} // namespace boost

#endif


code:
/* boost random/normal_distribution.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2010-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_NORMAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <istream>
#include <iosfwd>
#include <boost/assert.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/integer.hpp>
#include <boost/integer/integer_mask.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/type_traits/make_unsigned.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/integer_log2.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/exponential_distribution.hpp>
#include <boost/mpl/bool.hpp>

namespace boost {
namespace random {

namespace detail {

// tables for the ziggurat algorithm
template<class RealType>
struct normal_table {
    static const RealType table_x[129];
    static const RealType table_y[129];
};

template<class RealType>
const RealType normal_table<RealType>::table_x[129] = {
    3.7130862467403632609, 3.4426198558966521214, 3.2230849845786185446, 3.0832288582142137009,
    2.9786962526450169606, 2.8943440070186706210, 2.8231253505459664379, 2.7611693723841538514,
    2.7061135731187223371, 2.6564064112581924999, 2.6109722484286132035, 2.5690336259216391328,
    2.5300096723854666170, 2.4934545220919507609, 2.4590181774083500943, 2.4264206455302115930,
    2.3954342780074673425, 2.3658713701139875435, 2.3375752413355307354, 2.3104136836950021558,
    2.2842740596736568056, 2.2590595738653295251, 2.2346863955870569803, 2.2110814088747278106,
    2.1881804320720206093, 2.1659267937448407377, 2.1442701823562613518, 2.1231657086697899595,
    2.1025731351849988838, 2.0824562379877246441, 2.0627822745039633575, 2.0435215366506694976,
    2.0246469733729338782, 2.0061338699589668403, 1.9879595741230607243, 1.9701032608497132242,
    1.9525457295488889058, 1.9352692282919002011, 1.9182573008597320303, 1.9014946531003176140,
    1.8849670357028692380, 1.8686611409895420085, 1.8525645117230870617, 1.8366654602533840447,
    1.8209529965910050740, 1.8054167642140487420, 1.7900469825946189862, 1.7748343955807692457,
    1.7597702248942318749, 1.7448461281083765085, 1.7300541605582435350, 1.7153867407081165482,
    1.7008366185643009437, 1.6863968467734863258, 1.6720607540918522072, 1.6578219209482075462,
    1.6436741568569826489, 1.6296114794646783962, 1.6156280950371329644, 1.6017183802152770587,
    1.5878768648844007019, 1.5740982160167497219, 1.5603772223598406870, 1.5467087798535034608,
    1.5330878776675560787, 1.5195095847593707806, 1.5059690368565502602, 1.4924614237746154081,
    1.4789819769830978546, 1.4655259573357946276, 1.4520886428822164926, 1.4386653166774613138,
    1.4252512545068615734, 1.4118417124397602509, 1.3984319141236063517, 1.3850170377251486449,
    1.3715922024197322698, 1.3581524543224228739, 1.3446927517457130432, 1.3312079496576765017,
    1.3176927832013429910, 1.3041418501204215390, 1.2905495919178731508, 1.2769102735516997175,
    1.2632179614460282310, 1.2494664995643337480, 1.2356494832544811749, 1.2217602305309625678,
    1.2077917504067576028, 1.1937367078237721994, 1.1795873846544607035, 1.1653356361550469083,
    1.1509728421389760651, 1.1364898520030755352, 1.1218769225722540661, 1.1071236475235353980,
    1.0922188768965537614, 1.0771506248819376573, 1.0619059636836193998, 1.0464709007525802629,
    1.0308302360564555907, 1.0149673952392994716, 0.99886423348064351303, 0.98250080350276038481,
    0.96585507938813059489, 0.94890262549791195381, 0.93161619660135381056, 0.91396525100880177644,
    0.89591535256623852894, 0.87742742909771569142, 0.85845684317805086354, 0.83895221428120745572,
    0.81885390668331772331, 0.79809206062627480454, 0.77658398787614838598, 0.75423066443451007146,
    0.73091191062188128150, 0.70647961131360803456, 0.68074791864590421664, 0.65347863871504238702,
    0.62435859730908822111, 0.59296294244197797913, 0.55869217837551797140, 0.52065603872514491759,
    0.47743783725378787681, 0.42654798630330512490, 0.36287143102841830424, 0.27232086470466385065,
    0
};

template<class RealType>
const RealType normal_table<RealType>::table_y[129] = {
    0, 0.0026696290839025035092, 0.0055489952208164705392, 0.0086244844129304709682,
    0.011839478657982313715, 0.015167298010672042468, 0.018592102737165812650, 0.022103304616111592615,
    0.025693291936149616572, 0.029356317440253829618, 0.033087886146505155566, 0.036884388786968774128,
    0.040742868074790604632, 0.044660862200872429800, 0.048636295860284051878, 0.052667401903503169793,
    0.056752663481538584188, 0.060890770348566375972, 0.065080585213631873753, 0.069321117394180252601,
    0.073611501884754893389, 0.077950982514654714188, 0.082338898242957408243, 0.086774671895542968998,
    0.091257800827634710201, 0.09578784912257815216, 0.10036444102954554013, 0.10498725541035453978,
    0.10965602101581776100, 0.11437051244988827452, 0.11913054670871858767, 0.12393598020398174246,
    0.12878670619710396109, 0.13368265258464764118, 0.13862377998585103702, 0.14361008009193299469,
    0.14864157424369696566, 0.15371831220958657066, 0.15884037114093507813, 0.16400785468492774791,
    0.16922089223892475176, 0.17447963833240232295, 0.17978427212496211424, 0.18513499701071343216,
    0.19053204032091372112, 0.19597565311811041399, 0.20146611007620324118, 0.20700370944187380064,
    0.21258877307373610060, 0.21822164655637059599, 0.22390269938713388747, 0.22963232523430270355,
    0.23541094226572765600, 0.24123899354775131610, 0.24711694751469673582, 0.25304529850976585934,
    0.25902456739871074263, 0.26505530225816194029, 0.27113807914102527343, 0.27727350292189771153,
    0.28346220822601251779, 0.28970486044581049771, 0.29600215684985583659, 0.30235482778947976274,
    0.30876363800925192282, 0.31522938806815752222, 0.32175291587920862031, 0.32833509837615239609,
    0.33497685331697116147, 0.34167914123501368412, 0.34844296754987246935, 0.35526938485154714435,
    0.36215949537303321162, 0.36911445366827513952, 0.37613546951445442947, 0.38322381105988364587,
    0.39038080824138948916, 0.39760785649804255208, 0.40490642081148835099, 0.41227804010702462062,
    0.41972433205403823467, 0.42724699830956239880, 0.43484783025466189638, 0.44252871528024661483,
    0.45029164368692696086, 0.45813871627287196483, 0.46607215269457097924, 0.47409430069824960453,
    0.48220764633483869062, 0.49041482528932163741, 0.49871863547658432422, 0.50712205108130458951,
    0.51562823824987205196, 0.52424057267899279809, 0.53296265938998758838, 0.54179835503172412311,
    0.55075179312105527738, 0.55982741271069481791, 0.56902999107472161225, 0.57836468112670231279,
    0.58783705444182052571, 0.59745315095181228217, 0.60721953663260488551, 0.61714337082656248870,
    0.62723248525781456578, 0.63749547734314487428, 0.64794182111855080873, 0.65858200005865368016,
    0.66942766735770616891, 0.68049184100641433355, 0.69178914344603585279, 0.70333609902581741633,
    0.71515150742047704368, 0.72725691835450587793, 0.73967724368333814856, 0.75244155918570380145,
    0.76558417390923599480, 0.77914608594170316563, 0.79317701178385921053, 0.80773829469612111340,
    0.82290721139526200050, 0.83878360531064722379, 0.85550060788506428418, 0.87324304892685358879,
    0.89228165080230272301, 0.91304364799203805999, 0.93628268170837107547, 0.96359969315576759960,
    1
};

template<class Engine>
inline typename boost::make_unsigned<typename Engine::result_type>::type
generate_one_digit(Engine& eng, std::size_t bits)
{
    typedef typename Engine::result_type base_result;
    typedef typename boost::make_unsigned<base_result>::type base_unsigned;

    base_unsigned range =
        detail::subtract<base_result>()((eng.max)(), (eng.min)());
    base_unsigned y0_mask = (base_unsigned(2) << (bits - 1)) - 1;
    base_unsigned y0 = (range + 1) & ~y0_mask;
    base_unsigned u;
    do {
        u = detail::subtract<base_result>()(eng(), (eng.min)());
    } while(y0 != 0 && u > base_unsigned(y0 - 1));
    return u & y0_mask;
}

template<class RealType, std::size_t w, class Engine>
std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::true_)
{
    typedef typename Engine::result_type base_result;
    typedef typename boost::make_unsigned<base_result>::type base_unsigned;

    base_unsigned range =
        detail::subtract<base_result>()((eng.max)(), (eng.min)());

    std::size_t m =
        (range == (std::numeric_limits<base_unsigned>::max)()) ?
            std::numeric_limits<base_unsigned>::digits :
            detail::integer_log2(range + 1);

    int bucket = 0;
    // process as many full digits as possible into the int part
    for(std::size_t i = 0; i < w/m; ++i) {
        base_unsigned u = generate_one_digit(eng, m);
        bucket = (bucket << m) | u;
    }
    RealType r;

    const std::size_t digits = std::numeric_limits<RealType>::digits;
    {
        base_unsigned u = generate_one_digit(eng, m);
        base_unsigned mask = (base_unsigned(1) << (w%m)) - 1;
        bucket = (bucket << (w%m)) | (mask & u);
        const RealType mult = RealType(1)/RealType(base_unsigned(1) << (m - w%m));
        // zero out unused bits
        if (m - w%m > digits) {
            u &= ~(base_unsigned(1) << (m - digits));
        }
        r = RealType(u >> (w%m)) * mult;
    }
    for(std::size_t i = m - w%m; i + m < digits; ++i) {
        base_unsigned u = generate_one_digit(eng, m);
        r += u;
        r *= RealType(0.5)/RealType(base_unsigned(1) << (m - 1));
    }
    if (m - w%m < digits)
    {
        const std::size_t remaining = (digits - m + w%m) % m;
        base_unsigned u = generate_one_digit(eng, m);
        r += u & ((base_unsigned(2) << (remaining - 1)) - 1);
        const RealType mult = RealType(0.5)/RealType(base_unsigned(1) << (remaining - 1));
        r *= mult;
    }
    return std::make_pair(r, bucket);
}

template<class RealType, std::size_t w, class Engine>
inline std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::false_)
{
    int bucket = uniform_int_distribution<>(0, (1 << w) - 1)(eng);
    RealType r = uniform_01<RealType>()(eng);
    return std::make_pair(r, bucket);
}

template<class RealType, std::size_t w, class Engine>
inline std::pair<RealType, int> generate_int_float_pair(Engine& eng)
{
    typedef typename Engine::result_type base_result;
    return generate_int_float_pair<RealType, w>(eng,
        boost::is_integral<base_result>());
}

template<class RealType = double>
struct unit_normal_distribution
{
    template<class Engine>
    RealType operator()(Engine& eng) {
        const double * const table_x = normal_table<double>::table_x;
        const double * const table_y = normal_table<double>::table_y;
        for(;;) {
            std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
            int i = vals.second;
            int sign = (i & 1) * 2 - 1;
            i = i >> 1;
            RealType x = vals.first * RealType(table_x[i]);
            if(x < table_x[i + 1]) return x * sign;
            if(i == 0) return generate_tail(eng) * sign;
            RealType y = RealType(table_y[i]) + uniform_01<RealType>()(eng) * RealType(table_y[i + 1] - table_y[i]);
            if (y < f(x)) return x * sign;
        }
    }
    static RealType f(RealType x) {
        using std::exp;
        return exp(-x*x/2);
    }
    template<class Engine>
    RealType generate_tail(Engine& eng) {
        boost::random::exponential_distribution<RealType> exponential;
        const RealType tail_start = RealType(normal_table<double>::table_x[1]);
        for(;;) {
            RealType x = exponential(eng)/tail_start;
            RealType y = exponential(eng);
            if(2*y > x*x) return x + tail_start;
        }
    }
};

}

// deterministic Box-Muller method, uses trigonometric functions

/**
 * Instantiations of class template normal_distribution model a
 * \random_distribution. Such a distribution produces random numbers
 * @c x distributed with probability density function
 * \f$\displaystyle p(x) =
 *   \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
 * \f$,
 * where mean and sigma are the parameters of the distribution.
 */
template<class RealType = double>
class normal_distribution
{
public:
    typedef RealType input_type;
    typedef RealType result_type;

    class param_type {
    public:
        typedef normal_distribution distribution_type;

        /**
         * Constructs a @c param_type with a given mean and
         * standard deviation.
         *
         * Requires: sigma >= 0
         */
        explicit param_type(RealType mean_arg = RealType(0.0),
                            RealType sigma_arg = RealType(1.0))
          : _mean(mean_arg),
            _sigma(sigma_arg)
        {}

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

        /** Returns the standand deviation of the distribution. */
        RealType sigma() const { return _sigma; }

        /** Writes a @c param_type to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        { os << parm._mean << " " << parm._sigma ; return os; }

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

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

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

    private:
        RealType _mean;
        RealType _sigma;
    };

    /**
     * Constructs a @c normal_distribution object. @c mean and @c sigma are
     * the parameters for the distribution.
     *
     * Requires: sigma >= 0
     */
    explicit normal_distribution(const RealType& mean_arg = RealType(0.0),
                                 const RealType& sigma_arg = RealType(1.0))
      : _mean(mean_arg), _sigma(sigma_arg)
    {
        BOOST_ASSERT(_sigma >= RealType(0));
    }

    /**
     * Constructs a @c normal_distribution object from its parameters.
     */
    explicit normal_distribution(const param_type& parm)
      : _mean(parm.mean()), _sigma(parm.sigma())
    {}

    /**  Returns the mean of the distribution. */
    RealType mean() const { return _mean; }
    /** Returns the standard deviation of the distribution. */
    RealType sigma() const { return _sigma; }

    /** 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(_mean, _sigma); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _mean = parm.mean();
        _sigma = parm.sigma();
    }

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

    /**  Returns a normal variate. */
    template<class Engine>
    result_type operator()(Engine& eng)
    {
        detail::unit_normal_distribution<RealType> impl;
        return impl(eng) * _sigma + _mean;
    }

    /** Returns a normal variate with parameters specified by @c param. */
    template<class URNG>
    result_type operator()(URNG& urng, const param_type& parm)
    {
        return normal_distribution(parm)(urng);
    }

    /** Writes a @c normal_distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution, nd)
    {
        os << nd._mean << " " << nd._sigma;
        return os;
    }

    /** Reads a @c normal_distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution, nd)
    {
        is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
        return is;
    }

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

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

private:
    RealType _mean, _sigma;

};

} // namespace random

using random::normal_distribution;

} // namespace boost

#endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP


code:
/* boost random/piecewise_constant_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_CONSTANT_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_HPP_INCLUDED

#include <vector>
#include <numeric>
#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_constant_distribution models a \random_distribution.
 */
template<class RealType = double, class WeightType = double>
class piecewise_constant_distribution {
public:
    typedef std::size_t input_type;
    typedef RealType result_type;

    class param_type {
    public:

        typedef piecewise_constant_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(WeightType(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 the interval weights.
         * If there are less than two boundaries, then this is equivalent to
         * the default constructor and creates a single interval, [0, 1).
         *
         * The values of the interval boundaries must be strictly
         * increasing, and the number of weights must be one less than
         * 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();
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
                _weights.push_back(WeightType(1));
            } else {
                _weights.reserve(_intervals.size() - 1);
                for(std::size_t i = 0; i < _intervals.size() - 1; ++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.  Each
         * weight is determined by calling the function at the
         * midpoint of the corresponding interval.
         *
         * If the initializer_list contains less 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();
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
                _weights.push_back(WeightType(1));
            } else {
                _weights.reserve(_intervals.size() - 1);
                for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
                    RealType midpoint = (_intervals[i] + _intervals[i + 1]) / 2;
                    _weights.push_back(f(midpoint));
                }
            }
        }
#endif
        /**
         * Constructs a @c param_type object from Boost.Range
         * ranges holding the interval boundaries and the weights.  If
         * there are less 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 one less than 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) {
                _intervals.clear();
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
                _weights.push_back(WeightType(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 midpoints 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 + delta/2));
                _intervals.push_back(xmin + k*delta);
            }
            _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
         * over all the intervals of the distribution.
         */
        std::vector<RealType> densities() const
        {
            RealType sum = std::accumulate(_weights.begin(), _weights.end(),
                                             static_cast<RealType>(0));
            std::vector<RealType> result;
            result.reserve(_weights.size());
            for(std::size_t i = 0; i < _weights.size(); ++i) {
                RealType width = _intervals[i + 1] - _intervals[i];
                result.push_back(_weights[i] / (sum * width));
            }
            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<WeightType> 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_constant_distribution;

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

    /**
     * Creates a new @c piecewise_constant_distribution with
     * a single interval, [0, 1).
     */
    piecewise_constant_distribution()
    {
        _intervals.push_back(RealType(0));
        _intervals.push_back(RealType(1));
    }
    /**
     * Constructs a piecewise_constant_distribution from two iterator ranges
     * containing the interval boundaries and the interval weights.
     * If there are less than two boundaries, then this is equivalent to
     * the default constructor and creates a single interval, [0, 1).
     *
     * The values of the interval boundaries must be strictly
     * increasing, and the number of weights must be one less than
     * the number of interval boundaries.  If there are extra
     * weights, they are ignored.
     *
     * For example,
     *
     * @code
     * double intervals[] = { 0.0, 1.0, 4.0 };
     * double weights[] = { 1.0, 1.0 };
     * piecewise_constant_distribution<> dist(
     *     &intervals[0], &intervals[0] + 3, &weights[0]);
     * @endcode
     *
     * The distribution has a 50% chance of producing a
     * value between 0 and 1 and a 50% chance of producing
     * a value between 1 and 4.
     */
    template<class IntervalIter, class WeightIter>
    piecewise_constant_distribution(IntervalIter first_interval,
                                    IntervalIter last_interval,
                                    WeightIter first_weight)
      : _intervals(first_interval, last_interval)
    {
        if(_intervals.size() < 2) {
            _intervals.clear();
            _intervals.push_back(RealType(0));
            _intervals.push_back(RealType(1));
        } else {
            std::vector<WeightType> actual_weights;
            actual_weights.reserve(_intervals.size() - 1);
            for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
                actual_weights.push_back(*first_weight++);
            }
            typedef discrete_distribution<std::size_t, WeightType> bins_type;
            typename bins_type::param_type bins_param(actual_weights);
            _bins.param(bins_param);
        }
    }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    /**
     * Constructs a piecewise_constant_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
     * midpoint of the corresponding interval.
     *
     * If the initializer_list contains less 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_constant_distribution(std::initializer_list<T> il, F f)
      : _intervals(il.begin(), il.end())
    {
        if(_intervals.size() < 2) {
            _intervals.clear();
            _intervals.push_back(RealType(0));
            _intervals.push_back(RealType(1));
        } else {
            std::vector<WeightType> actual_weights;
            actual_weights.reserve(_intervals.size() - 1);
            for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
                RealType midpoint = (_intervals[i] + _intervals[i + 1]) / 2;
                actual_weights.push_back(f(midpoint));
            }
            typedef discrete_distribution<std::size_t, WeightType> bins_type;
            typename bins_type::param_type bins_param(actual_weights);
            _bins.param(bins_param);
        }
    }
#endif
    /**
     * Constructs a piecewise_constant_distribution from Boost.Range
     * ranges holding the interval boundaries and the weights.  If
     * there are less 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 one less than the number of
     * interval boundaries.
     */
    template<class IntervalsRange, class WeightsRange>
    piecewise_constant_distribution(const IntervalsRange& intervals_arg,
                                    const WeightsRange& weights_arg)
      : _bins(weights_arg),
        _intervals(boost::begin(intervals_arg), boost::end(intervals_arg))
    {
        if(_intervals.size() < 2) {
            _intervals.clear();
            _intervals.push_back(RealType(0));
            _intervals.push_back(RealType(1));
        }
    }
    /**
     * Constructs a piecewise_constant_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 midpoints of the
     * intervals.
     */
    template<class F>
    piecewise_constant_distribution(std::size_t nw,
                                    RealType xmin,
                                    RealType xmax,
                                    F f)
      : _bins(nw, xmin, xmax, f)
    {
        if(nw == 0) { nw = 1; }
        RealType delta = (xmax - xmin) / nw;
        _intervals.reserve(nw + 1);
        for(std::size_t i = 0; i < nw; ++i) {
            _intervals.push_back(xmin + i * delta);
        }
        _intervals.push_back(xmax);
    }
    /**
     * Constructs a piecewise_constant_distribution from its parameters.
     */
    explicit piecewise_constant_distribution(const param_type& parm)
      : _bins(parm._weights),
        _intervals(parm._intervals)
    {
    }

    /**
     * Returns a value distributed according to the parameters of the
     * piecewist_constant_distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng) const
    {
        std::size_t i = _bins(urng);
        return uniform_real<RealType>(_intervals[i], _intervals[i+1])(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_constant_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 density
     * over each interval.
     */
    std::vector<RealType> densities() const
    {
        std::vector<RealType> result(_bins.probabilities());
        for(std::size_t i = 0; i < result.size(); ++i) {
            result[i] /= (_intervals[i+1] - _intervals[i]);
        }
        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, _bins.probabilities());
    }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        std::vector<RealType> new_intervals(parm._intervals);
        typedef discrete_distribution<std::size_t, WeightType> bins_type;
        typename bins_type::param_type bins_param(parm._weights);
        _bins.param(bins_param);
        _intervals.swap(new_intervals);
    }

    /**
     * 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_constant_distribution, pcd)
    {
        os << pcd.param();
        return os;
    }

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

private:
    discrete_distribution<std::size_t, WeightType> _bins;
    std::vector<RealType> _intervals;
};

}
}

#endif