math.hpp

浏览该文件的文档.

#ifndef NEFORCE_CORE_NUMERIC_MATH_HPP__
#define NEFORCE_CORE_NUMERIC_MATH_HPP__

 
#include "NeForce/core/exception/exception.hpp"
#include "NeForce/core/numeric/numeric_types.hpp"
NEFORCE_BEGIN_NAMESPACE__
 
NEFORCE_BEGIN_CONSTANTS__

 
NEFORCE_INLINE17 constexpr decimal_t EULER = 2.71828182845904523536L; 
NEFORCE_INLINE17 constexpr decimal_t PI = 3.14159265358979323846L;    
NEFORCE_INLINE17 constexpr decimal_t PHI = 1.61803398874989484820L;   
NEFORCE_INLINE17 constexpr decimal_t SEMI_CIRCLE = 180.0L;            
NEFORCE_INLINE17 constexpr decimal_t CIRCLE = 360.0L;                 

NEFORCE_INLINE17 constexpr decimal_t TWO_PI_HI = 6.28318530717958647692L;
NEFORCE_INLINE17 constexpr decimal_t TWO_PI_LO = 2.44929359829470641435e-16L;

NEFORCE_INLINE17 constexpr decimal_t MACHINE_EPSILON = numeric_traits<decimal_t>::epsilon();
NEFORCE_INLINE17 constexpr decimal_t DEFAULT_TOLERANCE = 1e-12L;
NEFORCE_INLINE17 constexpr decimal_t LOOSE_TOLERANCE = 1e-9L;

NEFORCE_INLINE17 constexpr uint64_t FIBONACCI_LIST[] = {
        0,          1,          1,          2,          3,         5,         8,         13,        21,
        34,         55,         89,         144,        233,       377,       610,       987,       1597,
        2584,       4181,       6765,       10946,      17711,     28657,     46368,     75025,     121393,
        196418,     317811,     514229,     832040,     1346269,   2178309,   3524578,   5702887,   9227465,
        14930352,   24157817,   39088169,   63245986,   102334155, 165580141, 267914296, 433494437, 701408733,
        1134903170, 1836311903, 2971215073, 4807526976, 7778742049};
 
NEFORCE_INLINE17 constexpr uint32_t FIBONACCI_COUNT = extent_v<decltype(FIBONACCI_LIST)>; 
 // MathConstants
 
