diff options
Diffstat (limited to 'src/3rdparty/eigen/Eigen/src/Core/MathFunctionsImpl.h')
-rw-r--r-- | src/3rdparty/eigen/Eigen/src/Core/MathFunctionsImpl.h | 200 |
1 files changed, 200 insertions, 0 deletions
diff --git a/src/3rdparty/eigen/Eigen/src/Core/MathFunctionsImpl.h b/src/3rdparty/eigen/Eigen/src/Core/MathFunctionsImpl.h new file mode 100644 index 000000000..4eaaaa784 --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/Core/MathFunctionsImpl.h @@ -0,0 +1,200 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com) +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_MATHFUNCTIONSIMPL_H +#define EIGEN_MATHFUNCTIONSIMPL_H + +namespace Eigen { + +namespace internal { + +/** \internal \returns the hyperbolic tan of \a a (coeff-wise) + Doesn't do anything fancy, just a 13/6-degree rational interpolant which + is accurate up to a couple of ulps in the (approximate) range [-8, 8], + outside of which tanh(x) = +/-1 in single precision. The input is clamped + to the range [-c, c]. The value c is chosen as the smallest value where + the approximation evaluates to exactly 1. In the reange [-0.0004, 0.0004] + the approxmation tanh(x) ~= x is used for better accuracy as x tends to zero. + + This implementation works on both scalars and packets. +*/ +template<typename T> +T generic_fast_tanh_float(const T& a_x) +{ + // Clamp the inputs to the range [-c, c] +#ifdef EIGEN_VECTORIZE_FMA + const T plus_clamp = pset1<T>(7.99881172180175781f); + const T minus_clamp = pset1<T>(-7.99881172180175781f); +#else + const T plus_clamp = pset1<T>(7.90531110763549805f); + const T minus_clamp = pset1<T>(-7.90531110763549805f); +#endif + const T tiny = pset1<T>(0.0004f); + const T x = pmax(pmin(a_x, plus_clamp), minus_clamp); + const T tiny_mask = pcmp_lt(pabs(a_x), tiny); + // The monomial coefficients of the numerator polynomial (odd). + const T alpha_1 = pset1<T>(4.89352455891786e-03f); + const T alpha_3 = pset1<T>(6.37261928875436e-04f); + const T alpha_5 = pset1<T>(1.48572235717979e-05f); + const T alpha_7 = pset1<T>(5.12229709037114e-08f); + const T alpha_9 = pset1<T>(-8.60467152213735e-11f); + const T alpha_11 = pset1<T>(2.00018790482477e-13f); + const T alpha_13 = pset1<T>(-2.76076847742355e-16f); + + // The monomial coefficients of the denominator polynomial (even). + const T beta_0 = pset1<T>(4.89352518554385e-03f); + const T beta_2 = pset1<T>(2.26843463243900e-03f); + const T beta_4 = pset1<T>(1.18534705686654e-04f); + const T beta_6 = pset1<T>(1.19825839466702e-06f); + + // Since the polynomials are odd/even, we need x^2. + const T x2 = pmul(x, x); + + // Evaluate the numerator polynomial p. + T p = pmadd(x2, alpha_13, alpha_11); + p = pmadd(x2, p, alpha_9); + p = pmadd(x2, p, alpha_7); + p = pmadd(x2, p, alpha_5); + p = pmadd(x2, p, alpha_3); + p = pmadd(x2, p, alpha_1); + p = pmul(x, p); + + // Evaluate the denominator polynomial q. + T q = pmadd(x2, beta_6, beta_4); + q = pmadd(x2, q, beta_2); + q = pmadd(x2, q, beta_0); + + // Divide the numerator by the denominator. + return pselect(tiny_mask, x, pdiv(p, q)); +} + +template<typename RealScalar> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y) +{ + // IEEE IEC 6059 special cases. + if ((numext::isinf)(x) || (numext::isinf)(y)) + return NumTraits<RealScalar>::infinity(); + if ((numext::isnan)(x) || (numext::isnan)(y)) + return NumTraits<RealScalar>::quiet_NaN(); + + EIGEN_USING_STD(sqrt); + RealScalar p, qp; + p = numext::maxi(x,y); + if(p==RealScalar(0)) return RealScalar(0); + qp = numext::mini(y,x) / p; + return p * sqrt(RealScalar(1) + qp*qp); +} + +template<typename Scalar> +struct hypot_impl +{ + typedef typename NumTraits<Scalar>::Real RealScalar; + static EIGEN_DEVICE_FUNC + inline RealScalar run(const Scalar& x, const Scalar& y) + { + EIGEN_USING_STD(abs); + return positive_real_hypot<RealScalar>(abs(x), abs(y)); + } +}; + +// Generic complex sqrt implementation that correctly handles corner cases +// according to https://en.cppreference.com/w/cpp/numeric/complex/sqrt +template<typename T> +EIGEN_DEVICE_FUNC std::complex<T> complex_sqrt(const std::complex<T>& z) { + // Computes the principal sqrt of the input. + // + // For a complex square root of the number x + i*y. We want to find real + // numbers u and v such that + // (u + i*v)^2 = x + i*y <=> + // u^2 - v^2 + i*2*u*v = x + i*v. + // By equating the real and imaginary parts we get: + // u^2 - v^2 = x + // 2*u*v = y. + // + // For x >= 0, this has the numerically stable solution + // u = sqrt(0.5 * (x + sqrt(x^2 + y^2))) + // v = y / (2 * u) + // and for x < 0, + // v = sign(y) * sqrt(0.5 * (-x + sqrt(x^2 + y^2))) + // u = y / (2 * v) + // + // Letting w = sqrt(0.5 * (|x| + |z|)), + // if x == 0: u = w, v = sign(y) * w + // if x > 0: u = w, v = y / (2 * w) + // if x < 0: u = |y| / (2 * w), v = sign(y) * w + + const T x = numext::real(z); + const T y = numext::imag(z); + const T zero = T(0); + const T w = numext::sqrt(T(0.5) * (numext::abs(x) + numext::hypot(x, y))); + + return + (numext::isinf)(y) ? std::complex<T>(NumTraits<T>::infinity(), y) + : x == zero ? std::complex<T>(w, y < zero ? -w : w) + : x > zero ? std::complex<T>(w, y / (2 * w)) + : std::complex<T>(numext::abs(y) / (2 * w), y < zero ? -w : w ); +} + +// Generic complex rsqrt implementation. +template<typename T> +EIGEN_DEVICE_FUNC std::complex<T> complex_rsqrt(const std::complex<T>& z) { + // Computes the principal reciprocal sqrt of the input. + // + // For a complex reciprocal square root of the number z = x + i*y. We want to + // find real numbers u and v such that + // (u + i*v)^2 = 1 / (x + i*y) <=> + // u^2 - v^2 + i*2*u*v = x/|z|^2 - i*v/|z|^2. + // By equating the real and imaginary parts we get: + // u^2 - v^2 = x/|z|^2 + // 2*u*v = y/|z|^2. + // + // For x >= 0, this has the numerically stable solution + // u = sqrt(0.5 * (x + |z|)) / |z| + // v = -y / (2 * u * |z|) + // and for x < 0, + // v = -sign(y) * sqrt(0.5 * (-x + |z|)) / |z| + // u = -y / (2 * v * |z|) + // + // Letting w = sqrt(0.5 * (|x| + |z|)), + // if x == 0: u = w / |z|, v = -sign(y) * w / |z| + // if x > 0: u = w / |z|, v = -y / (2 * w * |z|) + // if x < 0: u = |y| / (2 * w * |z|), v = -sign(y) * w / |z| + + const T x = numext::real(z); + const T y = numext::imag(z); + const T zero = T(0); + + const T abs_z = numext::hypot(x, y); + const T w = numext::sqrt(T(0.5) * (numext::abs(x) + abs_z)); + const T woz = w / abs_z; + // Corner cases consistent with 1/sqrt(z) on gcc/clang. + return + abs_z == zero ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN()) + : ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero) + : x == zero ? std::complex<T>(woz, y < zero ? woz : -woz) + : x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z)) + : std::complex<T>(numext::abs(y) / (2 * w * abs_z), y < zero ? woz : -woz ); +} + +template<typename T> +EIGEN_DEVICE_FUNC std::complex<T> complex_log(const std::complex<T>& z) { + // Computes complex log. + T a = numext::abs(z); + EIGEN_USING_STD(atan2); + T b = atan2(z.imag(), z.real()); + return std::complex<T>(numext::log(a), b); +} + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_MATHFUNCTIONSIMPL_H |