mkl.hpp

Go to the documentation of this file.

//============================================================================
// MCKL/include/mckl/random/mkl.hpp
//----------------------------------------------------------------------------
// MCKL: Monte Carlo Kernel Library
//----------------------------------------------------------------------------
// Copyright (c) 2013-2018, Yan Zhou
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
//
//   Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//============================================================================

#ifndef MCKL_RANDOM_MKL_HPP
#define MCKL_RANDOM_MKL_HPP

#include <mckl/random/internal/common.hpp>
#include <mckl/random/uniform_real_distribution.hpp>
#include <mkl_vsl.h>

namespace mckl {

namespace internal {

#if MCKL_NO_RUNTIME_ASSERT

inline int mkl_error_check(int status, const char *, const char *)
{
    return status;
}

#else // MCKL_NO_RUNTIME_ASSERT

inline int mkl_error_check(int status, const char *cpp, const char *c)
{
    if (status == VSL_ERROR_OK)
        return status;

    std::string msg;
    msg += "**";
    msg += cpp;
    msg += "**";
    msg += " failed";
    msg += "; MKL function: ";
    msg += c;
    msg += "; Error code: ";
    msg += std::to_string(status);

    runtime_assert((status == VSL_ERROR_OK), msg.c_str());

    return status;
}

#endif // MCKL_NO_RUNTIME_ASSERT

} // namespace internal

class MKLStream
{
  public:
    explicit MKLStream(::VSLStreamStatePtr ptr = nullptr) : ptr_(nullptr)
    {
        reset(ptr);
    }

    MKLStream(MKL_INT brng, MKL_UINT seed) : ptr_(nullptr)
    {
        reset(brng, seed);
    }

    MKLStream(MKL_INT brng, MKL_INT n, unsigned *params) : ptr_(nullptr)
    {
        reset(brng, n, params);
    }

    MKLStream(const MKLStream &other) : ptr_(nullptr)
    {
        ::VSLStreamStatePtr ptr = nullptr;
        if (internal::mkl_error_check(::vslCopyStream(&ptr, other.ptr_),
                "MKLStream::MKLStream", "::vslCopyStream") == VSL_ERROR_OK) {
            reset(ptr);
        }
    }

    MKLStream &operator=(const MKLStream &other)
    {
        if (this != &other) {
            ::VSLStreamStatePtr ptr = nullptr;
            if (internal::mkl_error_check(::vslCopyStream(&ptr, other.ptr_),
                    "MKLStream::operator=", "::vslCopyStream") ==
                VSL_ERROR_OK) {
                reset(ptr);
            }
        }

        return *this;
    }

    MKLStream(MKLStream &&other) : ptr_(other.ptr_) { other.ptr_ = nullptr; }

    MKLStream &operator=(MKLStream &&other)
    {
        if (this != &other) {
            release();
            ptr_ = other.ptr_;
            other.ptr_ = nullptr;
        }

        return *this;
    }

    ~MKLStream() { release(); }

    int reset(::VSLStreamStatePtr ptr)
    {
        int status = release();
        ptr_ = ptr;

        return status;
    }

    int reset(MKL_INT brng, MKL_UINT seed)
    {
        ::VSLStreamStatePtr ptr = nullptr;
        int status =
            internal::mkl_error_check(::vslNewStream(&ptr, brng, seed),
                "MKLStream::reset", "::vslNewStream");
        if (status == VSL_ERROR_OK) {
            reset(ptr);
        }

        return status;
    }

    int reset(MKL_INT brng, MKL_INT n, unsigned *params)
    {
        ::VSLStreamStatePtr ptr = nullptr;
        int status =
            internal::mkl_error_check(::vslNewStreamEx(&ptr, brng, n, params),
                "MKLStream::reset", "::vslNewStreamEx");
        if (status == VSL_ERROR_OK) {
            reset(ptr);
        }

        return status;
    }

    int release()
    {
        if (ptr_ == nullptr) {
            return VSL_ERROR_OK;
        }

        int status = internal::mkl_error_check(::vslDeleteStream(&ptr_),
            "MKLStream::release", "::vslDeleteStream");
        ptr_ = nullptr;

        return status;
    }

    bool empty() const { return ptr_ == nullptr; }

    int save_f(const std::string &fname) const
    {
        return internal::mkl_error_check(::vslSaveStreamF(ptr_, fname.c_str()),
            "MKLStream::save_f", "::vslSaveStreamF");
    }

    int load_f(const std::string &fname)
    {
        ::VSLStreamStatePtr ptr = nullptr;
        int status =
            internal::mkl_error_check(::vslSaveStreamF(&ptr, fname.c_str()),
                "MKLStream::load_f", "::vslSaveStreamF");
        if (status == VSL_ERROR_OK) {
            reset(ptr);
        }

        return status;
    }

    int save_m(char *memptr) const
    {
        return internal::mkl_error_check(::vslSaveStreamM(ptr_, memptr),
            "MKLStream::save_m", "::vslSaveStreamM");
    }

    int load_m(const char *memptr)
    {
        ::VSLStreamStatePtr ptr = nullptr;
        int status = internal::mkl_error_check(::vslLoadStreamM(&ptr, memptr),
            "MKLStream::load_m", "::vslLoadStreamM");
        if (status == VSL_ERROR_OK) {
            reset(ptr);
        }

        return status;
    }

