math.hpp

浏览该文件的文档.

#ifndef MSTL_CORE_NUMERIC_MATH_HPP__
#define MSTL_CORE_NUMERIC_MATH_HPP__

 
#include "MSTL/core/exception/exception.hpp"
MSTL_BEGIN_NAMESPACE__

MSTL_BEGIN_CONSTANTS__

 
MSTL_INLINE17 constexpr decimal_t EULER = 2.718281828459045L;      
MSTL_INLINE17 constexpr decimal_t PI = 3.141592653589793L;        
MSTL_INLINE17 constexpr decimal_t PHI = 1.618033988749895L;       
MSTL_INLINE17 constexpr decimal_t SEMI_CIRCLE = 180.0;           
MSTL_INLINE17 constexpr decimal_t CIRCLE = 360.0;                
MSTL_INLINE17 constexpr decimal_t EPSILON = 1e-15L;              
 
MSTL_INLINE17 constexpr uint32_t TAYLOR_CONVERGENCE = 10000U;    
MSTL_INLINE17 constexpr decimal_t PRECISE_TOLERANCE = TAYLOR_CONVERGENCE * EPSILON;      
MSTL_INLINE17 constexpr decimal_t LOW_PRECISE_TOLERANCE = TAYLOR_CONVERGENCE * PRECISE_TOLERANCE;  

MSTL_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
};
 
MSTL_INLINE17 constexpr uint32_t FIBONACCI_COUNT = extent_v<decltype(FIBONACCI_LIST)>;   
 // MathConstants

MSTL_END_CONSTANTS__


MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 uint64_t fibonacci(const uint32_t n) {
    if (n < _CONSTANTS FIBONACCI_COUNT) {
        return _CONSTANTS FIBONACCI_LIST[n];
    }
    return fibonacci(n - 1) + fibonacci(n - 2);
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 uint64_t leonardo(const uint32_t n) {
    return 2 * fibonacci(n + 1) - 1;
}

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

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

template <typename T>
MSTL_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>
MSTL_CONST_FUNCTION constexpr enable_if_t<is_unsigned_v<T>, T>
absolute(const T x) noexcept {
    return x;
}

template <typename T>
MSTL_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>
MSTL_CONST_FUNCTION constexpr decltype(auto) sum(First first, Rests... args) {
    return first + _MSTL sum(args...);
}

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

template <typename T>
MSTL_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>
MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 T gcd(const T& m, const T& n) noexcept {
    T x = _MSTL absolute(m), y = _MSTL absolute(n);
    constexpr T zero = T(0);
    while (y != zero) {
        T t = x % y;
        x = y;
        y = t;
    }
    return x;
}

template <typename T>
MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 T lcm(const T& m, const T& n) noexcept {
    return m * n / _MSTL gcd(m, n);
}
 

template <typename T>
MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 T float_mod(const T x, const T y) {
    static_assert(is_arithmetic_v<T>, "arithmetic required");
    if (y == 0) throw_exception(math_exception("zero can not be dividend."));
    const T result = x - static_cast<make_integer_t<sizeof(T)>>(x / y) * y;
    return result;
}

template <typename T>
MSTL_CONST_FUNCTION MSTL_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;
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t exponential(const uint32_t n) noexcept {
    return power(_CONSTANTS EULER, n);
}

template <typename T>
MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 T logarithm_e(const T x) noexcept {
    static_assert(is_floating_point_v<T>, "floating point required");
    uint32_t N = _CONSTANTS TAYLOR_CONVERGENCE;
    const T a = (x - 1) / (x + 1);
    const T a_sqrt = a * a;
    T nk = 2 * N + 1;
    T y = 1.0 / nk;
    while (N--) {
        nk -= 2;
        y = 1.0 / nk + a_sqrt * y;
    }
    return 2.0 * a * y;
}

template <typename T>
MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 T logarithm(const T x, const uint32_t base) {
    const auto under = _MSTL logarithm_e(static_cast<float64_t>(base));
    if (under == 0) {
        throw_exception(math_exception("zero can not be dividend."));
    }
    return _MSTL logarithm_e(x) / under;
}

template <typename T>
MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 T logarithm_2(const T x) {
    return _MSTL logarithm(x, 2);
}

template <typename T>
MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 T logarithm_10(const T x) {
    return _MSTL logarithm(x, 10);
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t
square_root(const decimal_t x, const decimal_t precise = _CONSTANTS PRECISE_TOLERANCE) noexcept {
    decimal_t t = 0.0;
    decimal_t result = x;
    while (absolute(result - t) > precise) {
        t = result;
        result = 0.5 * (t + x / t);
    }
    return result;
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t
cube_root(const decimal_t x, const decimal_t precise = _CONSTANTS PRECISE_TOLERANCE) noexcept {
    decimal_t t = 0.0;
    decimal_t result = x;
    while (absolute(result - t) > precise) {
        t = result;
        result = (2 * t + x / (t * t)) / 3;
    }
    return result;
}

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 uint64_t factorial(const uint32_t n) noexcept {
    uint64_t h = 1;
    for (uint32_t i = 1; i <= n; i++) {
        h *= i;
    }
    return h;
}
 

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t floor_bit(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t times = power(10.0, bit);
    const auto int_part = x * times;
    if (x < 0 && x * times * 10 / 10.0 != int_part)
        return (int_part - 1) / times;
    return int_part / times;
}

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t ceil_bit(const decimal_t x, const uint32_t bit) noexcept {
    const decimal_t times = power(10.0, bit);
    const auto int_part = x * times;
    if (x > 0 && x * times * 10 / 10.0 != int_part)
        return (int_part + 1) / times;
    return int_part / times;
}

MSTL_CONST_FUNCTION MSTL_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);
}

MSTL_CONST_FUNCTION MSTL_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);
}
 

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t floor(const decimal_t x) noexcept {
    return (x >= 0) ? static_cast<int64_t>(x) : static_cast<int64_t>(x - 1);
}

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

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t ceil(const decimal_t x) noexcept {
    return (x >= 0) ? static_cast<int64_t>(x + 1) : static_cast<int64_t>(x);
}

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

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t round(const decimal_t x) noexcept {
    return (x >= 0) ? floor(x + 0.5) : ceil(x - 0.5);
}

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

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t truncate(const decimal_t x, const int bit) noexcept {
    return x < 0 ? ceil(x, bit) : floor(x, bit);
}

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t truncate(const decimal_t x) noexcept {
    return truncate(x, 0);
}
 

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 bool around_multiple(const decimal_t x, const decimal_t axis,
    const decimal_t toler = _CONSTANTS PRECISE_TOLERANCE) {
    if (absolute(axis) < _CONSTANTS PRECISE_TOLERANCE) {
        throw_exception(math_exception("Axis Cannot be 0"));
    }
    const decimal_t multi = _MSTL round(x / axis) * axis;
    return absolute(x - multi) < toler;
}

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 bool around_pi(const decimal_t x,
    const decimal_t toler = _CONSTANTS LOW_PRECISE_TOLERANCE) {
    return around_multiple(x, _CONSTANTS PI, toler);
}

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 bool around_zero(const decimal_t x,
    const decimal_t toler = _CONSTANTS PRECISE_TOLERANCE) {
    return absolute(x) < toler;
}
 

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t remainder(const decimal_t x, const decimal_t y) noexcept {
    return x - y * _MSTL round(x / y);
}

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

