MathUtils.h

Go to the documentation of this file.

#pragma once
 
 
#include "common/Assert.h"
#include "math/MathBasic.h"
#include <cmath>
#include <cstring>
#include <limits>
#include <utility>
 
NAMESPACE_SPH_BEGIN
 
template <typename T>
struct Eps;
 
template <>
struct Eps<float> {
    constexpr static float value = 1.e-6f;
};
template <>
struct Eps<double> {
    constexpr static double value = 1.e-12f;
};
 
constexpr Float EPS = Eps<Float>::value;
 
constexpr Float LARGE = 1.e20f;
 
constexpr Float INFTY = std::numeric_limits<Float>::max();
 
 
template <typename T>
INLINE T sqrtInv(const T& f) {
    int i;
    // cast to float
    float x2 = float(f) * 0.5f;
    float y = float(f);
    memcpy(&i, &y, sizeof(float));
    i = 0x5f3759df - (i >> 1);
    memcpy(&y, &i, sizeof(float));
    return y * (1.5f - (x2 * y * y));
}
 
template <typename T>
INLINE T sqrtApprox(const T f) {
    if (f == T(0)) {
        return T(0);
    }
    return T(1.f) / sqrtInv(f);
}
 
 
template <typename T>
INLINE constexpr T sqr(const T& f) noexcept {
    return f * f;
}
 
INLINE constexpr bool isPower2(const Size n) noexcept {
    return n >= 1 && (n & (n - 1)) == 0;
}
 
template <typename T>
INLINE T sqrt(const T f) {
    SPH_ASSERT(f >= T(0), f);
    return std::sqrt(f);
}
 
INLINE Float cbrt(const Float f) {
    return std::cbrt(f);
}
 
INLINE int positiveMod(const int i, const int n) {
    return (i % n + n) % n;
}
 
 
template <int N>
INLINE Float root(const Float f);
template <>
INLINE Float root<1>(const Float f) {
    return f;
}
template <>
INLINE Float root<2>(const Float f) {
    return sqrt(f);
}
template <>
INLINE Float root<3>(const Float f) {
    return cbrt(f);
}
template <>
INLINE Float root<4>(const Float f) {
    return sqrt(sqrt(f));
}
 
 
template <int N>
INLINE constexpr Float pow(const Float v);
 
template <>
INLINE constexpr Float pow<0>(const Float) {
    return 1._f;
}
template <>
INLINE constexpr Float pow<1>(const Float v) {
    return v;
}
template <>
INLINE constexpr Float pow<2>(const Float v) {
    return v * v;
}
template <>
INLINE constexpr Float pow<3>(const Float v) {
    return v * v * v;
}
template <>
INLINE constexpr Float pow<4>(const Float v) {
    return pow<2>(pow<2>(v));
}
template <>
INLINE constexpr Float pow<5>(const Float v) {
    return pow<2>(pow<2>(v)) * v;
}
template <>
INLINE constexpr Float pow<6>(const Float v) {
    return pow<3>(pow<2>(v));
}
template <>
INLINE constexpr Float pow<7>(const Float v) {
    return pow<6>(v) * v;
}
template <>
INLINE constexpr Float pow<8>(const Float v) {
    return pow<4>(pow<2>(v));
}
template <>
INLINE constexpr Float pow<-1>(const Float v) {
    return 1._f / v;
}
template <>
INLINE constexpr Float pow<-2>(const Float v) {
    return 1._f / (v * v);
}
template <>
INLINE constexpr Float pow<-3>(const Float v) {
    return 1._f / (v * v * v);
}
template <>
INLINE constexpr Float pow<-4>(const Float v) {
    return 1._f / (pow<2>(pow<2>(v)));
}
template <>
INLINE constexpr Float pow<-5>(const Float v) {
    return 1._f / (v * pow<2>(pow<2>(v)));
}
template <>
INLINE constexpr Float pow<-8>(const Float v) {
    return 1._f / (pow<2>(pow<4>(v)));
}
template <>
INLINE constexpr Float pow<-16>(const Float v) {
    return 1._f / (pow<2>(pow<8>(v)));
}
 
 
template <int N>
INLINE constexpr Size pow(const Size v);
 