    int get_size() const { return ::vslGetStreamSize(ptr_); }

    int leapfrog(MKL_INT k, MKL_INT nstreams)
    {
        return internal::mkl_error_check(
            ::vslLeapfrogStream(ptr_, k, nstreams), "MKLStream::leapfrog",
            "::vslLeapfrogStream");
    }

    int skip_ahead(long long nskip)
    {
        return internal::mkl_error_check(::vslSkipAheadStream(ptr_, nskip),
            "MKLStream::skip_ahead", "::vslSkipAheadStream");
    }

    int get_brng() const { return ::vslGetStreamStateBrng(ptr_); }

    static int get_num_reg_brngs() { return ::vslGetNumRegBrngs(); }

    static int get_brng_properties(
        MKL_INT brng, ::VSLBRngProperties *properties)
    {
        return internal::mkl_error_check(
            ::vslGetBrngProperties(brng, properties),
            "MKLStream::get_brng_properties", "::vslGetBrngProperties");
    }

    static bool has_leap_frog(MKL_INT brng)
    {
        MKLStream stream(brng, 1);

        return ::vslLeapfrogStream(stream.ptr_, 1, 2) == VSL_ERROR_OK;
    }

    static bool has_skip_ahead(MKL_INT brng)
    {
        MKLStream stream(brng, 1);

        return ::vslSkipAheadStream(stream.ptr_, 1) == VSL_ERROR_OK;
    }

    static bool has_uniform_bits32(
        MKL_INT brng, MKL_INT method = VSL_RNG_METHOD_UNIFORMBITS32_STD)
    {
        MKLStream stream(brng, 1);
        unsigned r;

        return ::viRngUniformBits32(method, stream.ptr_, 1, &r) ==
            VSL_ERROR_OK;
    }

    static bool has_uniform_bits64(
        MKL_INT brng, MKL_INT method = VSL_RNG_METHOD_UNIFORMBITS64_STD)
    {
        MKLStream stream(brng, 1);
        unsigned MKL_INT64 r;

        return ::viRngUniformBits64(method, stream.ptr_, 1, &r) ==
            VSL_ERROR_OK;
    }

    int uniform(MKL_INT n, float *r, float a, float b,
        MKL_INT method = VSL_RNG_METHOD_UNIFORM_STD)
    {
        return internal::mkl_error_check(
            ::vsRngUniform(method, ptr_, n, r, a, b), "MKLStream::uniform",
            "::vsRngUniform");
    }

    int uniform(MKL_INT n, double *r, double a, double b,
        MKL_INT method = VSL_RNG_METHOD_UNIFORM_STD)
    {
        return internal::mkl_error_check(
            ::vdRngUniform(method, ptr_, n, r, a, b), "MKLStream::uniform",
            "::vdRngUniform");
    }

    int gaussian(MKL_INT n, float *r, float a, float sigma,
        MKL_INT method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vsRngGaussian(method, ptr_, n, r, a, sigma),
            "MKLStream::gaussian", "::vsRngGaussian");
    }