MSTL_CONST_FUNCTION MSTL_CONSTEXPR14 decimal_t float_apart(decimal_t x, int64_t* int_ptr) noexcept {
    *int_ptr = static_cast<int64_t>(x);
    x -= *int_ptr;
    return x;
}
 

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t sine(decimal_t x) noexcept {
    decimal_t sign = 1.0;
    if (x < 0) {
        sign = -1.0;
        x = -x;
    }
     constexpr decimal_t twoPi = 2 * _CONSTANTS PI;
    x = x - twoPi * floor(x / twoPi);
 
    if (x > _CONSTANTS PI) {
        x -= _CONSTANTS PI;
        sign = -sign;
    }
    if (x > _CONSTANTS PI / 2) {
        x = _CONSTANTS PI - x;
    }
 
    decimal_t i = 1;
    int32_t neg = 1;
    decimal_t sum = 0;
    decimal_t idx = x;
    decimal_t fac = 1;
    decimal_t taylor = x;
    do {
        fac = fac * (i + 1) * (i + 2);
        idx *= x * x;
        neg = -neg;
        sum = idx / fac * neg;
        taylor += sum;
        i += 2;
    } while (absolute(sum) > _CONSTANTS EPSILON);
    const auto fin = sign * taylor;
    return around_zero(fin) ? 0 : fin;
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t cosine(const decimal_t x) noexcept {
    return sine(_CONSTANTS PI / 2.0 - x);
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t tangent(const decimal_t x) {
    const decimal_t multiple = (2 * _MSTL round((2 * x - _CONSTANTS PI)
        / (2 * _CONSTANTS PI)) + 1) * (_CONSTANTS PI / 2);
    if (absolute(x - multiple) < _CONSTANTS LOW_PRECISE_TOLERANCE) {
        throw_exception(math_exception("Tangent Range Exceeded"));
    }
    return sine(x) / cosine(x);
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t cotangent(const decimal_t x) {
    return 1 / tangent(x);
}

MSTL_BEGIN_INNER__

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t __arctangent_taylor(const decimal_t x) noexcept {
    const decimal_t x_sq = x * x;
    decimal_t term = x;
    decimal_t sum = x;
    decimal_t n = 1.0;
    while (absolute(term) > _CONSTANTS EPSILON) {
        term *= -x_sq;
        n += 2.0;
        sum += term / n;
    }
    return sum;
}
 
MSTL_END_INNER__

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t arctangent(const decimal_t x) noexcept {
    if (absolute(x) > 1) {
        const decimal_t sign = x > 0 ? 1.0 : -1.0;
        return sign * (_CONSTANTS PI / 2 - _INNER __arctangent_taylor(1 / absolute(x)));
    }
    return _INNER __arctangent_taylor(x);
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t arcsine(const decimal_t x) {
    if (_MSTL absolute(x) > 1) {
        throw_exception(math_exception("Arcsine Range Exceeded"));
    }
    return arctangent(x / square_root(1 - x * x));
}

MSTL_PURE_FUNCTION MSTL_CONSTEXPR14 decimal_t arccosine(const decimal_t x) {
    if (_MSTL absolute(x) > 1) {
        throw_exception(math_exception("Arccosine Range Exceeded"));
    }
    return _CONSTANTS PI / 2.0 - arcsine(x);
}
 // MathFunctions
 
MSTL_END_NAMESPACE__
#endif // MSTL_CORE_NUMERIC_MATH_HPP__