vmath.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/math/vmath.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013-2016, 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 VSMC_MATH_VMATH_HPP
#define VSMC_MATH_VMATH_HPP

#include <vsmc/internal/config.h>
#include <vsmc/math/constants.hpp>
#include <cmath>

#if VSMC_USE_MKL_VML
#include <mkl.h>
#endif

#define VSMC_DEFINE_MATH_VMATH_1(func, name)                                  \
    template <typename T>                                                     \
    inline void name(std::size_t n, const T *a, T *y)                         \
    {                                                                         \
        for (std::size_t i = 0; i != n; ++i)                                  \
            y[i] = func(a[i]);                                                \
    }

#define VSMC_DEFINE_MATH_VMATH_2(func, name)                                  \
    template <typename T>                                                     \
    inline void name(std::size_t n, const T *a, const T *b, T *y)             \
    {                                                                         \
        for (std::size_t i = 0; i != n; ++i)                                  \
            y[i] = func(a[i], b[i]);                                          \
    }

#define VSMC_DEFINE_MATH_VMATH_B(op, name)                                    \
    template <typename T>                                                     \
    inline void name(std::size_t n, const T *a, const T *b, T *y)             \
    {                                                                         \
        for (std::size_t i = 0; i != n; ++i)                                  \
            y[i] = a[i] op b[i];                                              \
    }

#define VSMC_DEFINE_MATH_VMATH_VS(op, name)                                   \
    template <typename T>                                                     \
    inline void name(std::size_t n, const T *a, T b, T *y)                    \
    {                                                                         \
        for (std::size_t i = 0; i != n; ++i)                                  \
            y[i] = a[i] op b;                                                 \
    }

#define VSMC_DEFINE_MATH_VMATH_SV(op, name)                                   \
    template <typename T>                                                     \
    inline void name(std::size_t n, T a, const T *b, T *y)                    \
    {                                                                         \
        for (std::size_t i = 0; i != n; ++i)                                  \
            y[i] = a op b[i];                                                 \
    }

#if VSMC_USE_MKL_VML

#define VSMC_DEFINE_MATH_VMATH_VML_1(func, name)                              \
    inline void name(std::size_t n, const float *a, float *y)                 \
    {                                                                         \
        ::vs##func(static_cast<MKL_INT>(n), a, y);                            \
    }                                                                         \
    inline void name(std::size_t n, const double *a, double *y)               \
    {                                                                         \
        ::vd##func(static_cast<MKL_INT>(n), a, y);                            \
    }

#define VSMC_DEFINE_MATH_VMATH_VML_2(func, name)                              \
    inline void name(std::size_t n, const float *a, const float *b, float *y) \
    {                                                                         \
        ::vs##func(static_cast<MKL_INT>(n), a, b, y);                         \
    }                                                                         \
    inline void name(                                                         \
        std::size_t n, const double *a, const double *b, double *y)           \
    {                                                                         \
        ::vd##func(static_cast<MKL_INT>(n), a, b, y);                         \
    }