    int gaussian(MKL_INT n, double *r, double a, double sigma,
        MKL_INT method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vdRngGaussian(method, ptr_, n, r, a, sigma),
            "MKLStream::gaussian", "::vdRngGaussian");
    }

    int gaussian_mv(MKL_INT n, float *r, MKL_INT dimen, MKL_INT mstorage,
        const float *a, const float *t,
        MKL_INT method = VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vsRngGaussianMV(method, ptr_, n, r, dimen, mstorage, a, t),
            "MKLStream::gaussian_mv", "::vsRngGaussianMV");
    }

    int gaussian_mv(MKL_INT n, double *r, MKL_INT dimen, MKL_INT mstorage,
        const double *a, const double *t,
        MKL_INT method = VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vdRngGaussianMV(method, ptr_, n, r, dimen, mstorage, a, t),
            "MKLStream::gaussian_mv", "::vdRngGaussianMV");
    }

    int exponential(MKL_INT n, float *r, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_EXPONENTIAL_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngExponential(method, ptr_, n, r, a, beta),
            "MKLStream::exponential", "::vsRngExponential");
    }

    int exponential(MKL_INT n, double *r, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_EXPONENTIAL_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngExponential(method, ptr_, n, r, a, beta),
            "MKLStream::exponential", "::vdRngExponential");
    }

    int laplace(MKL_INT n, float *r, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_LAPLACE_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngLaplace(method, ptr_, n, r, a, beta), "MKLStream::laplace",
            "::vsRngLaplace");
    }

    int laplace(MKL_INT n, double *r, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_LAPLACE_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngLaplace(method, ptr_, n, r, a, beta), "MKLStream::laplace",
            "::vdRngLaplace");
    }

    int weibull(MKL_INT n, float *r, float alpha, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_WEIBULL_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngWeibull(method, ptr_, n, r, alpha, a, beta),
            "MKLStream::weibull", "::vsRngWeibull");
    }

    int weibull(MKL_INT n, double *r, double alpha, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_WEIBULL_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngWeibull(method, ptr_, n, r, alpha, a, beta),
            "MKLStream::weibull", "::vdRngWeibull");
    }

    int cauchy(MKL_INT n, float *r, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_CAUCHY_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngCauchy(method, ptr_, n, r, a, beta), "MKLStream::cauchy",
            "::vsRngCauchy");
    }

    int cauchy(MKL_INT n, double *r, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_CAUCHY_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngCauchy(method, ptr_, n, r, a, beta), "MKLStream::cauchy",
            "::vdRngCauchy");
    }

    int rayleigh(MKL_INT n, float *r, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_RAYLEIGH_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngRayleigh(method, ptr_, n, r, a, beta),
            "MKLStream::rayleigh", "::vsRngRayleigh");
    }

    int rayleigh(MKL_INT n, double *r, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_RAYLEIGH_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngRayleigh(method, ptr_, n, r, a, beta),
            "MKLStream::rayleigh", "::vdRngRayleigh");
    }

    int lognormal(MKL_INT n, float *r, float a, float sigma, float b,
        float beta, MKL_INT method = VSL_RNG_METHOD_LOGNORMAL_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vsRngLognormal(method, ptr_, n, r, a, sigma, b, beta),
            "MKLStream::lognormal", "::vsRngLognormal");
    }

    int lognormal(MKL_INT n, double *r, double a, double sigma, double b,
        double beta, MKL_INT method = VSL_RNG_METHOD_LOGNORMAL_BOXMULLER2)
    {
        return internal::mkl_error_check(
            ::vdRngLognormal(method, ptr_, n, r, a, sigma, b, beta),
            "MKLStream::lognormal", "::vdRngLognormal");
    }

    int gumbel(MKL_INT n, float *r, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_GUMBEL_ICDF)
    {
        return internal::mkl_error_check(
            ::vsRngGumbel(method, ptr_, n, r, a, beta), "MKLStream::gumbel",
            "::vsRngGumbel");
    }

    int gumbel(MKL_INT n, double *r, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_GUMBEL_ICDF)
    {
        return internal::mkl_error_check(
            ::vdRngGumbel(method, ptr_, n, r, a, beta), "MKLStream::gumbel",
            "::vdRngGumbel");
    }

    int gamma(MKL_INT n, float *r, float alpha, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_GAMMA_GNORM)
    {
        return internal::mkl_error_check(
            ::vsRngGamma(method, ptr_, n, r, alpha, a, beta),
            "MKLStream::gamma", "::vsRngGamma");
    }

    int gamma(MKL_INT n, double *r, double alpha, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_GAMMA_GNORM)
    {
        return internal::mkl_error_check(
            ::vdRngGamma(method, ptr_, n, r, alpha, a, beta),
            "MKLStream::gamma", "::vdRngGamma");
    }

    int beta(MKL_INT n, float *r, float p, float q, float a, float beta,
        MKL_INT method = VSL_RNG_METHOD_BETA_CJA)
    {
        return internal::mkl_error_check(
            ::vsRngBeta(method, ptr_, n, r, p, q, a, beta), "MKLStream::beta",
            "::vsRngBeta");
    }

    int beta(MKL_INT n, double *r, double p, double q, double a, double beta,
        MKL_INT method = VSL_RNG_METHOD_BETA_CJA)
    {
        return internal::mkl_error_check(
            ::vdRngBeta(method, ptr_, n, r, p, q, a, beta), "MKLStream::beta",
            "::vdRngBeta");
    }

    int uniform(MKL_INT n, int *r, int a, int b,
        MKL_INT method = VSL_RNG_METHOD_UNIFORM_STD)
    {
        return internal::mkl_error_check(
            ::viRngUniform(method, ptr_, n, r, a, b), "MKLStream::uniform",
            "::viRngUniform");
    }

    int uniform_bits(MKL_INT n, unsigned *r,
        MKL_INT method = VSL_RNG_METHOD_UNIFORMBITS_STD)
    {
        return internal::mkl_error_check(
            ::viRngUniformBits(method, ptr_, n, r), "MKLStream::uniform_bits",
            "::viRngUniformBits");
    }

    int uniform_bits32(MKL_INT n, unsigned *r,
        MKL_INT method = VSL_RNG_METHOD_UNIFORMBITS32_STD)
    {
        return internal::mkl_error_check(
            ::viRngUniformBits32(method, ptr_, n, r),
            "MKLStream::uniform_bits32", "::viRngUniformBits32");
    }

    int uniform_bits64(MKL_INT n, unsigned MKL_INT64 *r,
        MKL_INT method = VSL_RNG_METHOD_UNIFORMBITS64_STD)
    {
        return internal::mkl_error_check(
            ::viRngUniformBits64(method, ptr_, n, r),
            "MKLStream::uniform_bits64", "::viRngUniformBits64");
    }

    int bernoulli(MKL_INT n, int *r, double p,
        MKL_INT method = VSL_RNG_METHOD_BERNOULLI_ICDF)
    {
        return internal::mkl_error_check(
            ::viRngBernoulli(method, ptr_, n, r, p), "MKLStream::bernoulli",
            "::viRngBernoulli");
    }

    int geometric(MKL_INT n, int *r, double p,
        MKL_INT method = VSL_RNG_METHOD_GEOMETRIC_ICDF)
    {
        return internal::mkl_error_check(
            ::viRngGeometric(method, ptr_, n, r, p), "MKLStream::geometric",
            "::viRngGeometric");
    }

    int binomial(MKL_INT n, int *r, int ntrial, double p,
        MKL_INT method = VSL_RNG_METHOD_BINOMIAL_BTPE)
    {
        return internal::mkl_error_check(
            ::viRngBinomial(method, ptr_, n, r, ntrial, p),
            "MKLStream::binomial", "::viRngBinomial");
    }

    int hypergeometric(MKL_INT n, int *r, int l, int s, int m,
        MKL_INT method = VSL_RNG_METHOD_HYPERGEOMETRIC_H2PE)
    {
        return internal::mkl_error_check(
            ::viRngHypergeometric(method, ptr_, n, r, l, s, m),
            "MKLStream::hypergeometric", "::viRngHypergeometric");
    }

    int poisson(MKL_INT n, int *r, double lambda,
        MKL_INT method = VSL_RNG_METHOD_POISSON_POISNORM)
    {
        return internal::mkl_error_check(
            ::viRngPoisson(method, ptr_, n, r, lambda), "MKLStream::poisson",
            "::viRngPoisson");
    }

    int poisson_v(MKL_INT n, int *r, const double *lambda,
        MKL_INT method = VSL_RNG_METHOD_POISSONV_POISNORM)
    {
        return internal::mkl_error_check(
            ::viRngPoissonV(method, ptr_, n, r, lambda),
            "MKLStream::poisson_v", "::viRngPoissonV");
    }

    int neg_binomial(MKL_INT n, int *r, double a, double p,
        MKL_INT method = VSL_RNG_METHOD_NEGBINOMIAL_NBAR)
    {
        return internal::mkl_error_check(
            ::viRngNegbinomial(method, ptr_, n, r, a, p),
            "MKLStream::neg_binomial", "::viRngNegbinomial");
    }

  private:
    ::VSLStreamStatePtr ptr_;
}; // class MKLStream

