real_analysis.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_COMMON_H
#define THEORETICA_COMMON_H
 
#include "./core_traits.h"
#include "./constants.h"
#include "./error.h"
 
 
namespace theoretica {
 
 
    TH_CONSTEXPR inline real identity(real x) {
        return x;
    }
 
 
    template<typename Type, typename = std::enable_if<is_real_type<Type>::value>>
    TH_CONSTEXPR inline Type conjugate(Type x) {
        return x;
    }
 
 
    TH_CONSTEXPR inline real square(real x) {
        return x * x;
    }
 
 
    TH_CONSTEXPR inline real cube(real x) {
        return x * x * x;
    }
 
 
    template<typename UnsignedIntType = uint64_t>
    inline UnsignedIntType isqrt(UnsignedIntType n) {
 
        // Upper bound
        UnsignedIntType upper = n + 1;
        
        // Lower bound
        UnsignedIntType lower = 0;
        
        // Carry for safe long division
        UnsignedIntType c = 0;
 
        while(lower != upper - 1) {
 
            // Compute carry for long division by 2
            c = ((lower % 2 != 0) && (upper % 2 != 0)) ? 1 : 0;
 
            // Safer division by 2 for big numbers
            const UnsignedIntType m = (lower >> 1) + (upper >> 1) + c;
 
            // Using division instead of multiplication avoids
            // overflows which would remove significant bits
            const UnsignedIntType q = n / m;
 
            if(m > q)
                upper = m;
            else if(m < q)
                lower = m;
            else
                return m;
        }
 
        return lower;
    }
 
 
    template<typename UnsignedIntType = uint64_t>
    inline UnsignedIntType icbrt(UnsignedIntType n) {
 
        // Upper bound
        UnsignedIntType upper = n + 1;
        
        // Lower bound
        UnsignedIntType lower = 0;
        
        // Carry for safe long division
        UnsignedIntType c = 0;
 
        while(lower != upper - 1) {
 
            // Compute carry for long division by 2
            c = ((lower % 2 != 0) && (upper % 2 != 0)) ? 1 : 0;
 
            // Safer division by 2 for big numbers
            const UnsignedIntType m = (lower >> 1) + (upper >> 1) + c;
 
            // Using division instead of multiplication avoids
            // overflows which would remove significant bits
            const UnsignedIntType q = (n / m) / m;
 
            if(m > q)
                upper = m;
            else if(m < q)
                lower = m;
            else
                return m;
        }
 
        return lower;
    }
 
 
    inline real sqrt(real x) {
 
        if(x < 0) {
            TH_MATH_ERROR("sqrt", x, MathError::OutOfDomain);
            return nan();
        }
 
#ifdef THEORETICA_X86
 
    #ifdef MSVC_ASM
        __asm {
            fld x
            fsqrt
        }
    #else
        asm("fsqrt" : "+t"(x));
    #endif
 
        return x;
#else
 
        if(x < 1) {
 
            if(x == 0)
                return 0;
 
            // Approximate sqrt(x) between 0 and 1
            // The root of the inverse is the inverse of the root
            // !!! Possible precision problems with smaller numbers
            return 1.0 / sqrt(1.0 / x);
        }
 
        // Approximate sqrt(x) using Newton-Raphson
        real y = x;
        unsigned int i = 0;
 
        while((square(y) - x) > OPTIMIZATION_TOL && i < OPTIMIZATION_NEWTON_ITER) {
            y = (y + x / y) / 2.0;
            i++;
        }
 
        return y;
#endif
 
    }
 
 
    inline real cbrt(real x) {
 
        if(x < 1) {
 
            if(x == 0)
                return 0;
 
            // cbrt(x) is odd
            if(x < 0)
                return -cbrt(-x);
 
            // Approximate cbrt between 0 and 1
            // The root of the inverse is the inverse of the root
            // !!! Possible precision problems with smaller numbers
            return 1.0 / cbrt(1.0 / x);
        }
 
        // Approximate cbrt(x) using Newton-Raphson
        real y = x;
        unsigned int i = 0;
 
        while((cube(y) - x) > OPTIMIZATION_TOL && i < OPTIMIZATION_NEWTON_ITER) {
            y = (y * 2.0 + x / (y * y)) / 3.0;
            i++;
        }
 
        return y;
    }
 
 
    inline real abs(real x) {
 
#ifdef THEORETICA_X86
 
    #ifdef MSVC_ASM
        __asm {
            fld x
            fabs
        }
    #else
        asm("fabs" : "+t"(x));
    #endif
 
        return x;
#else
        return x >= 0 ? x : -x;
#endif
 
    }
 
 
    inline int sgn(real x) {
        return (x > 0) ? 1 : (x < 0 ? -1 : 0);
    }
 
 
    TH_CONSTEXPR inline int floor(real x) {
 
        if(x < 0 && x > -1)
            return -1;
 
        // Compute the biggest integer number
        // that is smaller than x
        return x - (int(x) % 1);
    }
 
 
    inline real fract(real x) {
        return abs(x - floor(x));
    }
 
 
    inline real max(real x, real y) {
        
        #ifdef THEORETICA_BRANCHLESS
            return (x + y + abs(x - y)) / 2.0;
        #else
            return x > y ? x : y;
        #endif
    }
 
