arch/AVX/MathFunctions.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 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_MATH_FUNCTIONS_AVX_H
0011 #define EIGEN_MATH_FUNCTIONS_AVX_H
0012
0013 /* The sin and cos functions of this file are loosely derived from
0014  * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
0015  */
0016
0017 namespace Eigen {
0018
0019 namespace internal {
0020
0021 template <>
0022 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0023 psin<Packet8f>(const Packet8f& _x) {
0024   return psin_float(_x);
0025 }
0026
0027 template <>
0028 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0029 pcos<Packet8f>(const Packet8f& _x) {
0030   return pcos_float(_x);
0031 }
0032
0033 template <>
0034 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0035 plog<Packet8f>(const Packet8f& _x) {
0036   return plog_float(_x);
0037 }
0038
0039 template <>
0040 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d
0041 plog<Packet4d>(const Packet4d& _x) {
0042   return plog_double(_x);
0043 }
0044
0045 template <>
0046 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0047 plog2<Packet8f>(const Packet8f& _x) {
0048   return plog2_float(_x);
0049 }
0050
0051 template <>
0052 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d
0053 plog2<Packet4d>(const Packet4d& _x) {
0054   return plog2_double(_x);
0055 }
0056
0057 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0058 Packet8f plog1p<Packet8f>(const Packet8f& _x) {
0059   return generic_plog1p(_x);
0060 }
0061
0062 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0063 Packet8f pexpm1<Packet8f>(const Packet8f& _x) {
0064   return generic_expm1(_x);
0065 }
0066
0067 // Exponential function. Works by writing "x = m*log(2) + r" where
0068 // "m = floor(x/log(2)+1/2)" and "r" is the remainder. The result is then
0069 // "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1).
0070 template <>
0071 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0072 pexp<Packet8f>(const Packet8f& _x) {
0073   return pexp_float(_x);
0074 }
0075
0076 // Hyperbolic Tangent function.
0077 template <>
0078 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
0079 ptanh<Packet8f>(const Packet8f& _x) {
0080   return internal::generic_fast_tanh_float(_x);
0081 }
0082
0083 // Exponential function for doubles.
0084 template <>
0085 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d
0086 pexp<Packet4d>(const Packet4d& _x) {
0087   return pexp_double(_x);
0088 }
0089
0090 // Functions for sqrt.
0091 // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step
0092 // of Newton's method, at a cost of 1-2 bits of precision as opposed to the
0093 // exact solution. It does not handle +inf, or denormalized numbers correctly.
0094 // The main advantage of this approach is not just speed, but also the fact that
0095 // it can be inlined and pipelined with other computations, further reducing its
0096 // effective latency. This is similar to Quake3's fast inverse square root.
0097 // For detail see here: http://www.beyond3d.com/content/articles/8/
0098 #if EIGEN_FAST_MATH
0099 template <>
0100 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0101 Packet8f psqrt<Packet8f>(const Packet8f& _x) {
0102   Packet8f minus_half_x = pmul(_x, pset1<Packet8f>(-0.5f));
0103   Packet8f denormal_mask = pandnot(
0104       pcmp_lt(_x, pset1<Packet8f>((std::numeric_limits<float>::min)())),
0105       pcmp_lt(_x, pzero(_x)));
0106
0107   // Compute approximate reciprocal sqrt.
0108   Packet8f x = _mm256_rsqrt_ps(_x);
0109   // Do a single step of Newton's iteration.
0110   x = pmul(x, pmadd(minus_half_x, pmul(x,x), pset1<Packet8f>(1.5f)));
0111   // Flush results for denormals to zero.
0112   return pandnot(pmul(_x,x), denormal_mask);
0113 }
0114
0115 #else
0116
0117 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0118 Packet8f psqrt<Packet8f>(const Packet8f& _x) {
0119   return _mm256_sqrt_ps(_x);
0120 }
0121
0122 #endif
0123
0124 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0125 Packet4d psqrt<Packet4d>(const Packet4d& _x) {
0126   return _mm256_sqrt_pd(_x);
0127 }
0128
0129 #if EIGEN_FAST_MATH
0130 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0131 Packet8f prsqrt<Packet8f>(const Packet8f& _x) {
0132   _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inf, 0x7f800000);
0133   _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f);
0134   _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f);
0135   _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000);
0136
0137   Packet8f neg_half = pmul(_x, p8f_minus_half);
0138
0139   // select only the inverse sqrt of positive normal inputs (denormals are
0140   // flushed to zero and cause infs as well).
0141   Packet8f lt_min_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_LT_OQ);
0142   Packet8f inf_mask =  _mm256_cmp_ps(_x, p8f_inf, _CMP_EQ_OQ);
0143   Packet8f not_normal_finite_mask = _mm256_or_ps(lt_min_mask, inf_mask);
0144
0145   // Compute an approximate result using the rsqrt intrinsic.
0146   Packet8f y_approx = _mm256_rsqrt_ps(_x);
0147
0148   // Do a single step of Newton-Raphson iteration to improve the approximation.
0149   // This uses the formula y_{n+1} = y_n * (1.5 - y_n * (0.5 * x) * y_n).
0150   // It is essential to evaluate the inner term like this because forming
0151   // y_n^2 may over- or underflow.
0152   Packet8f y_newton = pmul(y_approx, pmadd(y_approx, pmul(neg_half, y_approx), p8f_one_point_five));
0153
0154   // Select the result of the Newton-Raphson step for positive normal arguments.
0155   // For other arguments, choose the output of the intrinsic. This will
0156   // return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(x) = +inf if
0157   // x is zero or a positive denormalized float (equivalent to flushing positive
0158   // denormalized inputs to zero).
0159   return pselect<Packet8f>(not_normal_finite_mask, y_approx, y_newton);
0160 }
0161
0162 #else
0163 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0164 Packet8f prsqrt<Packet8f>(const Packet8f& _x) {
0165   _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f);
0166   return _mm256_div_ps(p8f_one, _mm256_sqrt_ps(_x));
0167 }
0168 #endif
0169
0170 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0171 Packet4d prsqrt<Packet4d>(const Packet4d& _x) {
0172   _EIGEN_DECLARE_CONST_Packet4d(one, 1.0);
0173   return _mm256_div_pd(p4d_one, _mm256_sqrt_pd(_x));
0174 }
0175
0176 F16_PACKET_FUNCTION(Packet8f, Packet8h, psin)
0177 F16_PACKET_FUNCTION(Packet8f, Packet8h, pcos)
0178 F16_PACKET_FUNCTION(Packet8f, Packet8h, plog)
0179 F16_PACKET_FUNCTION(Packet8f, Packet8h, plog2)
0180 F16_PACKET_FUNCTION(Packet8f, Packet8h, plog1p)
0181 F16_PACKET_FUNCTION(Packet8f, Packet8h, pexpm1)
0182 F16_PACKET_FUNCTION(Packet8f, Packet8h, pexp)
0183 F16_PACKET_FUNCTION(Packet8f, Packet8h, ptanh)
0184 F16_PACKET_FUNCTION(Packet8f, Packet8h, psqrt)
0185 F16_PACKET_FUNCTION(Packet8f, Packet8h, prsqrt)
0186
0187 template <>
0188 EIGEN_STRONG_INLINE Packet8h pfrexp(const Packet8h& a, Packet8h& exponent) {
0189   Packet8f fexponent;
0190   const Packet8h out = float2half(pfrexp<Packet8f>(half2float(a), fexponent));
0191   exponent = float2half(fexponent);
0192   return out;
0193 }
0194
0195 template <>
0196 EIGEN_STRONG_INLINE Packet8h pldexp(const Packet8h& a, const Packet8h& exponent) {
0197   return float2half(pldexp<Packet8f>(half2float(a), half2float(exponent)));
0198 }
0199
0200 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, psin)
0201 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pcos)
0202 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, plog)
0203 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, plog2)
0204 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, plog1p)
0205 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pexpm1)
0206 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pexp)
0207 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, ptanh)
0208 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, psqrt)
0209 BF16_PACKET_FUNCTION(Packet8f, Packet8bf, prsqrt)
0210
0211 template <>
0212 EIGEN_STRONG_INLINE Packet8bf pfrexp(const Packet8bf& a, Packet8bf& exponent) {
0213   Packet8f fexponent;
0214   const Packet8bf out = F32ToBf16(pfrexp<Packet8f>(Bf16ToF32(a), fexponent));
0215   exponent = F32ToBf16(fexponent);
0216   return out;
0217 }
0218
0219 template <>
0220 EIGEN_STRONG_INLINE Packet8bf pldexp(const Packet8bf& a, const Packet8bf& exponent) {
0221   return F32ToBf16(pldexp<Packet8f>(Bf16ToF32(a), Bf16ToF32(exponent)));
0222 }
0223
0224 }  // end namespace internal
0225
0226 }  // end namespace Eigen
0227
0228 #endif  // EIGEN_MATH_FUNCTIONS_AVX_H