inline bool operator==(const MKLStream &stream1, const MKLStream &stream2)
{
    int brng1 = stream1.get_brng();
    int brng2 = stream2.get_brng();
    if (brng1 != brng2) {
        return false;
    }
    if (brng1 == VSL_BRNG_NONDETERM) {
        return false;
    }

    std::size_t n = static_cast<std::size_t>(stream1.get_size());
    Vector<char> s1(n);
    Vector<char> s2(n);
    stream1.save_m(s1.data());
    stream2.save_m(s2.data());
    if (s1 != s2) {
        return false;
    }

    return true;
}

inline bool operator!=(const MKLStream &stream1, const MKLStream &stream2)
{
    return !(stream1 == stream2);
}

template <typename CharT, typename Traits>
inline std::basic_ostream<CharT, Traits> &operator<<(
    std::basic_ostream<CharT, Traits> &os, const MKLStream &stream)
{
    if (!os) {
        return os;
    }

    std::size_t n = static_cast<std::size_t>(stream.get_size());
    std::size_t m = sizeof(std::uintmax_t);
    if (n % m != 0) {
        n += m - n % m;
    }
    n /= m;
    Vector<std::uintmax_t> s(n);
    stream.save_m(reinterpret_cast<char *>(s.data()));

    os << stream.get_brng() << ' ';
    os << s;

    return os;
}

template <typename CharT, typename Traits>
inline std::basic_istream<CharT, Traits> &operator>>(
    std::basic_istream<CharT, Traits> &is, MKLStream &stream)
{
    if (!is) {
        return is;
    }

    MKL_INT brng;
    Vector<std::uintmax_t> s;
    is >> std::ws >> brng;
    is >> std::ws >> s;

    if (is) {
        MKLStream tmp;
        tmp.load_m(reinterpret_cast<const char *>(s.data()));
        stream = std::move(tmp);
    }

    return is;
}

namespace internal {

template <MKL_INT BRNG>
class MKLMaxOffset : public std::integral_constant<MKL_INT, 0>
{
}; // class MKLMaxOffset;

template <>
class MKLMaxOffset<VSL_BRNG_MT2203>
    : public std::integral_constant<MKL_INT, 6024>
{
}; // class MKLMaxOffset

template <>
class MKLMaxOffset<VSL_BRNG_WH> : public std::integral_constant<MKL_INT, 273>
{
}; // class MKLMaxOffset

template <MKL_INT BRNG, MKL_INT MaxOffset = MKLMaxOffset<BRNG>::value>
class MKLOffset
{
  public:
    static MKL_INT eval(MKL_INT offset) { return BRNG + offset % MaxOffset; }
}; // class MKLOffset

template <MKL_INT BRNG>
class MKLOffset<BRNG, 0>
{
  public:
    static MKL_INT eval(MKL_INT) { return BRNG; }
}; // class MKLOffset

template <MKL_INT, int>
class MKLUniformBits;

template <MKL_INT BRNG>
class MKLUniformBits<BRNG, 32>
{
  public:
    using result_type = unsigned;

    static constexpr result_type min()
    {
        return std::numeric_limits<result_type>::min();
    }

    static constexpr result_type max()
    {
        return std::numeric_limits<result_type>::max();
    }

    static void eval(MKLStream &stream, MKL_INT n, result_type *r)
    {
        stream.uniform_bits32(n, r);
    }
}; // class MKLUniformBits

template <MKL_INT BRNG>
class MKLUniformBits<BRNG, 64>
{
  public:
    using result_type = unsigned MKL_INT64;

    static constexpr result_type min()
    {
        return std::numeric_limits<result_type>::min();
    }

    static constexpr result_type max()
    {
        return std::numeric_limits<result_type>::max();
    }