    template<typename T>
    inline T max(T x, T y) {
        return x > y ? x : y;
    }
 
 
    inline real min(real x, real y) {
 
        #ifdef THEORETICA_BRANCHLESS
            return (x + y - abs(x - y)) / 2.0;
        #else
            return x > y ? y : x;
        #endif
    }
 
    template<typename T>
    inline T min(T x, T y) {
        return x > y ? y : x;
    }
 
 
    inline real clamp(real x, real a, real b) {
        return x > b ? b : (x < a ? a : x);
    }
 
 
    template<typename T>
    inline T clamp(T x, T a, T b) {
        return x > b ? b : (x < a ? a : x);
    }
 
 
    // x86 instruction wrappers
 
#ifdef THEORETICA_X86
 
    inline real fyl2x(real x, real y) {
 
        #ifdef MSVC_ASM
 
            __asm {
                fld y
                fld x
                fyl2x
            }
 
        #else
            double r;
            asm ("fyl2x" : "=t"(r) : "0"(x), "u"(y) : "st(1)");
            return r;
        #endif
    }
 
 
    inline real f2xm1(real x) {
 
        #ifdef MSVC_ASM
 
            __asm {
                fld x
                f2xm1
            }
 
        #else
            asm("f2xm1" : "+t"(x));
            return x;
        #endif
    }
 
#endif
 