namespace vsmc
{

VSMC_DEFINE_MATH_VMATH_VML_2(Add, add)
VSMC_DEFINE_MATH_VMATH_VML_2(Sub, sub)
VSMC_DEFINE_MATH_VMATH_VML_1(Sqr, sqr)
VSMC_DEFINE_MATH_VMATH_VML_2(Mul, mul)
VSMC_DEFINE_MATH_VMATH_VML_1(Abs, abs)
inline void linear_frac(std::size_t n, const float *a, const float *b,
    float beta_a, float beta_b, float mu_a, float mu_b, float *y)
{
    ::vsLinearFrac(
        static_cast<MKL_INT>(n), a, b, beta_a, beta_b, mu_a, mu_b, y);
}
inline void linear_frac(std::size_t n, const double *a, const double *b,
    double beta_a, double beta_b, double mu_a, double mu_b, double *y)
{
    ::vdLinearFrac(
        static_cast<MKL_INT>(n), a, b, beta_a, beta_b, mu_a, mu_b, y);
}

VSMC_DEFINE_MATH_VMATH_VML_1(Inv, inv)
VSMC_DEFINE_MATH_VMATH_VML_2(Div, div)
VSMC_DEFINE_MATH_VMATH_VML_1(Sqrt, sqrt)
VSMC_DEFINE_MATH_VMATH_VML_1(InvSqrt, invsqrt)
VSMC_DEFINE_MATH_VMATH_VML_1(Cbrt, cbrt)
VSMC_DEFINE_MATH_VMATH_VML_1(InvCbrt, invcbrt)
VSMC_DEFINE_MATH_VMATH_VML_1(Pow2o3, pow2o3)
VSMC_DEFINE_MATH_VMATH_VML_1(Pow3o2, pow3o2)
VSMC_DEFINE_MATH_VMATH_VML_2(Pow, pow)
inline void pow(std::size_t n, const float *a, float b, float *y)
{
    ::vsPowx(static_cast<MKL_INT>(n), a, b, y);
}
inline void pow(std::size_t n, const double *a, double b, double *y)
{
    ::vdPowx(static_cast<MKL_INT>(n), a, b, y);
}
VSMC_DEFINE_MATH_VMATH_VML_2(Hypot, hypot)

VSMC_DEFINE_MATH_VMATH_VML_1(Exp, exp)
VSMC_DEFINE_MATH_VMATH_VML_1(Expm1, expm1)
VSMC_DEFINE_MATH_VMATH_VML_1(Ln, log)
VSMC_DEFINE_MATH_VMATH_VML_1(Log10, log10)
VSMC_DEFINE_MATH_VMATH_VML_1(Log1p, log1p)

VSMC_DEFINE_MATH_VMATH_VML_1(Cos, cos)
VSMC_DEFINE_MATH_VMATH_VML_1(Sin, sin)
inline void sincos(std::size_t n, const float *a, float *y, float *z)
{
    ::vsSinCos(static_cast<MKL_INT>(n), a, y, z);
}
inline void sincos(std::size_t n, const double *a, double *y, double *z)
{
    ::vdSinCos(static_cast<MKL_INT>(n), a, y, z);
}
VSMC_DEFINE_MATH_VMATH_VML_1(Tan, tan)
VSMC_DEFINE_MATH_VMATH_VML_1(Acos, acos)
VSMC_DEFINE_MATH_VMATH_VML_1(Asin, asin)
VSMC_DEFINE_MATH_VMATH_VML_1(Atan, atan)
VSMC_DEFINE_MATH_VMATH_VML_2(Atan2, atan2)

VSMC_DEFINE_MATH_VMATH_VML_1(Cosh, cosh)
VSMC_DEFINE_MATH_VMATH_VML_1(Sinh, sinh)
VSMC_DEFINE_MATH_VMATH_VML_1(Tanh, tanh)
VSMC_DEFINE_MATH_VMATH_VML_1(Acosh, acosh)
VSMC_DEFINE_MATH_VMATH_VML_1(Asinh, asinh)
VSMC_DEFINE_MATH_VMATH_VML_1(Atanh, atanh)

VSMC_DEFINE_MATH_VMATH_VML_1(Erf, erf)
VSMC_DEFINE_MATH_VMATH_VML_1(Erfc, erfc)
VSMC_DEFINE_MATH_VMATH_VML_1(CdfNorm, cdfnorm)
VSMC_DEFINE_MATH_VMATH_VML_1(ErfInv, erfinv)
VSMC_DEFINE_MATH_VMATH_VML_1(ErfcInv, erfcinv)
VSMC_DEFINE_MATH_VMATH_VML_1(CdfNormInv, cdfnorminv)
VSMC_DEFINE_MATH_VMATH_VML_1(LGamma, lgamma)
VSMC_DEFINE_MATH_VMATH_VML_1(TGamma, tgamm)

} // namespace vsmc

#endif // VSMC_USE_MKL_VML

namespace vsmc
{


VSMC_DEFINE_MATH_VMATH_B(+, add)

VSMC_DEFINE_MATH_VMATH_VS(+, add)

VSMC_DEFINE_MATH_VMATH_SV(+, add)

VSMC_DEFINE_MATH_VMATH_B(-, sub)

VSMC_DEFINE_MATH_VMATH_VS(-, sub)

VSMC_DEFINE_MATH_VMATH_SV(-, sub)

template <typename T>
inline void sqr(std::size_t n, const T *a, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a[i] * a[i];
}

VSMC_DEFINE_MATH_VMATH_B(*, mul)

VSMC_DEFINE_MATH_VMATH_VS(*, mul)

VSMC_DEFINE_MATH_VMATH_SV(*, mul)

VSMC_DEFINE_MATH_VMATH_1(std::abs, abs)

template <typename T>
inline void linear_frac(std::size_t n, const T *a, const T *b, T beta_a,
    T beta_b, T mu_a, T mu_b, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        for (std::size_t j = 0; j != k; ++j)
            y[j] = beta_a * a[j] + mu_a;
        for (std::size_t j = 0; j != k; ++j)
            y[j] /= beta_b * b[j] + mu_b;
    }
    for (std::size_t i = 0; i != l; ++i)
        y[i] = beta_a * a[i] + mu_a;
    for (std::size_t i = 0; i != l; ++i)
        y[i] /= beta_b * b[i] + mu_b;
}

template <typename T>
inline void fma(std::size_t n, const T *a, const T *b, const T *c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a[i] * b[i] + c[i];
}

template <typename T>
inline void fma(std::size_t n, const T *a, const T *b, T c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a[i] * b[i] + c;
}

template <typename T>
inline void fma(std::size_t n, const T *a, T b, const T *c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a[i] * b + c[i];
}

template <typename T>
inline void fma(std::size_t n, const T *a, T b, T c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a[i] * b + c;
}

template <typename T>
inline void fma(std::size_t n, T a, const T *b, const T *c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a * b[i] + c[i];
}

template <typename T>
inline void fma(std::size_t n, T a, const T *b, T c, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = a * b[i] + c;
}