    static void eval(MKLStream &stream, MKL_INT n, result_type *r)
    {
        stream.uniform_bits64(n, r);
    }
}; // class MKLUniformBits

} // namespace internal

template <MKL_INT BRNG, int Bits>
class MKLEngine
{
    static_assert(Bits == 32 || Bits == 64,
        "**MKLEngine** used with bits other than 32, or 64");

  public:
    using result_type =
        typename internal::MKLUniformBits<BRNG, Bits>::result_type;

  private:
    template <typename T>
    using is_seed_seq = internal::is_seed_seq<T, MKLEngine<BRNG, Bits>>;

  public:
    explicit MKLEngine(result_type s = 1) : index_(M_) { seed(s); }

    template <typename SeedSeq>
    explicit MKLEngine(SeedSeq &seq,
        std::enable_if_t<is_seed_seq<SeedSeq>::value> * = nullptr)
        : index_(M_)
    {
        seed(seq);
    }

    MKLEngine(MKL_INT offset, result_type s) : index_(M_)
    {
        static_assert(internal::MKLMaxOffset<BRNG>::value > 0,
            "**MKLEngine** does not support offseting");

        seed(offset, s);
    }

    template <typename SeedSeq>
    explicit MKLEngine(MKL_INT offset, SeedSeq &seq,
        std::enable_if_t<is_seed_seq<SeedSeq>::value> * = nullptr)
        : index_(M_)
    {
        static_assert(internal::MKLMaxOffset<BRNG>::value > 0,
            "**MKLEngine** does not support offseting");

        seed(offset, seq);
    }

    void seed(result_type s)
    {
        s %= static_cast<result_type>(std::numeric_limits<MKL_UINT>::max());
        MKL_INT brng = stream_.empty() ? BRNG : stream_.get_brng();
        stream_.reset(brng, static_cast<MKL_UINT>(s));
        index_ = M_;
    }

    template <typename SeedSeq>
    void seed(SeedSeq &seq,
        std::enable_if_t<is_seed_seq<SeedSeq>::value> * = nullptr)
    {
        MKL_INT brng = stream_.empty() ? BRNG : stream_.get_brng();
        Vector<MKL_UINT> params;
        MKL_INT n = seed_params(brng, seq, params);
        stream_.reset(brng, n, params.data());
        index_ = M_;
    }

    void seed(MKL_INT offset, result_type s)
    {
        static_assert(internal::MKLMaxOffset<BRNG>::value > 0,
            "**MKLEngine** does not support offseting");

        s %= static_cast<result_type>(std::numeric_limits<MKL_UINT>::max());
        MKL_INT brng = internal::MKLOffset<BRNG>::eval(offset);
        stream_.reset(brng, static_cast<MKL_UINT>(s));
        index_ = M_;
    }

    template <typename SeedSeq>
    void seed(MKL_INT offset, SeedSeq &seq,
        std::enable_if_t<is_seed_seq<SeedSeq>::value> * = nullptr)
    {
        static_assert(internal::MKLMaxOffset<BRNG>::value > 0,
            "**MKLEngine** does not support offseting");

        MKL_INT brng = internal::MKLOffset<BRNG>::eval(offset);
        Vector<unsigned> params;
        MKL_INT n = seed_params(brng, seq, params);
        stream_.reset(brng, n, params.data());
        index_ = M_;
    }

    result_type operator()()
    {
        if (index_ == M_) {
            generate();
            index_ = 0;
        }

        return result_[index_++];
    }

    void operator()(std::size_t n, result_type *r)
    {
        internal::size_check<MKL_INT>(n, "MKLEngine::operator()");

        const std::size_t remain = M_ - index_;

        if (n <= remain) {
            std::memcpy(r, result_.data() + index_, sizeof(result_type) * n);
            index_ += n;
            return;
        }

        std::memcpy(r, result_.data() + index_, sizeof(result_type) * remain);
        r += remain;
        n -= remain;
        index_ = M_;

        generate(n, r);
    }

    std::size_t discard()
    {
        const std::size_t remain = M_ - index_;
        index_ = M_;

        return remain;
    }

    void discard(long long nskip)
    {

        if (nskip <= 0) {
            return;
        }

        const long long remain = static_cast<long long>(M_ - index_);
        if (nskip <= remain) {
            index_ += static_cast<std::size_t>(nskip);
            return;
        }
        nskip -= remain;
        index_ = M_;

        ::VSLBRngProperties properties;
        MKLStream::get_brng_properties(stream_.get_brng(), &properties);
        const int bits = properties.NBits;
        const long long M = static_cast<long long>(M_);
        long long m = nskip / M * M;
        if (Bits >= bits) {
            m *= Bits / bits + (Bits % bits == 0 ? 0 : 1);
        }
        switch (stream_.get_brng()) {
            case VSL_BRNG_ARS5:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_PHILOX4X32X10:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_MCG31:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_MRG32K3A:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_MCG59:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_WH:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_MT19937:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_SFMT19937:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_SOBOL:
                stream_.skip_ahead(m);
                break;
            case VSL_BRNG_NIEDERR:
                stream_.skip_ahead(m);
                break;
            default:
                while (nskip >= M) {
                    generate();
                    nskip -= M;
                }
        };
        generate();
        index_ = static_cast<std::size_t>(nskip % M);
    }