    inline real log2(real x) {
 
        if(x <= 0) {
 
            if(x == 0) {
                TH_MATH_ERROR("log2", x, MathError::OutOfRange);
                return -inf();
            }
 
            TH_MATH_ERROR("log2", x, MathError::OutOfDomain);
            return nan();
        }
 
#ifdef THEORETICA_X86
 
        // Approximate the binary logarithm of x by
        // exploiting x86 Assembly instructions
        return fyl2x(x, 1.0);
#else
 
        // Domain reduction to [1, +inf)
        if (x < 1.0)
            return -log2(1.0 / x);
 
        // Compute the biggest power of 2 so that x <= 2^i
        unsigned int i = 0;
        while(x > (uint64_t(1) << i))
            i++;
 
        // Domain reduction to [1, 2]
        x /= (uint64_t(1) << i);
 
        // Use the Taylor expansion of the logarithm
        // ln(1 + x) = \sum_k^n (-1)^(k+1) x^k / k
        real log_z = 0;
        const real z = x - 1;
 
        // Exact powers of 2 don't need further computation
        if(abs(z) > MACH_EPSILON) {
 
            int sign = 1;
            real pow_z = z;
            log_z = z;
 
            for (int i = 2; i <= 24; ++i) {
                
                pow_z *= z;
                sign *= -1;
 
                log_z += sign * pow_z / i;
            }
        }
 
        return i + (log_z / LN2);
#endif
    }
 
 
    inline real log10(real x) {
 
        if(x <= 0) {
 
            if(x == 0) {
                TH_MATH_ERROR("log10", x, MathError::OutOfRange);
                return -inf();
            }
 
            TH_MATH_ERROR("log10", x, MathError::OutOfDomain);
            return nan();
        }
 
#ifdef THEORETICA_X86
 
        // Approximate the binary logarithm of x by
        // exploiting x86 Assembly instructions
        return fyl2x(x, 1.0 / LOG210);
#else
        return log2(x) / LOG210;
#endif
    }
 
 
    inline real ln(real x) {
 
        if(x <= 0) {
 
            if(x == 0) {
                TH_MATH_ERROR("ln", x, MathError::OutOfRange);
                return -inf();
            }
 
            TH_MATH_ERROR("ln", x, MathError::OutOfDomain);
            return nan();
        }
 
#ifdef THEORETICA_X86
 
        // Approximate the binary logarithm of x by
        // exploiting x86 Assembly instructions
        return fyl2x(x, 1.0 / LOG2E);
#else
        return log2(x) / LOG2E;
#endif
    }
 
 
    template<typename UnsignedIntType = uint64_t>
    inline UnsignedIntType ilog2(UnsignedIntType x) {
 
        if(x == 0) {
            TH_MATH_ERROR("ilog2", x, MathError::OutOfRange);
            return std::numeric_limits<UnsignedIntType>::max();
        }
 
        UnsignedIntType bit = 0;
 
        // Find the highest set bit
        for (int i = (sizeof(UnsignedIntType) * 8 - 1); i > 0; --i) {
            if(x & ((UnsignedIntType) 1 << i)) {
                bit = i;
                break;
            }
        }
 
        return bit;
    }
 
 
    template<typename UnsignedIntType = uint64_t>
    inline UnsignedIntType pad2(UnsignedIntType x) {
 
        UnsignedIntType bit = 0;
 
        for (int i = (sizeof(UnsignedIntType) * 8 - 1); i > 0; --i) {
            if(x & ((UnsignedIntType) 1 << i)) {
 
                // Find the highest set bit
                bit = i;
 
                // Check if x is a power of 2
                for (int j = 0; j < i; ++j)
                    if(x & ((UnsignedIntType) 1 << j))
                        return (1 << (bit + 1));
            }
        }
 
        return (1 << bit);
    }
 
 
    template<typename T = real>
    TH_CONSTEXPR inline T pow(T x, int n) {
 
        if(n > 0) {
 
            T res = x;
            T x_sqr = x * x;
            int i = 1;
 
            // Self-multiply up to biggest power of 2
            for (; i < (n / 2); i *= 2)
                res = res * res;
 
            // Multiply by x^2 for remaining even powers
            for (; i < (n - 1); i += 2)
                res = res * x_sqr;
 
            // Multiply for remaining powers
            for (; i < n; ++i)
                res = res * x;
 
            return res;
 
        } else if(n < 0) {
            return T(1.0) / pow(x, -n);
        } else {
            return 1;
        }
    }
 
 
    template<typename T = real>
    TH_CONSTEXPR inline T ipow(T x, unsigned int n, T neutral_element = T(1.0)) {
 
        if(n == 0)
            return neutral_element;
 
        T res = x;
        T x_sqr = x * x;
        unsigned int i = 1;
 
        // Self-multiply up to biggest power of 2
        for (; i <= (n / 2); i *= 2)
            res = res * res;
 
        // Multiply by x^2 for remaining even powers
        for (; i < (n - 1); i += 2)
            res = res * x_sqr;
 
        // Multiply for remaining powers
        for (; i < n; ++i)
            res = res * x;
 
        return res;
    }
 
 
    template<typename IntType = uint64_t>
    TH_CONSTEXPR inline IntType fact(unsigned int n) {
 
        IntType res = 1;
        for (int i = n; i > 1; --i)
            res *= i;
 
        return res;
    }
 
 
    template<typename T = uint64_t>
    TH_CONSTEXPR inline T falling_fact(T x, unsigned int n) {
 
        T res = 1;
        for (unsigned int i = 0; i < n; i++)
            res *= (x - i);
 
        return res;
    }
 
 
    template<typename T = uint64_t>
    TH_CONSTEXPR inline T rising_fact(T x, unsigned int n) {
 
        T res = 1;
        for (unsigned int i = 0; i < n; i++)
            res *= (x + i);
 
        return res;
    }
 
 
    template<typename IntType = unsigned long long int>
    TH_CONSTEXPR inline IntType double_fact(unsigned int n) {
 
        IntType res = 1;
 
        for (int i = n; i > 1; i -= 2)
            res *= i;
 
        return res;
    }
 
 
#ifdef THEORETICA_X86
 