template <typename T>
inline void fma(std::size_t n, T a, T b, const T *c, T *y)
{
    add(n, a * b, c, y);
}



template <typename T>
inline void inv(std::size_t n, const T *a, T *y)
{
    const T one = static_cast<T>(1);
    for (std::size_t i = 0; i != n; ++i)
        y[i] = one / a[i];
}

VSMC_DEFINE_MATH_VMATH_B(/, div)

VSMC_DEFINE_MATH_VMATH_VS(/, div)

VSMC_DEFINE_MATH_VMATH_SV(/, div)

VSMC_DEFINE_MATH_VMATH_1(std::sqrt, sqrt)

template <typename T>
inline void invsqrt(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        sqrt(k, a, y);
        inv(k, y, y);
    }
    sqrt(l, a, y);
    inv(l, y, y);
}

VSMC_DEFINE_MATH_VMATH_1(std::cbrt, cbrt)

template <typename T>
inline void invcbrt(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        cbrt(k, a, y);
        inv(k, y, y);
    }
    cbrt(l, a, y);
    inv(l, y, y);
}

template <typename T>
inline void pow2o3(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        cbrt(k, a, y);
        sqr(k, y, y);
    }
    cbrt(l, a, y);
    sqr(l, y, y);
}

template <typename T>
inline void pow3o2(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        sqrt(k, a, y);
        for (std::size_t j = 0; j != k; ++j)
            y[j] = y[j] * y[j] * y[j];
    }
    sqrt(l, a, y);
    for (std::size_t i = 0; i != l; ++i)
        y[i] = y[i] * y[i] * y[i];
}

VSMC_DEFINE_MATH_VMATH_2(std::pow, pow)

template <typename T>
inline void pow(std::size_t n, const T *a, T b, T *y)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = std::pow(a[i], b);
}

VSMC_DEFINE_MATH_VMATH_2(std::hypot, hypot)



VSMC_DEFINE_MATH_VMATH_1(std::exp, exp)

VSMC_DEFINE_MATH_VMATH_1(std::exp2, exp2)

template <typename T>
inline void exp10(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        mul(k, const_ln_10<T>(), a, y);
        exp(k, y, y);
    }
    mul(l, const_ln_10<T>(), a, y);
    exp(l, y, y);
}

VSMC_DEFINE_MATH_VMATH_1(std::expm1, expm1)

VSMC_DEFINE_MATH_VMATH_1(std::log, log)

VSMC_DEFINE_MATH_VMATH_1(std::log2, log2)

VSMC_DEFINE_MATH_VMATH_1(std::log10, log10)

VSMC_DEFINE_MATH_VMATH_1(std::log1p, log1p)


VSMC_DEFINE_MATH_VMATH_1(std::cos, cos)

VSMC_DEFINE_MATH_VMATH_1(std::sin, sin)

template <typename T>
inline void sincos(std::size_t n, const T *a, T *y, T *z)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k, z += k) {
        sin(k, a, y);
        cos(k, a, z);
    }
    sin(l, a, y);
    cos(l, a, z);
}

VSMC_DEFINE_MATH_VMATH_1(std::tan, tan)

VSMC_DEFINE_MATH_VMATH_1(std::acos, acos)

VSMC_DEFINE_MATH_VMATH_1(std::asin, asin)

VSMC_DEFINE_MATH_VMATH_1(std::atan, atan)

VSMC_DEFINE_MATH_VMATH_2(std::atan2, atan2)



VSMC_DEFINE_MATH_VMATH_1(std::cosh, cosh)

VSMC_DEFINE_MATH_VMATH_1(std::sinh, sinh)

VSMC_DEFINE_MATH_VMATH_1(std::tanh, tanh)

VSMC_DEFINE_MATH_VMATH_1(std::acosh, acosh)

VSMC_DEFINE_MATH_VMATH_1(std::asinh, asinh)

VSMC_DEFINE_MATH_VMATH_1(std::atanh, atanh)



VSMC_DEFINE_MATH_VMATH_1(std::erf, erf)

VSMC_DEFINE_MATH_VMATH_1(std::erfc, erfc)

template <typename T>
inline void cdfnorm(std::size_t n, const T *a, T *y)
{
    const std::size_t k = 1024;
    const std::size_t m = n / k;
    const std::size_t l = n % k;
    for (std::size_t i = 0; i != m; ++i, a += k, y += k) {
        mul(k, 1 / const_sqrt_2<T>(), a, y);
        erf(k, y, y);
        fma(k, static_cast<T>(0.5), static_cast<T>(0.5), y, y);
    }
    mul(l, 1 / const_sqrt_2<T>(), a, y);
    erf(l, y, y);
    fma(l, static_cast<T>(0.5), static_cast<T>(0.5), y, y);
}

VSMC_DEFINE_MATH_VMATH_1(std::lgamma, lgamma)

VSMC_DEFINE_MATH_VMATH_1(std::tgamma, tgamma)

} // namespace vsmc

#endif // VSMC_MATH_VMATH_HPP