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/basic.hpp>
#include <vsmc/math/constants.hpp>

#if VSMC_USE_MKL_VML
#include <mkl_vml.h>
#elif VSMC_USE_ACCELERATE
#include <Accelerate/Accelerate.h>
#endif

#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)                 \
    {                                                                         \
        internal::size_check<MKL_INT>(n, #name);                              \
        ::vs##func(static_cast<MKL_INT>(n), a, y);                            \
    }                                                                         \
    inline void name(std::size_t n, const double *a, double *y)               \
    {                                                                         \
        internal::size_check<MKL_INT>(n, #name);                              \
        ::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) \
    {                                                                         \
        internal::size_check<MKL_INT>(n, #name);                              \
        ::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)           \
    {                                                                         \
        internal::size_check<MKL_INT>(n, #name);                              \
        ::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)
{
    internal::size_check<MKL_INT>(n, "linear_frac");
    ::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)
{
    internal::size_check<MKL_INT>(n, "linear_frac");
    ::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)
{
    internal::size_check<MKL_INT>(n, "pow");
    ::vsPowx(static_cast<MKL_INT>(n), a, b, y);
}
inline void pow(std::size_t n, const double *a, double b, double *y)
{
    internal::size_check<MKL_INT>(n, "pow");
    ::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)
{
    internal::size_check<MKL_INT>(n, "sincos");
    ::vsSinCos(static_cast<MKL_INT>(n), a, y, z);
}
inline void sincos(std::size_t n, const double *a, double *y, double *z)
{
    internal::size_check<MKL_INT>(n, "sincos");
    ::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)

VSMC_DEFINE_MATH_VMATH_VML_1(Floor, floor)
VSMC_DEFINE_MATH_VMATH_VML_1(Ceil, ceil)
VSMC_DEFINE_MATH_VMATH_VML_1(Trunc, trunc)
VSMC_DEFINE_MATH_VMATH_VML_1(Round, round)
inline void modf(std::size_t n, const float *a, float *y, float *z)
{
    internal::size_check<MKL_INT>(n, "modf");
    ::vsModf(static_cast<MKL_INT>(n), a, y, z);
}
inline void modf(std::size_t n, const double *a, double *y, double *z)
{
    internal::size_check<MKL_INT>(n, "modf");
    ::vdModf(static_cast<MKL_INT>(n), a, y, z);
}

} // namespace vsmc

#elif VSMC_USE_ACCELERATE

#define VSMC_DEFINE_MATH_VMATH_VFORCE_1(func, name)                           \
    inline void name(std::size_t n, const float *a, float *y)                 \
    {                                                                         \
        internal::size_check<int>(n, #name);                                  \
        int m = static_cast<int>(n);                                          \
        ::vv##func##f(y, a, &m);                                              \
    }                                                                         \
    inline void name(std::size_t n, const double *a, double *y)               \
    {                                                                         \
        internal::size_check<int>(n, #name);                                  \
        int m = static_cast<int>(n);                                          \
        ::vv##func(y, a, &m);                                                 \
    }

#define VSMC_DEFINE_MATH_VMATH_VFORCE_2(func, name)                           \
    inline void name(std::size_t n, const float *a, const float *b, float *y) \
    {                                                                         \
        internal::size_check<int>(n, #name);                                  \
        int m = static_cast<int>(n);                                          \
        ::vv##func##f(y, a, b, &m);                                           \
    }                                                                         \
    inline void name(                                                         \
        std::size_t n, const double *a, const double *b, double *y)           \
    {                                                                         \
        internal::size_check<int>(n, #name);                                  \
        int m = static_cast<int>(n);                                          \
        ::vv##func(y, a, b, &m);                                              \
    }

namespace vsmc
{

VSMC_DEFINE_MATH_VMATH_VFORCE_1(fabs, abs)

VSMC_DEFINE_MATH_VMATH_VFORCE_2(div, div)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(sqrt, sqrt)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(cbrt, cbrt)
VSMC_DEFINE_MATH_VMATH_VFORCE_2(pow, pow)

VSMC_DEFINE_MATH_VMATH_VFORCE_1(exp, exp)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(expm1, expm1)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(log, log)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(log10, log10)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(log1p, log1p)

VSMC_DEFINE_MATH_VMATH_VFORCE_1(cos, cos)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(sin, sin)
inline void sincos(std::size_t n, const float *a, float *y, float *z)
{
    internal::size_check<int>(n, "sincos");
    const int m = static_cast<int>(n);
    ::vvsincosf(y, z, a, &m);
}
inline void sincos(std::size_t n, const double *a, double *y, double *z)
{
    internal::size_check<int>(n, "sincos");
    const int m = static_cast<int>(n);
    ::vvsincos(y, z, a, &m);
}
VSMC_DEFINE_MATH_VMATH_VFORCE_1(tan, tan)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(acos, acos)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(asin, asin)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(atan, atan)
VSMC_DEFINE_MATH_VMATH_VFORCE_2(atan2, atan2)

VSMC_DEFINE_MATH_VMATH_VFORCE_1(cosh, cosh)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(sinh, sinh)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(tanh, tanh)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(acosh, acosh)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(asinh, asinh)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(atanh, atanh)

VSMC_DEFINE_MATH_VMATH_VFORCE_1(floor, floor)
VSMC_DEFINE_MATH_VMATH_VFORCE_1(ceil, ceil)

} // namespace vsmc

#endif // VSMC_USE_MKL_VML

#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];                                                 \
    }

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 = internal::BufferSize<T>::value;
    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)
{
    if (internal::is_zero(c)) {
        mul(n, a, b, y);
    } else {
        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)
{
    if (internal::is_one(b)) {
        add(n, a, c, y);
    } else {
        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)
{
    if (internal::is_one(b) && internal::is_zero(c)) {
        std::copy_n(a, n, y);
    } else if (internal::is_one(b)) {
        add(n, a, c, y);
    } else if (internal::is_zero(c)) {
        mul(n, a, b, y);
    } else {
        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)
{
    if (internal::is_one(a)) {
        add(n, b, c, y);
    } else {
        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)
{
    if (internal::is_zero(c)) {
        mul(n, a, b, y);
    } else {
        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)
{
    for (std::size_t i = 0; i != n; ++i)
        y[i] = static_cast<T>(1) / 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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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 = internal::BufferSize<T>::value;
    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, static_cast<T>(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, static_cast<T>(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)

VSMC_DEFINE_MATH_VMATH_1(std::floor, floor)

VSMC_DEFINE_MATH_VMATH_1(std::ceil, ceil)

VSMC_DEFINE_MATH_VMATH_1(std::trunc, trunc)

VSMC_DEFINE_MATH_VMATH_1(std::round, round)

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


} // namespace vsmc

#endif // VSMC_MATH_VMATH_HPP