    inline real exp_x86_norm(real x) {
 
        // e^x is computed as 2^(x / ln2)
        return square(f2xm1(x / (2 * LN2)) + 1);
    }
 
#endif
 
 
    inline real exp(real x) {
 
        // Domain reduction to [0, +inf]
        if(x < 0)
            return 1.0 / exp(-x);
 
        const real fract_x = fract(x);
        const int floor_x = floor(x);
 
        // Taylor series expansion
        // Compute e^floor(x) * e^fract(x)
        
        real res = 1;
        real s_n = 1;
 
        for (int i = 1; i <= CORE_TAYLOR_ORDER; ++i) {
 
            // Recurrence formula to improve
            // numerical stability and performance
            s_n *= fract_x / (i * 4);
            res += s_n;
        }
 
        // The fractional part is divided by 4 to improve convergence
        const real sqr_r = res * res;
        return pow(E, floor_x) * sqr_r * sqr_r;
    }
 
 
    inline real expm1(real x) {
 
        if(abs(x) > 0.001)
            return exp(x) - 1;
 
        real res = 0;
        real s_n = 1;
 
        for (int i = 1; i <= 4; ++i) {
            s_n *= x / i;
            res += s_n;
        }
 
        return res;
    }
 
 
    inline real powf(real x, real a) {
 
        if(a < 0)
            return 1.0 / exp(-a * ln(abs(x)) * sgn(x));
 
        // x^a = e^(a * ln(x))
        return exp(a * ln(abs(x)) * sgn(x));
    }
 
 
    inline real root(real x, int n) {
 
        if(((n % 2 == 0) && (x < 0)) || (n == 0)) {
            TH_MATH_ERROR("root", n, MathError::OutOfDomain);
            return nan();
        }
 
        if(n < 0)
            return 1.0 / root(x, -n);
 
        // Trivial cases
        if(n == 1)
            return x;
 
        if(n == 2)
            return sqrt(x);
 
        if(n == 3)
            return cbrt(x);
 
        if(x < 1) {
 
            if(x == 0)
                return 0;
 
            // Approximate root between 0 and 1
            // The root of the inverse is the inverse of the root
            // !!! Possible precision problems with smaller numbers
            return 1.0 / root(1.0 / x, n);
        }
 
        // Approximate n-th root using Newton-Raphson.
        // If fast exponentials and logarithms are available,
        // use a first calculation to speed up convergence.
#ifdef THEORETICA_X86
        real y = exp(ln(x) / n);
#else
        real y = x;
#endif
 
        unsigned int i = 0;
 
        while(i < OPTIMIZATION_NEWTON_ITER) {
 
            const real y_pow = pow(y, n - 1);
 
            if((y_pow * y - x) < OPTIMIZATION_TOL)
                break;
 
            y = (y * (n - 1) + x / y_pow) / (real) n;
            i++;
        }
 
        if(i >= OPTIMIZATION_NEWTON_ITER) {
            TH_MATH_ERROR("root", i, MathError::NoConvergence);
            return nan();
        }
 
        return y;
    }
 
 
    inline real sin(real x) {
 
#ifdef THEORETICA_X86
 
    #ifdef MSVC_ASM
        __asm {
            fld x
            fsin
        }
    #else
        asm("fsin" : "+t"(x));
    #endif
 
        return x;
 
#else
 
        // Clamp x between -2PI and 2PI
        if(abs(x) >= TAU)
            x -= floor(x / TAU) * TAU;
 
        // Domain reduction to [-PI, PI]
        if(x > PI)
            x = PI - x;
        else if(x < -PI)
            x = -PI - x;
 
        // Compute series with recurrence formula
        real res = x;
        real s = x;
        const real sqr_x = x * x;
 
        for (int i = 1; i < 16; ++i) {
            s = s * -sqr_x / (4 * i * i + 2 * i);
            res += s;
        }
 
        return res;
#endif
    }
 
 
    inline real cos(real x) {
 
#ifdef THEORETICA_X86
 
        #ifdef MSVC_ASM
        __asm {
            fld x
            fcos
        }
        #else
        asm("fcos" : "+t"(x));
        #endif
 
        return x;
 
#else
        return sin(PI2 - x);
#endif
    }
    
 
    inline real tan(real x) {
 
#if defined(THEORETICA_X86) && !defined(MSVC_ASM)
 
        real s, c;
 
        asm ("fsincos" : "=t"(c), "=u"(s) : "0"(x));
 
        if(abs(c) < MACH_EPSILON) {
            TH_MATH_ERROR("tan", c, MathError::DivByZero);
            return nan();
        }
 
        return s / c;
#else
 
        // Reflection
        if(x < 0)
            return -tan(-x);
 
        // Domain reduction to [0, PI]
        x -= floor(x / PI) * PI;
 
        // Domain reduction to [0, PI / 4]
        if(x > (PI / 4)) {
            const real t = tan(x - PI / 4.0);
            return (1 + t) / (1 - t);
        }
 
        const real s = sin(x);
        const real c = cos(x);
 
        if(abs(c) < MACH_EPSILON) {
            TH_MATH_ERROR("tan", c, MathError::DivByZero);
            return nan();
        }
 
        return s / c;
#endif
    }
 
 
    inline real cot(real x) {
 
        real s, c;
 
#if defined(THEORETICA_X86) && !defined(MSVC_ASM)
        
        asm ("fsincos" : "=t"(c), "=u"(s) : "0"(x));
#else
        s = sin(x);
        c = cos(x);
#endif
 
        if(abs(s) < MACH_EPSILON) {
            TH_MATH_ERROR("cot", s, MathError::DivByZero);
            return nan();
        }
 
        return c / s;
    }
 
 
    inline real atan(real x) {
 
        // Odd symmetry: atan(-x) = -atan(x)
        const int sign = (x > 0) ? 1 : -1;
        x = abs(x);
 
        // Reduce to [0, 1]
        if (x > 1.0)
            return sign * (PI2 - atan(1.0 / x));
 
        // Around x = 0, atan(x) = x
        if (x < MACH_EPSILON)
            return sign * x;
 
        // Domain reduction using the half-argument formula:
        // atan(x) = 2 * atan(x / (1 + sqrt(1 + x^2)))
        constexpr real TAN_PI12 = 0.26794919243112270647;
        constexpr int MAX_ITER = 5;
        int red = 0;
 
        while (x > TAN_PI12 && red < MAX_ITER) {
 
            // Half-argument formula
            x = x / (1.0 + sqrt(1.0 + x * x));
            red++;
        }
 
        // Adaptive Taylor series with atan(x) = sum_n (-1)^n x^(2n+1) / (2n+1)
        const real x2 = x * x;
        real power = x;
        real sum = x;
 
        // Use Taylor series up to 15 terms or until convergence
        for (int n = 1; n <= 15; ++n) {
 
            // Update power
            power *= -x2;
            
            // Current term with alternating sign built into power
            const real term = power / (2 * n + 1);
            sum += term;
            
            // Convergence criterion
            if (abs(term) < MACH_EPSILON * abs(sum))
                break;
        }
 
        // Undo argument reductions
        for (int i = 0; i < red; ++i)
            sum *= 2.0;
 
        return sign * sum;
    }
 
 
    inline real asin(real x) {
 
        if(abs(x) > 1) {
            TH_MATH_ERROR("asin", x, MathError::OutOfDomain);
            return nan();
        }
 
        return atan(x / sqrt(1 - x * x));
    }
 
 
    inline real acos(real x) {
 
        if(abs(x) > 1) {
            TH_MATH_ERROR("acos", x, MathError::OutOfDomain);
            return nan();
        }
 
        if(x < 0)
            return atan(sqrt(1 - x * x) / x) + PI;
        else
            return atan(sqrt(1 - x * x) / x);
    }
 
 
    inline real atan2(real y, real x) {
 
        if(abs(x) < MACH_EPSILON) {
 
            if(abs(y) < MACH_EPSILON) {
                TH_MATH_ERROR("atan2", y, MathError::OutOfDomain);
                return nan();
            }
 
            return y >= 0 ? PI2 : -PI2;
        }
 
        if(x > 0)
            return atan(y / x);
        else
            return y >= 0 ? atan(y / x) + PI : atan(y / x) - PI;
    }
 
 
    inline real sinh(real x) {
        real exp_x = exp(x);
        return (exp_x - 1.0 / exp_x) / 2.0;
    }
 
 
    inline real cosh(real x) {
        real exp_x = exp(x);
        return (exp_x + 1.0 / exp_x) / 2.0;
    }
 
 
    inline real tanh(real x) {
        real exp_x = exp(x);
        return (exp_x - 1.0 / exp_x) / (exp_x + 1.0 / exp_x);
    }
 
 
    inline real coth(real x) {
        real exp_x = exp(x);
        return (exp_x + 1.0 / exp_x) / (exp_x - 1.0 / exp_x);
    }
 
 
    inline real asinh(real x) {
        return ln(x + sqrt(square(x) + 1));
    }
 
 
    inline real acosh(real x) {
 
        if(x < 1) {
            TH_MATH_ERROR("acosh", x, MathError::OutOfDomain);
            return nan();
        }
 
        return ln(x + sqrt(square(x) - 1));
    }
 
 
    inline real atanh(real x) {
 
        if(x < -1 || x > 1) {
            TH_MATH_ERROR("atanh", x, MathError::OutOfDomain);
            return nan();
        }
 
        return 0.5 * ln((x + 1) / (1 - x));
    }
 
 
    inline real sigmoid(real x) {
        return 1.0 / (1.0 + 1.0 / exp(x));
    }
 
 
    inline real sinc(real x) {
 
        if(abs(x) <= MACH_EPSILON)
            return 1;
 
        return sin(PI * x) / (PI * x);
    }
 
 
    inline real heaviside(real x) {
 
        if(abs(x) < MACH_EPSILON)
            return 0.5;
 
        return x > 0 ? 1 : 0;
    }
 
 
    template<typename IntType = unsigned long long int>
    TH_CONSTEXPR inline IntType binomial_coeff(unsigned int n, unsigned int m) {
 
        if(n < m) {
            TH_MATH_ERROR("binomial_coeff", n, MathError::ImpossibleOperation);
            return 0;
        }
 
        // TO-DO Check out of range
 
        IntType res = 1;
 
        for (unsigned int i = n; i > m; --i)
            res *= i;
 
        return res / fact(n - m);
    }
 
 
    template <
        typename Integer1, typename Integer2,
        typename = std::enable_if_t<std::is_integral<Integer1>::value>,
        typename = std::enable_if_t<std::is_integral<Integer2>::value>
    >
    Integer1 gcd(Integer1 a, Integer2 b) {
 
        if (b == 0)
            return a > 0 ? a : -a;
 
        return gcd(b, a % b);
    }
 
 
    TH_CONSTEXPR inline real radians(real degrees) {
        return degrees * DEG2RAD;
    }
 
 
    TH_CONSTEXPR inline real degrees(real radians) {
        return radians * RAD2DEG;
    }
 
 
    template<typename T>
    TH_CONSTEXPR inline T kronecker_delta(T i, T j) {
        return i == j ? 1 : 0;
    }
 
 
    template<typename IntType = unsigned long long int>
    TH_CONSTEXPR inline IntType catalan(unsigned int n) {
        return binomial_coeff(2 * n, n) / (n + 1);
    }
 
 
    template <
        typename Type = real,
        disable_vector<Type> = true
    >
    TH_CONSTEXPR inline Type make_error() {
        return Type(nan());
    }
 
}
 
#endif