    static constexpr result_type min()
    {
        return internal::MKLUniformBits<BRNG, Bits>::min();
    }

    static constexpr result_type max()
    {
        return internal::MKLUniformBits<BRNG, Bits>::max();
    }

    MKLStream &stream() { return stream_; }
    const MKLStream &stream() const { return stream_; }

    friend bool operator==(
        const MKLEngine<BRNG, Bits> &eng1, const MKLEngine<BRNG, Bits> &eng2)
    {
        if (eng1.stream_ != eng2.stream_) {
            return false;
        }
        if (eng1.result_ != eng2.result_) {
            return false;
        }
        if (eng1.index_ != eng2.index_) {
            return false;
        }
        return true;
    }

    friend bool operator!=(
        const MKLEngine<BRNG, Bits> &eng1, const MKLEngine<BRNG, Bits> &eng2)
    {
        return !(eng1 == eng2);
    }

    template <typename CharT, typename Traits>
    friend std::basic_ostream<CharT, Traits> &operator<<(
        std::basic_ostream<CharT, Traits> &os,
        const MKLEngine<BRNG, Bits> &eng)
    {
        if (!os) {
            return os;
        }

        os << eng.stream_ << ' ';
        os << eng.result_ << ' ';
        os << eng.index_;

        return os;
    }

    template <typename CharT, typename Traits>
    friend std::basic_istream<CharT, Traits> &operator>>(
        std::basic_istream<CharT, Traits> &is, MKLEngine<BRNG, Bits> &eng)
    {
        if (!is) {
            return is;
        }

        MKLStream stream;
        std::array<result_type, M_> result;
        unsigned index;
        is >> std::ws >> stream;
        is >> std::ws >> result;
        is >> std::ws >> index;

        if (is) {
            eng.stream_ = std::move(stream);
            eng.result_ = result;
            eng.index_ = index;
        }

        return is;
    }

  private:
    static constexpr std::size_t M_ = 256;

    std::array<result_type, M_> result_;
    MKLStream stream_;
    std::size_t index_;

    void generate()
    {
        internal::MKLUniformBits<BRNG, Bits>::eval(
            stream_, static_cast<MKL_INT>(M_), result_.data());
    }

    void generate(std::size_t n, result_type *r)
    {
        internal::MKLUniformBits<BRNG, Bits>::eval(
            stream_, static_cast<MKL_INT>(n), r);
    }

    template <typename SeedSeq>
    MKL_INT seed_params(MKL_INT brng, SeedSeq &seq, Vector<unsigned> &params)
    {
        ::VSLBRngProperties properties;
        MKLStream::get_brng_properties(brng, &properties);
        MKL_INT n = properties.NSeeds;
        params.resize(static_cast<std::size_t>(n));
        seq.generate(params.begin(), params.end());

        return n;
    }
}; // class MKLEngine

template <MKL_INT BRNG, int Bits>
inline void rand(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    typename MKLEngine<BRNG, Bits>::result_type *r)
{
    rng(n, r);
}

using MKL_ARS5 = MKLEngine<VSL_BRNG_ARS5, 32>;

using MKL_ARS5_64 = MKLEngine<VSL_BRNG_ARS5, 64>;

using MKL_PHILOX4X32X10 = MKLEngine<VSL_BRNG_PHILOX4X32X10, 32>;

using MKL_PHILOX4X32X10_64 = MKLEngine<VSL_BRNG_PHILOX4X32X10, 64>;

using MKL_MCG59 = MKLEngine<VSL_BRNG_MCG59, 32>;

using MKL_MCG59_64 = MKLEngine<VSL_BRNG_MCG59, 64>;

using MKL_MT19937 = MKLEngine<VSL_BRNG_MT19937, 32>;

using MKL_MT19937_64 = MKLEngine<VSL_BRNG_MT19937, 64>;

using MKL_MT2203 = MKLEngine<VSL_BRNG_MT2203, 32>;

using MKL_MT2203_64 = MKLEngine<VSL_BRNG_MT2203, 64>;

using MKL_SFMT19937 = MKLEngine<VSL_BRNG_SFMT19937, 32>;

using MKL_SFMT19937_64 = MKLEngine<VSL_BRNG_SFMT19937, 64>;

using MKL_NONDETERM = MKLEngine<VSL_BRNG_NONDETERM, 32>;

using MKL_NONDETERM_64 = MKLEngine<VSL_BRNG_NONDETERM, 64>;