template <>
INLINE constexpr Size pow<0>(const Size) {
    return 1;
}
template <>
INLINE constexpr Size pow<1>(const Size v) {
    return v;
}
template <>
INLINE constexpr Size pow<2>(const Size v) {
    return v * v;
}
template <>
INLINE constexpr Size pow<3>(const Size v) {
    return v * v * v;
}
template <>
INLINE constexpr Size pow<4>(const Size v) {
    return pow<2>(pow<2>(v));
}
template <>
INLINE constexpr Size pow<5>(const Size v) {
    return pow<2>(pow<2>(v)) * v;
}
template <>
INLINE constexpr Size pow<6>(const Size v) {
    return pow<3>(pow<2>(v));
}
 
 
 
template <typename T>
INLINE T pow(const T value, const T power) {
    return std::pow(value, power);
}
 
INLINE Float powFastest(const Float value, const Float power) {
    SPH_ASSERT(value > 0._f && power > 0._f, value, power);
    // https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
    // this type-punning through union is undefined behavior according to C++ standard, but most compilers
    // behave expectably ...
    union {
        double d;
        int x[2];
    } u = { value };
    u.x[1] = (int)(power * (u.x[1] - 1072632447) + 1072632447);
    u.x[0] = 0;
    return u.d;
}
 
INLINE Float powFast(Float value, const Float power) {
    SPH_ASSERT(value > 0._f && power > 0._f, value, power);
    // https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
    int e = (int)power;
    union {
        double d;
        int x[2];
    } u = { value };
    u.x[1] = (int)((power - e) * (u.x[1] - 1072632447) + 1072632447);
    u.x[0] = 0;
 
    double r = 1.0;
    while (e) {
        if (e & 1) {
            r *= value;
        }
        value *= value;
        e >>= 1;
    }
    return r * u.d;
}
 
template <typename T>
INLINE T exp(const T f) {
    return ::exp(f);
}
 
template <typename T>
INLINE auto abs(const T& f) {
    return ::fabs(f);
}
 
template <>
INLINE auto abs(const int& f) {
    return int(::abs(f));
}
 
template <>
INLINE auto abs(const Size& f) {
    return Size(f);
}
 
template <typename T>
INLINE T cos(const T f) {
    return std::cos(f);
}
 
template <typename T>
INLINE T sin(const T f) {
    return std::sin(f);
}
 
template <typename T>
INLINE T tan(const T f) {
    return std::tan(f);
}
 
template <typename T>
INLINE T acos(const T f) {
    return std::acos(f);
}
 
template <typename T>
INLINE T asin(const T f) {
    return std::asin(f);
}
 
template <typename T>
INLINE T atan(const T f) {
    return std::atan(f);
}
 
template <typename T>
INLINE T atan2(const T y, const T x) {
    return std::atan2(y, x);
}
 
template <typename T>
INLINE T exp10(const T f) {
    return std::pow(10._f, f);
}
 
template <typename T>
INLINE T log10(const T f) {
    return std::log10(f);
}
 
template <typename T>
INLINE int sgn(const T val) {
    return (T(0) < val) - (val < T(0));
}
 
template <typename T, typename TAmount>
INLINE T lerp(const T v1, const T v2, const TAmount amount) {
    return v1 * (TAmount(1) - amount) + v2 * amount;
}
 
template <typename T>
INLINE auto floor(const T& f) {
    return std::floor(f);
}
 
template <typename T>
INLINE auto ceil(const T& f) {
    return std::ceil(f);
}
 
template <typename T>
INLINE auto round(const T& f) {
    return std::round(f);
}
 
constexpr Float PI = 3.14159265358979323846264338327950288419716939937510582097_f;
 
constexpr Float PI_INV = 1._f / PI;
 
constexpr Float E = 2.718281828459045235360287471352662497757247093699959574967_f;
 
constexpr Float SQRT_3 = 1.732050807568877293527446341505872366942805253810380628055_f;
 
constexpr Float DEG_TO_RAD = PI / 180._f;
 
constexpr Float RAD_TO_DEG = 180._f / PI;
 
 
INLINE constexpr Float sphereVolume(const Float radius) {
    return 1.3333333333333333333333_f * PI * pow<3>(radius);
}
 
INLINE constexpr Float sphereSurfaceArea(const Float radius) {
    return 4._f * PI * pow<2>(radius);
}
 
INLINE bool almostEqual(const Float& f1, const Float& f2, const Float& eps = EPS) {
    return abs(f1 - f2) <= eps * (1._f + max(abs(f1), abs(f2)));
}
 
NAMESPACE_SPH_END