NEFORCE_END_CONSTANTS__

 
NEFORCE_CONST_FUNCTION constexpr int64_t safe_trunc(const decimal_t x) noexcept {
    constexpr decimal_t max_int64 = static_cast<decimal_t>(numeric_traits<int64_t>::max());
    constexpr decimal_t min_int64 = static_cast<decimal_t>(numeric_traits<int64_t>::min());
    if (is_nan(x) || is_infinity(x) || x > max_int64 || x < min_int64) {
        return 0;
    }
    return static_cast<int64_t>(x);
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 uint64_t fibonacci(const uint32_t n) noexcept {
    if (n < constants::FIBONACCI_COUNT) {
        return constants::FIBONACCI_LIST[n];
    }
    uint64_t a = constants::FIBONACCI_LIST[constants::FIBONACCI_COUNT - 2];
    uint64_t b = constants::FIBONACCI_LIST[constants::FIBONACCI_COUNT - 1];
    for (uint32_t i = constants::FIBONACCI_COUNT; i <= n; ++i) {
        const uint64_t c = a + b;
        a = b;
        b = c;
    }
    return b;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 uint64_t leonardo(const uint32_t n) noexcept {
    return 2 * fibonacci(n + 1) - 1;
}

template <typename T>
NEFORCE_PURE_FUNCTION constexpr T angular2radian(const T angular) noexcept {
    static_assert(is_arithmetic_v<T>, "arithmetic required");
    return angular * static_cast<T>(constants::PI) / static_cast<T>(constants::SEMI_CIRCLE);
}

template <typename T>
NEFORCE_PURE_FUNCTION constexpr T radian2angular(const T radian) noexcept {
    static_assert(is_arithmetic_v<T>, "arithmetic required");
    return radian * (static_cast<T>(constants::SEMI_CIRCLE) / static_cast<T>(constants::PI));
}

template <typename T>
NEFORCE_CONST_FUNCTION constexpr enable_if_t<is_signed_v<T>, T> absolute(const T x) noexcept {
    return x > T(0) ? x : -x;
}

 
template <typename T>
NEFORCE_CONST_FUNCTION constexpr enable_if_t<is_unsigned_v<T>, T> absolute(const T x) noexcept {
    return x;
}

template <typename T>
NEFORCE_CONST_FUNCTION constexpr const T& sum(const T& x) noexcept {
    return x;
}

template <typename First, typename... Rests, enable_if_t<(sizeof...(Rests) > 0), int> = 0>
NEFORCE_CONST_FUNCTION constexpr decltype(auto) sum(First first, Rests... args) {
    return first + _NEFORCE sum(args...);
}

template <typename... Args, enable_if_t<(sizeof...(Args) > 0), int> = 0>
NEFORCE_CONST_FUNCTION constexpr decltype(auto) average(Args... args) {
    return _NEFORCE sum(args...) / sizeof...(Args);
}

template <typename T>
NEFORCE_CONSTEXPR14 int sign(const T& value) noexcept {
    static_assert(is_arithmetic_v<T>, "arithmetic required");
    constexpr T zero = T(0);
    if (value > zero) {
        return 1;
    }
    if (value < zero) {
        return -1;
    }
    return 0;
}

template <typename T>
NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 T gcd(const T& m, const T& n) noexcept {
    T x = _NEFORCE absolute(m), y = _NEFORCE absolute(n);
    constexpr T zero = T(0);
    while (y != zero) {
        T t = x % y;
        x = y;
        y = t;
    }
    return x;
}

template <typename T>
NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 T lcm(const T& m, const T& n) noexcept {
    return (m / _NEFORCE gcd(m, n)) * n;
}

template <typename T>
NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 enable_if_t<is_floating_point_v<T>, T> mod(const T x, const T y) {
    if (_NEFORCE is_nan(x) || _NEFORCE is_nan(y)) {
        return numeric_traits<T>::quiet_nan();
    }
    if (y == 0) {
        NEFORCE_THROW_EXCEPTION(math_exception("zero can not be dividend."));
    }
    if (_NEFORCE is_infinity(x) || _NEFORCE is_infinity(y)) {
        return numeric_traits<T>::quiet_nan();
    }
    return x - static_cast<make_integer_t<sizeof(T)>>(x / y) * y;
}

template <typename T>
NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 enable_if_t<is_integral_v<T>, T> mod(const T x, const T y) {
    if (y == 0) {
        NEFORCE_THROW_EXCEPTION(math_exception("zero can not be dividend."));
    }
    return x % y;
}

template <typename T>
NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 T power(const T& x, uint32_t n) noexcept {
    static_assert(is_arithmetic_v<T>, "arithmetic required");
    if (n == 0) {
        return 1;
    }
    T result(1);
    T base = x;
    while (n > 0) {
        if (n % 2 == 1) {
            result *= base;
        }
        base *= base;
        n /= 2;
    }
    return result;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t exponential(const uint32_t n) noexcept {
    return power(constants::EULER, n);
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t logarithm_e(const decimal_t x) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (x < 0.0L) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
    if (x == 0.0L) {
        return -numeric_traits<decimal_t>::infinity();
    }
    if (is_infinity(x) && x > 0) {
        return x;
    }
 
    int64_t k = 0;
    decimal_t m = x;
 
    if (m >= 2.0L) {
        while (m >= 2.0L) {
            m *= 0.5L;
            ++k;
        }
    } else if (m < 1.0L) {
        while (m < 1.0L) {
            m *= 2.0L;
            --k;
        }
    }
 
    const decimal_t a = (m - 1.0L) / (m + 1.0L);
    const decimal_t a2 = a * a;
    decimal_t term = a;
    decimal_t s = a;
    decimal_t n = 1.0L;
 
    while (absolute(term) > constants::MACHINE_EPSILON * absolute(s)) {
        term *= a2;
        n += 2.0L;
        s += term / n;
    }
 
    constexpr decimal_t LN2 = 0.69314718055994530941723212145817656807L;
    return 2.0L * s + static_cast<decimal_t>(k) * LN2;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t logarithm(const decimal_t x, const uint32_t base) {
    const decimal_t under = logarithm_e(base);
    if (under == 0.0L) {
        NEFORCE_THROW_EXCEPTION(math_exception("zero can not be dividend."));
    }
    return logarithm_e(x) / under;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t logarithm_2(const decimal_t x) { return logarithm(x, 2); }

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t logarithm_10(const decimal_t x) { return logarithm(x, 10); }

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t
square_root(const decimal_t x, const decimal_t precise = constants::DEFAULT_TOLERANCE) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (x < 0.0L) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
    if (x == 0.0L) {
        return x;
    }
    if (is_infinity(x) && x > 0) {
        return x;
    }
 
    decimal_t guess = x * 0.5L;
    decimal_t prev;
    do {
        prev = guess;
        guess = 0.5L * (prev + x / prev);
    } while (absolute(guess - prev) > precise * absolute(guess));
    return guess;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t
cube_root(const decimal_t x, const decimal_t precise = constants::DEFAULT_TOLERANCE) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (is_infinity(x)) {
        return x;
    }
    if (x == 0.0L) {
        return x;
    }
 
    const bool negative = x < 0.0L;
    const decimal_t v = negative ? -x : x;
 
    decimal_t guess = v > 1.0L ? v / 3.0L : v;
    decimal_t prev;
    do {
        prev = guess;
        guess = (2.0L * prev + v / (prev * prev)) / 3.0L;
    } while (absolute(guess - prev) > precise * absolute(guess));
 
    return negative ? -guess : guess;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 uint64_t factorial(const uint32_t n) noexcept {
    uint64_t result = 1;
    for (uint32_t i = 2; i <= n; ++i) {
        result *= i;
    }
    return result;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t floor(const decimal_t x) noexcept {
    if (!is_finite(x)) {
        return x;
    }
    const int64_t i = safe_trunc(x);
    const decimal_t di = static_cast<decimal_t>(i);
    return (x < 0.0L && x != di) ? di - 1.0L : di;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t floor(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t factor = power(10.0L, bit);
    return floor(x * factor) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t ceil(const decimal_t x) noexcept {
    if (!is_finite(x)) {
        return x;
    }
    const int64_t i = safe_trunc(x);
    const decimal_t di = static_cast<decimal_t>(i);
    return (x > 0.0L && x != di) ? di + 1.0L : di;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t ceil(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t factor = power(10.0L, bit);
    return ceil(x * factor) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t round(const decimal_t x) noexcept {
    if (!is_finite(x)) {
        return x;
    }
 
    const decimal_t int_part = floor(absolute(x));
    const decimal_t frac = absolute(x) - int_part;
    if (frac < 0.5L) {
        return x < 0 ? -int_part : int_part;
    }
    if (frac > 0.5L) {
        return x < 0 ? -(int_part + 1) : (int_part + 1);
    }
 
    if (static_cast<int64_t>(int_part) % 2 == 0) {
        return x < 0 ? -int_part : int_part;
    } else {
        return x < 0 ? -(int_part + 1) : (int_part + 1);
    }
}
 

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t floor_bit(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t factor = power(10.0L, bit);
    return floor(x * factor) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t ceil_bit(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t factor = power(10.0L, bit);
    return ceil(x * factor) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t round_bit(const decimal_t x, const uint32_t bit) noexcept {
    return x < 0 ? ceil_bit(x - 0.5 / power(10.0, bit), bit) : floor_bit(x + 0.5 / power(10.0, bit), bit);
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t truncate_bit(const decimal_t x, const uint32_t bit) noexcept {
    return x < 0 ? ceil_bit(x, bit) : floor_bit(x, bit);
}
 

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t round(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t factor = power(10.0L, bit);
    return round(x * factor) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t truncate(const decimal_t x, const int bit) noexcept {
    const decimal_t factor = power(10.0L, static_cast<uint32_t>(absolute(bit)));
    return (x < 0 ? ceil(x * factor, 0) : floor(x * factor, 0)) / factor;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t truncate(const decimal_t x) noexcept {
    if (!is_finite(x)) {
        return x;
    }
    return static_cast<decimal_t>(safe_trunc(x));
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 bool around_multiple(const decimal_t x, const decimal_t axis,
                                                                const decimal_t toler = constants::DEFAULT_TOLERANCE) {
    if (absolute(axis) < constants::MACHINE_EPSILON) {
        NEFORCE_THROW_EXCEPTION(math_exception("Axis Cannot be 0"));
    }
    const decimal_t multi = round(x / axis) * axis;
    return absolute(x - multi) < toler;
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 bool around_pi(const decimal_t x,
                                                          const decimal_t toler = constants::LOOSE_TOLERANCE) {
    return around_multiple(x, constants::PI, toler);
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 bool
around_zero(const decimal_t x, const decimal_t toler = constants::LOOSE_TOLERANCE) noexcept {
    return absolute(x) < toler;
}
 

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t remainder(const decimal_t x, const decimal_t y) noexcept {
    if (is_nan(x) || is_nan(y) || is_infinity(x) || y == 0.0L) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
    return x - y * round(x / y);
}

NEFORCE_CONST_FUNCTION constexpr decimal_t float_part(const decimal_t x) noexcept {
    return x - static_cast<decimal_t>(safe_trunc(x));
}

NEFORCE_CONST_FUNCTION NEFORCE_CONSTEXPR14 decimal_t float_apart(decimal_t x, int64_t* int_ptr) noexcept {
    *int_ptr = safe_trunc(x);
    return x - static_cast<decimal_t>(*int_ptr);
}

NEFORCE_BEGIN_INNER__
 
NEFORCE_CONSTEXPR14 void reduce_arg_sincos(decimal_t& x, int& quadrant) noexcept {
    if (x < 0.0L) {
        x = -x;
        quadrant = 3;
    } else {
        quadrant = 0;
    }
 
    constexpr decimal_t HALF_PI_HI = 1.57079632679489661923L;
    constexpr decimal_t HALF_PI_LO = 6.12323399573676603587e-17L;
    constexpr decimal_t INV_HALF_PI = 0.63661977236758134308L;
 
    const decimal_t kd = round(x * INV_HALF_PI);
    const int64_t k = static_cast<int64_t>(kd);
    quadrant = (quadrant + (k & 3)) & 3;
 
    decimal_t t = kd * HALF_PI_HI;
    x = x - t;
    t = kd * HALF_PI_LO;
    x = x - t;
}
 
NEFORCE_CONSTEXPR14 void sincos_taylor(const decimal_t x, decimal_t& sin_val, decimal_t& cos_val) noexcept {
    const decimal_t x2 = x * x;
 
    decimal_t term_sin = x;
    decimal_t sum_sin = x;
 
    decimal_t term_cos = 1.0L;
    decimal_t sum_cos = 1.0L;
 
    decimal_t i = 1.0L;
    bool sin_done = false, cos_done = false;
 
    while (!sin_done || !cos_done) {
        if (!sin_done) {
            term_sin = -term_sin * x2 / ((i + 1.0L) * (i + 2.0L));
            sum_sin += term_sin;
            if (absolute(term_sin) <= constants::MACHINE_EPSILON * absolute(sum_sin)) {
                sin_done = true;
            }
        }
        if (!cos_done) {
            term_cos = -term_cos * x2 / (i * (i + 1.0L));
            sum_cos += term_cos;
            if (absolute(term_cos) <= constants::MACHINE_EPSILON * absolute(sum_cos)) {
                cos_done = true;
            }
        }
        i += 2.0L;
    }
 
    sin_val = sum_sin;
    cos_val = sum_cos;
}
 
NEFORCE_END_INNER__

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t sine(decimal_t x) noexcept {
    if (is_nan(x) || is_infinity(x)) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
 
    int quadrant;
    inner::reduce_arg_sincos(x, quadrant);
 
    decimal_t sin_x, cos_x;
    inner::sincos_taylor(x, sin_x, cos_x);
 
    switch (quadrant) {
        case 0:
            return sin_x;
        case 1:
            return cos_x;
        case 2:
            return -sin_x;
        default:
            return -cos_x;
    }
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t cosine(decimal_t x) noexcept {
    if (is_nan(x) || is_infinity(x)) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
 
    int quadrant;
    inner::reduce_arg_sincos(x, quadrant);
    decimal_t sin_x, cos_x;
    inner::sincos_taylor(x, sin_x, cos_x);
 
    switch (quadrant) {
        case 0:
            return cos_x;
        case 1:
            return -sin_x;
        case 2:
            return -cos_x;
        default:
            return sin_x;
    }
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t tangent(decimal_t x) noexcept {
    if (is_nan(x) || is_infinity(x)) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
 
    int quadrant;
    inner::reduce_arg_sincos(x, quadrant);
    decimal_t sin_x, cos_x;
    inner::sincos_taylor(x, sin_x, cos_x);
 
    decimal_t s, c;
    switch (quadrant) {
        case 0:
            s = sin_x;
            c = cos_x;
            break;
        case 1:
            s = cos_x;
            c = -sin_x;
            break;
        case 2:
            s = -sin_x;
            c = -cos_x;
            break;
        default:
            s = -cos_x;
            c = sin_x;
            break;
    }
 
    if (around_zero(c, constants::LOOSE_TOLERANCE)) {
        constexpr auto inf = numeric_traits<decimal_t>::infinity();
        return (s > 0) ? inf : -inf;
    }
    return s / c;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t cotangent(const decimal_t x) { return 1.0L / tangent(x); }

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t arctangent(decimal_t x) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (is_infinity(x)) {
        return (x > 0) ? constants::PI / 2 : -constants::PI / 2;
    }
 
    bool negative = x < 0.0L;
    if (negative) {
        x = -x;
    }
 
    bool recip = false;
    if (x > 1.0L) {
        x = 1.0L / x;
        recip = true;
    }
 
    x = x / (1.0L + square_root(1.0L + x * x));
    x = x / (1.0L + square_root(1.0L + x * x));
 
    const decimal_t x2 = x * x;
    decimal_t term = x, s = x, n = 1.0L;
    while (absolute(term) > constants::MACHINE_EPSILON * absolute(s)) {
        term *= -x2;
        n += 2.0L;
        s += term / n;
    }
 
    decimal_t result = 4.0L * s;
    if (recip) {
        result = constants::PI / 2 - result;
    }
    if (negative) {
        result = -result;
    }
    return result;
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t arcsine(const decimal_t x) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (absolute(x) > 1.0L) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
    if (x == 1.0L) {
        return constants::PI / 2;
    }
    if (x == -1.0L) {
        return -constants::PI / 2;
    }
    return arctangent(x / square_root(1.0L - x * x));
}

NEFORCE_PURE_FUNCTION NEFORCE_CONSTEXPR14 decimal_t arccosine(const decimal_t x) noexcept {
    if (is_nan(x)) {
        return x;
    }
    if (absolute(x) > 1.0L) {
        return numeric_traits<decimal_t>::quiet_nan();
    }
    if (x == 1.0L) {
        return 0.0L;
    }
    if (x == -1.0L) {
        return constants::PI;
    }
    return constants::PI / 2 - arcsine(x);
}
 // MathFunctions
 
NEFORCE_END_NAMESPACE__
#endif // NEFORCE_CORE_NUMERIC_MATH_HPP__