namespace internal {

template <typename RNGType>
class MKLStreamState
{
  public:
    unsigned reserved1[2];
    unsigned reserved2[2];
    RNGType rng;
}; // class MKLStreamState

template <typename RNGType>
inline constexpr int mkl_nseeds()
{
    return sizeof(SeedType<RNGType>) < sizeof(unsigned) ?
        1 :
        sizeof(SeedType<RNGType>) / sizeof(unsigned);
}

template <typename RNGType>
inline int mkl_init(RNGType &rng, int n, const unsigned *param, std::true_type)
{
    if (n == 0) {
        new (static_cast<void *>(&rng)) RNGType();
    } else {
        new (static_cast<void *>(&rng))
            RNGType(static_cast<typename RNGType::result_type>(param[0]));
    }

    return 0;
}

template <typename RNGType>
inline int mkl_init(
    int method, ::VSLStreamStatePtr stream, int n, const unsigned *param)
{
    RNGType &rng = (*static_cast<MKLStreamState<RNGType> *>(stream)).rng;

    if (method == VSL_INIT_METHOD_STANDARD) {
        if (n == 0) {
            new (static_cast<void *>(&rng)) RNGType();
        } else {
            constexpr std::size_t ns = mkl_nseeds<RNGType>();
            union {
                SeedType<RNGType> seed;
                std::array<unsigned, ns> useed;
            } buf;
            std::copy_n(param, std::min(static_cast<std::size_t>(n), ns),
                buf.useed.data());
            new (static_cast<void *>(&rng)) RNGType(buf.seed);
        }
    }

    if (method == VSL_INIT_METHOD_LEAPFROG) {
        return VSL_RNG_ERROR_LEAPFROG_UNSUPPORTED;
    }

    if (method == VSL_INIT_METHOD_SKIPAHEAD) {
        rng.discard(static_cast<unsigned>(n));
    }

    return 0;
}

template <typename RNGType, typename RealType>
inline int mkl_uniform_real(
    ::VSLStreamStatePtr stream, int n, RealType *r, RealType a, RealType b)
{
    RNGType &rng = (*static_cast<MKLStreamState<RNGType> *>(stream)).rng;
    uniform_real_distribution(rng, static_cast<std::size_t>(n), r, a, b);

    return 0;
}

template <typename RNGType>
inline int mkl_uniform_int(::VSLStreamStatePtr stream, int n, unsigned *r)
{
    RNGType &rng = (*static_cast<MKLStreamState<RNGType> *>(stream)).rng;
    uniform_bits_distribution(rng, static_cast<std::size_t>(n), r);

    return 0;
}

} // namespace internal

template <typename RNGType>
inline int mkl_brng()
{
    static ::VSLBRngProperties properties = {
        sizeof(internal::MKLStreamState<RNGType>),
        internal::mkl_nseeds<RNGType>(), 1, 4, 32, internal::mkl_init<RNGType>,
        internal::mkl_uniform_real<RNGType, float>,
        internal::mkl_uniform_real<RNGType, double>,
        internal::mkl_uniform_int<RNGType>};
    static int brng = ::vslRegisterBrng(&properties);

    return brng;
}

#if MCKL_USE_MKL_VSL

template <MKL_INT BRNG, int Bits>
inline void beta_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    float *r, float alpha, float beta)
{
    internal::size_check<MKL_INT>(n, "beta_distribution");
    rng.stream().beta(static_cast<MKL_INT>(n), r, alpha, beta, 0, 1);
}

template <MKL_INT BRNG, int Bits>
inline void beta_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    double *r, double alpha, double beta)
{
    internal::size_check<MKL_INT>(n, "beta_distribution");
    rng.stream().beta(static_cast<MKL_INT>(n), r, alpha, beta, 0, 1);
}

template <MKL_INT BRNG, int Bits>
inline void cauchy_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float a, float b)
{
    internal::size_check<MKL_INT>(n, "cauchy_distribution");
    rng.stream().cauchy(static_cast<MKL_INT>(n), r, a, b);
}

template <MKL_INT BRNG, int Bits>
inline void cauchy_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double a, double b)
{
    internal::size_check<MKL_INT>(n, "cauchy_distribution");
    rng.stream().cauchy(static_cast<MKL_INT>(n), r, a, b);
}

template <MKL_INT BRNG, int Bits>
inline void exponential_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float lambda)
{
    internal::size_check<MKL_INT>(n, "exponential_distribution");
    rng.stream().exponential(static_cast<MKL_INT>(n), r, 0, 1 / lambda);
}

template <MKL_INT BRNG, int Bits>
inline void exponential_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double lambda)
{
    internal::size_check<MKL_INT>(n, "exponential_distribution");
    rng.stream().exponential(static_cast<MKL_INT>(n), r, 0, 1 / lambda);
}

template <MKL_INT BRNG, int Bits>
inline void extreme_value_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float a, float b)
{
    internal::size_check<MKL_INT>(n, "extreme_value_distribution");
    rng.stream().gumbel(static_cast<MKL_INT>(n), r, a, b);
    sub(n, 2 * a, r, r);
}

template <MKL_INT BRNG, int Bits>
inline void extreme_value_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double a, double b)
{
    internal::size_check<MKL_INT>(n, "extreme_value_distribution");
    rng.stream().gumbel(static_cast<MKL_INT>(n), r, a, b);
    sub(n, 2 * a, r, r);
}

template <MKL_INT BRNG, int Bits>
inline void gamma_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    float *r, float alpha, float beta)
{
    internal::size_check<MKL_INT>(n, "gamma_distribution");
    rng.stream().gamma(static_cast<MKL_INT>(n), r, alpha, 0, beta);
}

template <MKL_INT BRNG, int Bits>
inline void gamma_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    double *r, double alpha, double beta)
{
    internal::size_check<MKL_INT>(n, "gamma_distribution");
    rng.stream().gamma(static_cast<MKL_INT>(n), r, alpha, 0, beta);
}

template <MKL_INT BRNG, int Bits>
inline void laplace_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    float *r, float location, float scale)
{
    internal::size_check<MKL_INT>(n, "lapace_distribution");
    rng.stream().laplace(static_cast<MKL_INT>(n), r, location, scale);
}

