src/Core/MathFunctionsImpl.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com)
0005 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_MATHFUNCTIONSIMPL_H
0012 #define EIGEN_MATHFUNCTIONSIMPL_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 /** \internal \returns the hyperbolic tan of \a a (coeff-wise)
0019     Doesn't do anything fancy, just a 13/6-degree rational interpolant which
0020     is accurate up to a couple of ulps in the (approximate) range [-8, 8],
0021     outside of which tanh(x) = +/-1 in single precision. The input is clamped
0022     to the range [-c, c]. The value c is chosen as the smallest value where
0023     the approximation evaluates to exactly 1. In the reange [-0.0004, 0.0004]
0024     the approxmation tanh(x) ~= x is used for better accuracy as x tends to zero.
0025
0026     This implementation works on both scalars and packets.
0027 */
0028 template<typename T>
0029 T generic_fast_tanh_float(const T& a_x)
0030 {
0031   // Clamp the inputs to the range [-c, c]
0032 #ifdef EIGEN_VECTORIZE_FMA
0033   const T plus_clamp = pset1<T>(7.99881172180175781f);
0034   const T minus_clamp = pset1<T>(-7.99881172180175781f);
0035 #else
0036   const T plus_clamp = pset1<T>(7.90531110763549805f);
0037   const T minus_clamp = pset1<T>(-7.90531110763549805f);
0038 #endif
0039   const T tiny = pset1<T>(0.0004f);
0040   const T x = pmax(pmin(a_x, plus_clamp), minus_clamp);
0041   const T tiny_mask = pcmp_lt(pabs(a_x), tiny);
0042   // The monomial coefficients of the numerator polynomial (odd).
0043   const T alpha_1 = pset1<T>(4.89352455891786e-03f);
0044   const T alpha_3 = pset1<T>(6.37261928875436e-04f);
0045   const T alpha_5 = pset1<T>(1.48572235717979e-05f);
0046   const T alpha_7 = pset1<T>(5.12229709037114e-08f);
0047   const T alpha_9 = pset1<T>(-8.60467152213735e-11f);
0048   const T alpha_11 = pset1<T>(2.00018790482477e-13f);
0049   const T alpha_13 = pset1<T>(-2.76076847742355e-16f);
0050
0051   // The monomial coefficients of the denominator polynomial (even).
0052   const T beta_0 = pset1<T>(4.89352518554385e-03f);
0053   const T beta_2 = pset1<T>(2.26843463243900e-03f);
0054   const T beta_4 = pset1<T>(1.18534705686654e-04f);
0055   const T beta_6 = pset1<T>(1.19825839466702e-06f);
0056
0057   // Since the polynomials are odd/even, we need x^2.
0058   const T x2 = pmul(x, x);
0059
0060   // Evaluate the numerator polynomial p.
0061   T p = pmadd(x2, alpha_13, alpha_11);
0062   p = pmadd(x2, p, alpha_9);
0063   p = pmadd(x2, p, alpha_7);
0064   p = pmadd(x2, p, alpha_5);
0065   p = pmadd(x2, p, alpha_3);
0066   p = pmadd(x2, p, alpha_1);
0067   p = pmul(x, p);
0068
0069   // Evaluate the denominator polynomial q.
0070   T q = pmadd(x2, beta_6, beta_4);
0071   q = pmadd(x2, q, beta_2);
0072   q = pmadd(x2, q, beta_0);
0073
0074   // Divide the numerator by the denominator.
0075   return pselect(tiny_mask, x, pdiv(p, q));
0076 }
0077
0078 template<typename RealScalar>
0079 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0080 RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y)
0081 {
0082   // IEEE IEC 6059 special cases.
0083   if ((numext::isinf)(x) || (numext::isinf)(y))
0084     return NumTraits<RealScalar>::infinity();
0085   if ((numext::isnan)(x) || (numext::isnan)(y))
0086     return NumTraits<RealScalar>::quiet_NaN();
0087
0088   EIGEN_USING_STD(sqrt);
0089   RealScalar p, qp;
0090   p = numext::maxi(x,y);
0091   if(p==RealScalar(0)) return RealScalar(0);
0092   qp = numext::mini(y,x) / p;
0093   return p * sqrt(RealScalar(1) + qp*qp);
0094 }
0095
0096 template<typename Scalar>
0097 struct hypot_impl
0098 {
0099   typedef typename NumTraits<Scalar>::Real RealScalar;
0100   static EIGEN_DEVICE_FUNC
0101   inline RealScalar run(const Scalar& x, const Scalar& y)
0102   {
0103     EIGEN_USING_STD(abs);
0104     return positive_real_hypot<RealScalar>(abs(x), abs(y));
0105   }
0106 };
0107
0108 // Generic complex sqrt implementation that correctly handles corner cases
0109 // according to https://en.cppreference.com/w/cpp/numeric/complex/sqrt
0110 template<typename T>
0111 EIGEN_DEVICE_FUNC std::complex<T> complex_sqrt(const std::complex<T>& z) {
0112   // Computes the principal sqrt of the input.
0113   //
0114   // For a complex square root of the number x + i*y. We want to find real
0115   // numbers u and v such that
0116   //    (u + i*v)^2 = x + i*y  <=>
0117   //    u^2 - v^2 + i*2*u*v = x + i*v.
0118   // By equating the real and imaginary parts we get:
0119   //    u^2 - v^2 = x
0120   //    2*u*v = y.
0121   //
0122   // For x >= 0, this has the numerically stable solution
0123   //    u = sqrt(0.5 * (x + sqrt(x^2 + y^2)))
0124   //    v = y / (2 * u)
0125   // and for x < 0,
0126   //    v = sign(y) * sqrt(0.5 * (-x + sqrt(x^2 + y^2)))
0127   //    u = y / (2 * v)
0128   //
0129   // Letting w = sqrt(0.5 * (|x| + |z|)),
0130   //   if x == 0: u = w, v = sign(y) * w
0131   //   if x > 0:  u = w, v = y / (2 * w)
0132   //   if x < 0:  u = |y| / (2 * w), v = sign(y) * w
0133
0134   const T x = numext::real(z);
0135   const T y = numext::imag(z);
0136   const T zero = T(0);
0137   const T w = numext::sqrt(T(0.5) * (numext::abs(x) + numext::hypot(x, y)));
0138
0139   return
0140     (numext::isinf)(y) ? std::complex<T>(NumTraits<T>::infinity(), y)
0141       : x == zero ? std::complex<T>(w, y < zero ? -w : w)
0142       : x > zero ? std::complex<T>(w, y / (2 * w))
0143       : std::complex<T>(numext::abs(y) / (2 * w), y < zero ? -w : w );
0144 }
0145
0146 // Generic complex rsqrt implementation.
0147 template<typename T>
0148 EIGEN_DEVICE_FUNC std::complex<T> complex_rsqrt(const std::complex<T>& z) {
0149   // Computes the principal reciprocal sqrt of the input.
0150   //
0151   // For a complex reciprocal square root of the number z = x + i*y. We want to
0152   // find real numbers u and v such that
0153   //    (u + i*v)^2 = 1 / (x + i*y)  <=>
0154   //    u^2 - v^2 + i*2*u*v = x/|z|^2 - i*v/|z|^2.
0155   // By equating the real and imaginary parts we get:
0156   //    u^2 - v^2 = x/|z|^2
0157   //    2*u*v = y/|z|^2.
0158   //
0159   // For x >= 0, this has the numerically stable solution
0160   //    u = sqrt(0.5 * (x + |z|)) / |z|
0161   //    v = -y / (2 * u * |z|)
0162   // and for x < 0,
0163   //    v = -sign(y) * sqrt(0.5 * (-x + |z|)) / |z|
0164   //    u = -y / (2 * v * |z|)
0165   //
0166   // Letting w = sqrt(0.5 * (|x| + |z|)),
0167   //   if x == 0: u = w / |z|, v = -sign(y) * w / |z|
0168   //   if x > 0:  u = w / |z|, v = -y / (2 * w * |z|)
0169   //   if x < 0:  u = |y| / (2 * w * |z|), v = -sign(y) * w / |z|
0170
0171   const T x = numext::real(z);
0172   const T y = numext::imag(z);
0173   const T zero = T(0);
0174
0175   const T abs_z = numext::hypot(x, y);
0176   const T w = numext::sqrt(T(0.5) * (numext::abs(x) + abs_z));
0177   const T woz = w / abs_z;
0178   // Corner cases consistent with 1/sqrt(z) on gcc/clang.
0179   return
0180     abs_z == zero ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN())
0181       : ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero)
0182       : x == zero ? std::complex<T>(woz, y < zero ? woz : -woz)
0183       : x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z))
0184       : std::complex<T>(numext::abs(y) / (2 * w * abs_z), y < zero ? woz : -woz );
0185 }
0186
0187 template<typename T>
0188 EIGEN_DEVICE_FUNC std::complex<T> complex_log(const std::complex<T>& z) {
0189   // Computes complex log.
0190   T a = numext::abs(z);
0191   EIGEN_USING_STD(atan2);
0192   T b = atan2(z.imag(), z.real());
0193   return std::complex<T>(numext::log(a), b);
0194 }
0195
0196 } // end namespace internal
0197
0198 } // end namespace Eigen
0199
0200 #endif // EIGEN_MATHFUNCTIONSIMPL_H