template <MKL_INT BRNG, int Bits>
inline void laplace_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    double *r, double location, double scale)
{
    internal::size_check<MKL_INT>(n, "lapace_distribution");
    rng.stream().laplace(static_cast<MKL_INT>(n), r, location, scale);
}

template <MKL_INT BRNG, int Bits>
inline void lognormal_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float m, float s)
{
    internal::size_check<MKL_INT>(n, "lognormal_distribution");
    rng.stream().lognormal(static_cast<MKL_INT>(n), r, m, s, 0, 1);
}

template <MKL_INT BRNG, int Bits>
inline void lognormal_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double m, double s)
{
    internal::size_check<MKL_INT>(n, "lognormal_distribution");
    rng.stream().lognormal(static_cast<MKL_INT>(n), r, m, s, 0, 1);
}

template <MKL_INT BRNG, int Bits>
inline void normal_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    float *r, float mean, float stddev)
{
    internal::size_check<MKL_INT>(n, "normal_distribution");
    rng.stream().gaussian(static_cast<MKL_INT>(n), r, mean, stddev);
}

template <MKL_INT BRNG, int Bits>
inline void normal_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    double *r, double mean, double stddev)
{
    internal::size_check<MKL_INT>(n, "normal_distribution");
    rng.stream().gaussian(static_cast<MKL_INT>(n), r, mean, stddev);
}

template <MKL_INT BRNG, int Bits>
inline void normal_mv_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    float *r, std::size_t m, const float *mean, const float *chol)
{
    internal::size_check<MKL_INT>(n, "normal_mv_distribution");
    internal::size_check<MKL_INT>(m, "normal_mv_distribution");
    rng.stream().gaussian_mv(static_cast<MKL_INT>(n), r,
        static_cast<MKL_INT>(m), VSL_MATRIX_STORAGE_PACKED, mean, chol);
}

template <MKL_INT BRNG, int Bits>
inline void normal_mv_distribution(MKLEngine<BRNG, Bits> &rng, std::size_t n,
    double *r, std::size_t m, const double *mean, const double *chol)
{
    internal::size_check<MKL_INT>(n, "normal_mv_distribution");
    internal::size_check<MKL_INT>(m, "normal_mv_distribution");
    rng.stream().gaussian_mv(static_cast<MKL_INT>(n), r,
        static_cast<MKL_INT>(m), VSL_MATRIX_STORAGE_PACKED, mean, chol);
}

template <MKL_INT BRNG, int Bits>
inline void rayleigh_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float sigma)
{
    internal::size_check<MKL_INT>(n, "rayleigh_distribution");
    rng.stream().rayleigh(
        static_cast<MKL_INT>(n), r, 0, const_sqrt_2<float>() * sigma);
}

template <MKL_INT BRNG, int Bits>
inline void rayleigh_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double sigma)
{
    internal::size_check<MKL_INT>(n, "rayleigh_distribution");
    rng.stream().rayleigh(
        static_cast<MKL_INT>(n), r, 0, const_sqrt_2<double>() * sigma);
}

template <MKL_INT BRNG, int Bits>
inline void uniform_real_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float a, float b)
{
    internal::size_check<MKL_INT>(n, "uniform_real_distribution");
    rng.stream().uniform(static_cast<MKL_INT>(n), r, a, b);
}

template <MKL_INT BRNG, int Bits>
inline void uniform_real_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double a, double b)
{
    internal::size_check<MKL_INT>(n, "uniform_real_distribution");
    rng.stream().uniform(static_cast<MKL_INT>(n), r, a, b);
}

template <MKL_INT BRNG, int Bits>
inline void weibull_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, float *r, float a, float b)
{
    internal::size_check<MKL_INT>(n, "weibull_distribution");
    rng.stream().weibull(static_cast<MKL_INT>(n), r, a, 0, b);
}

template <MKL_INT BRNG, int Bits>
inline void weibull_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, double *r, double a, double b)
{
    internal::size_check<MKL_INT>(n, "weibull_distribution");
    rng.stream().weibull(static_cast<MKL_INT>(n), r, a, 0, b);
}

template <MKL_INT BRNG, int Bits>
inline void bernoulli_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, int *r, double p)
{
    internal::size_check<MKL_INT>(n, "bernoulli_distribution");
    rng.stream().bernoulli(static_cast<MKL_INT>(n), r, p);
}

template <MKL_INT BRNG, int Bits>
inline void geometric_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, int *r, double p)
{
    internal::size_check<MKL_INT>(n, "geometric_distribution");
    rng.stream().geometric(static_cast<MKL_INT>(n), r, p);
}

template <MKL_INT BRNG, int Bits>
inline void uniform_int_distribution(
    MKLEngine<BRNG, Bits> &rng, std::size_t n, int *r, int a, int b)
{
    internal::size_check<MKL_INT>(n, "uniform_int_distribution");
    if (b < std::numeric_limits<int>::max()) {
        rng.stream().uniform(static_cast<MKL_INT>(n), r, a, b + 1);
    } else if (a > 0) {
        rng.stream().uniform(static_cast<MKL_INT>(n), r, a - 1, b);
        add(n, 1, r, r);
    } else {
        rng.stream().uniform_bits32(
            static_cast<MKL_INT>(n), reinterpret_cast<unsigned *>(r));
    }
}

#endif // MCKL_USE_MKL_VSL

} // namespace mckl

#endif // MCKL_RANDOM_MKL_HPP