arch/AVX512/MathFunctions.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2016 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 THIRD_PARTY_EIGEN3_EIGEN_SRC_CORE_ARCH_AVX512_MATHFUNCTIONS_H_
0011 #define THIRD_PARTY_EIGEN3_EIGEN_SRC_CORE_ARCH_AVX512_MATHFUNCTIONS_H_
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 // Disable the code for older versions of gcc that don't support many of the required avx512 instrinsics.
0018 #if EIGEN_GNUC_AT_LEAST(5, 3) || EIGEN_COMP_CLANG  || EIGEN_COMP_MSVC >= 1923
0019
0020 #define _EIGEN_DECLARE_CONST_Packet16f(NAME, X) \
0021   const Packet16f p16f_##NAME = pset1<Packet16f>(X)
0022
0023 #define _EIGEN_DECLARE_CONST_Packet16f_FROM_INT(NAME, X) \
0024   const Packet16f p16f_##NAME =  preinterpret<Packet16f,Packet16i>(pset1<Packet16i>(X))
0025
0026 #define _EIGEN_DECLARE_CONST_Packet8d(NAME, X) \
0027   const Packet8d p8d_##NAME = pset1<Packet8d>(X)
0028
0029 #define _EIGEN_DECLARE_CONST_Packet8d_FROM_INT64(NAME, X) \
0030   const Packet8d p8d_##NAME = _mm512_castsi512_pd(_mm512_set1_epi64(X))
0031
0032 #define _EIGEN_DECLARE_CONST_Packet16bf(NAME, X) \
0033   const Packet16bf p16bf_##NAME = pset1<Packet16bf>(X)
0034
0035 #define _EIGEN_DECLARE_CONST_Packet16bf_FROM_INT(NAME, X) \
0036   const Packet16bf p16bf_##NAME =  preinterpret<Packet16bf,Packet16i>(pset1<Packet16i>(X))
0037
0038 template <>
0039 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0040 plog<Packet16f>(const Packet16f& _x) {
0041   return plog_float(_x);
0042 }
0043
0044 template <>
0045 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8d
0046 plog<Packet8d>(const Packet8d& _x) {
0047   return plog_double(_x);
0048 }
0049
0050 F16_PACKET_FUNCTION(Packet16f, Packet16h, plog)
0051 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, plog)
0052
0053 template <>
0054 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0055 plog2<Packet16f>(const Packet16f& _x) {
0056   return plog2_float(_x);
0057 }
0058
0059 template <>
0060 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8d
0061 plog2<Packet8d>(const Packet8d& _x) {
0062   return plog2_double(_x);
0063 }
0064
0065 F16_PACKET_FUNCTION(Packet16f, Packet16h, plog2)
0066 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, plog2)
0067
0068 // Exponential function. Works by writing "x = m*log(2) + r" where
0069 // "m = floor(x/log(2)+1/2)" and "r" is the remainder. The result is then
0070 // "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1).
0071 template <>
0072 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0073 pexp<Packet16f>(const Packet16f& _x) {
0074   _EIGEN_DECLARE_CONST_Packet16f(1, 1.0f);
0075   _EIGEN_DECLARE_CONST_Packet16f(half, 0.5f);
0076   _EIGEN_DECLARE_CONST_Packet16f(127, 127.0f);
0077
0078   _EIGEN_DECLARE_CONST_Packet16f(exp_hi, 88.3762626647950f);
0079   _EIGEN_DECLARE_CONST_Packet16f(exp_lo, -88.3762626647949f);
0080
0081   _EIGEN_DECLARE_CONST_Packet16f(cephes_LOG2EF, 1.44269504088896341f);
0082
0083   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p0, 1.9875691500E-4f);
0084   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p1, 1.3981999507E-3f);
0085   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p2, 8.3334519073E-3f);
0086   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p3, 4.1665795894E-2f);
0087   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p4, 1.6666665459E-1f);
0088   _EIGEN_DECLARE_CONST_Packet16f(cephes_exp_p5, 5.0000001201E-1f);
0089
0090   // Clamp x.
0091   Packet16f x = pmax(pmin(_x, p16f_exp_hi), p16f_exp_lo);
0092
0093   // Express exp(x) as exp(m*ln(2) + r), start by extracting
0094   // m = floor(x/ln(2) + 0.5).
0095   Packet16f m = _mm512_floor_ps(pmadd(x, p16f_cephes_LOG2EF, p16f_half));
0096
0097   // Get r = x - m*ln(2). Note that we can do this without losing more than one
0098   // ulp precision due to the FMA instruction.
0099   _EIGEN_DECLARE_CONST_Packet16f(nln2, -0.6931471805599453f);
0100   Packet16f r = _mm512_fmadd_ps(m, p16f_nln2, x);
0101   Packet16f r2 = pmul(r, r);
0102   Packet16f r3 = pmul(r2, r);
0103
0104   // Evaluate the polynomial approximant,improved by instruction-level parallelism.
0105   Packet16f y, y1, y2;
0106   y  = pmadd(p16f_cephes_exp_p0, r, p16f_cephes_exp_p1);
0107   y1 = pmadd(p16f_cephes_exp_p3, r, p16f_cephes_exp_p4);
0108   y2 = padd(r, p16f_1);
0109   y  = pmadd(y, r, p16f_cephes_exp_p2);
0110   y1 = pmadd(y1, r, p16f_cephes_exp_p5);
0111   y  = pmadd(y, r3, y1);
0112   y  = pmadd(y, r2, y2);
0113
0114   // Build emm0 = 2^m.
0115   Packet16i emm0 = _mm512_cvttps_epi32(padd(m, p16f_127));
0116   emm0 = _mm512_slli_epi32(emm0, 23);
0117
0118   // Return 2^m * exp(r).
0119   return pmax(pmul(y, _mm512_castsi512_ps(emm0)), _x);
0120 }
0121
0122 template <>
0123 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8d
0124 pexp<Packet8d>(const Packet8d& _x) {
0125   return pexp_double(_x);
0126 }
0127
0128 F16_PACKET_FUNCTION(Packet16f, Packet16h, pexp)
0129 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pexp)
0130
0131 template <>
0132 EIGEN_STRONG_INLINE Packet16h pfrexp(const Packet16h& a, Packet16h& exponent) {
0133   Packet16f fexponent;
0134   const Packet16h out = float2half(pfrexp<Packet16f>(half2float(a), fexponent));
0135   exponent = float2half(fexponent);
0136   return out;
0137 }
0138
0139 template <>
0140 EIGEN_STRONG_INLINE Packet16h pldexp(const Packet16h& a, const Packet16h& exponent) {
0141   return float2half(pldexp<Packet16f>(half2float(a), half2float(exponent)));
0142 }
0143
0144 template <>
0145 EIGEN_STRONG_INLINE Packet16bf pfrexp(const Packet16bf& a, Packet16bf& exponent) {
0146   Packet16f fexponent;
0147   const Packet16bf out = F32ToBf16(pfrexp<Packet16f>(Bf16ToF32(a), fexponent));
0148   exponent = F32ToBf16(fexponent);
0149   return out;
0150 }
0151
0152 template <>
0153 EIGEN_STRONG_INLINE Packet16bf pldexp(const Packet16bf& a, const Packet16bf& exponent) {
0154   return F32ToBf16(pldexp<Packet16f>(Bf16ToF32(a), Bf16ToF32(exponent)));
0155 }
0156
0157 // Functions for sqrt.
0158 // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step
0159 // of Newton's method, at a cost of 1-2 bits of precision as opposed to the
0160 // exact solution. The main advantage of this approach is not just speed, but
0161 // also the fact that it can be inlined and pipelined with other computations,
0162 // further reducing its effective latency.
0163 #if EIGEN_FAST_MATH
0164 template <>
0165 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0166 psqrt<Packet16f>(const Packet16f& _x) {
0167   Packet16f neg_half = pmul(_x, pset1<Packet16f>(-.5f));
0168   __mmask16 denormal_mask = _mm512_kand(
0169       _mm512_cmp_ps_mask(_x, pset1<Packet16f>((std::numeric_limits<float>::min)()),
0170                         _CMP_LT_OQ),
0171       _mm512_cmp_ps_mask(_x, _mm512_setzero_ps(), _CMP_GE_OQ));
0172
0173   Packet16f x = _mm512_rsqrt14_ps(_x);
0174
0175   // Do a single step of Newton's iteration.
0176   x = pmul(x, pmadd(neg_half, pmul(x, x), pset1<Packet16f>(1.5f)));
0177
0178   // Flush results for denormals to zero.
0179   return _mm512_mask_blend_ps(denormal_mask, pmul(_x,x), _mm512_setzero_ps());
0180 }
0181
0182 template <>
0183 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8d
0184 psqrt<Packet8d>(const Packet8d& _x) {
0185   Packet8d neg_half = pmul(_x, pset1<Packet8d>(-.5));
0186   __mmask16 denormal_mask = _mm512_kand(
0187       _mm512_cmp_pd_mask(_x, pset1<Packet8d>((std::numeric_limits<double>::min)()),
0188                         _CMP_LT_OQ),
0189       _mm512_cmp_pd_mask(_x, _mm512_setzero_pd(), _CMP_GE_OQ));
0190
0191   Packet8d x = _mm512_rsqrt14_pd(_x);
0192
0193   // Do a single step of Newton's iteration.
0194   x = pmul(x, pmadd(neg_half, pmul(x, x), pset1<Packet8d>(1.5)));
0195
0196   // Do a second step of Newton's iteration.
0197   x = pmul(x, pmadd(neg_half, pmul(x, x), pset1<Packet8d>(1.5)));
0198
0199   return _mm512_mask_blend_pd(denormal_mask, pmul(_x,x), _mm512_setzero_pd());
0200 }
0201 #else
0202 template <>
0203 EIGEN_STRONG_INLINE Packet16f psqrt<Packet16f>(const Packet16f& x) {
0204   return _mm512_sqrt_ps(x);
0205 }
0206
0207 template <>
0208 EIGEN_STRONG_INLINE Packet8d psqrt<Packet8d>(const Packet8d& x) {
0209   return _mm512_sqrt_pd(x);
0210 }
0211 #endif
0212
0213 F16_PACKET_FUNCTION(Packet16f, Packet16h, psqrt)
0214 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, psqrt)
0215
0216 // prsqrt for float.
0217 #if defined(EIGEN_VECTORIZE_AVX512ER)
0218
0219 template <>
0220 EIGEN_STRONG_INLINE Packet16f prsqrt<Packet16f>(const Packet16f& x) {
0221   return _mm512_rsqrt28_ps(x);
0222 }
0223 #elif EIGEN_FAST_MATH
0224
0225 template <>
0226 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0227 prsqrt<Packet16f>(const Packet16f& _x) {
0228   _EIGEN_DECLARE_CONST_Packet16f_FROM_INT(inf, 0x7f800000);
0229   _EIGEN_DECLARE_CONST_Packet16f(one_point_five, 1.5f);
0230   _EIGEN_DECLARE_CONST_Packet16f(minus_half, -0.5f);
0231
0232   Packet16f neg_half = pmul(_x, p16f_minus_half);
0233
0234   // Identity infinite, negative and denormal arguments.
0235   __mmask16 inf_mask = _mm512_cmp_ps_mask(_x, p16f_inf, _CMP_EQ_OQ);
0236   __mmask16 not_pos_mask = _mm512_cmp_ps_mask(_x, _mm512_setzero_ps(), _CMP_LE_OQ);
0237   __mmask16 not_finite_pos_mask = not_pos_mask | inf_mask;
0238
0239   // Compute an approximate result using the rsqrt intrinsic, forcing +inf
0240   // for denormals for consistency with AVX and SSE implementations.
0241   Packet16f y_approx = _mm512_rsqrt14_ps(_x);
0242
0243   // Do a single step of Newton-Raphson iteration to improve the approximation.
0244   // This uses the formula y_{n+1} = y_n * (1.5 - y_n * (0.5 * x) * y_n).
0245   // It is essential to evaluate the inner term like this because forming
0246   // y_n^2 may over- or underflow.
0247   Packet16f y_newton = pmul(y_approx, pmadd(y_approx, pmul(neg_half, y_approx), p16f_one_point_five));
0248
0249   // Select the result of the Newton-Raphson step for positive finite arguments.
0250   // For other arguments, choose the output of the intrinsic. This will
0251   // return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(0) = +inf.
0252   return _mm512_mask_blend_ps(not_finite_pos_mask, y_newton, y_approx);
0253 }
0254 #else
0255
0256 template <>
0257 EIGEN_STRONG_INLINE Packet16f prsqrt<Packet16f>(const Packet16f& x) {
0258   _EIGEN_DECLARE_CONST_Packet16f(one, 1.0f);
0259   return _mm512_div_ps(p16f_one, _mm512_sqrt_ps(x));
0260 }
0261 #endif
0262
0263 F16_PACKET_FUNCTION(Packet16f, Packet16h, prsqrt)
0264 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, prsqrt)
0265
0266 // prsqrt for double.
0267 #if EIGEN_FAST_MATH
0268 template <>
0269 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8d
0270 prsqrt<Packet8d>(const Packet8d& _x) {
0271   _EIGEN_DECLARE_CONST_Packet8d(one_point_five, 1.5);
0272   _EIGEN_DECLARE_CONST_Packet8d(minus_half, -0.5);
0273   _EIGEN_DECLARE_CONST_Packet8d_FROM_INT64(inf, 0x7ff0000000000000LL);
0274
0275   Packet8d neg_half = pmul(_x, p8d_minus_half);
0276
0277   // Identity infinite, negative and denormal arguments.
0278   __mmask8 inf_mask = _mm512_cmp_pd_mask(_x, p8d_inf, _CMP_EQ_OQ);
0279   __mmask8 not_pos_mask = _mm512_cmp_pd_mask(_x, _mm512_setzero_pd(), _CMP_LE_OQ);
0280   __mmask8 not_finite_pos_mask = not_pos_mask | inf_mask;
0281
0282   // Compute an approximate result using the rsqrt intrinsic, forcing +inf
0283   // for denormals for consistency with AVX and SSE implementations.
0284 #if defined(EIGEN_VECTORIZE_AVX512ER)
0285   Packet8d y_approx = _mm512_rsqrt28_pd(_x);
0286 #else
0287   Packet8d y_approx = _mm512_rsqrt14_pd(_x);
0288 #endif
0289   // Do one or two steps of Newton-Raphson's to improve the approximation, depending on the
0290   // starting accuracy (either 2^-14 or 2^-28, depending on whether AVX512ER is available).
0291   // The Newton-Raphson algorithm has quadratic convergence and roughly doubles the number
0292   // of correct digits for each step.
0293   // This uses the formula y_{n+1} = y_n * (1.5 - y_n * (0.5 * x) * y_n).
0294   // It is essential to evaluate the inner term like this because forming
0295   // y_n^2 may over- or underflow.
0296   Packet8d y_newton = pmul(y_approx, pmadd(neg_half, pmul(y_approx, y_approx), p8d_one_point_five));
0297 #if !defined(EIGEN_VECTORIZE_AVX512ER)
0298   y_newton = pmul(y_newton, pmadd(y_newton, pmul(neg_half, y_newton), p8d_one_point_five));
0299 #endif
0300   // Select the result of the Newton-Raphson step for positive finite arguments.
0301   // For other arguments, choose the output of the intrinsic. This will
0302   // return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(0) = +inf.
0303   return _mm512_mask_blend_pd(not_finite_pos_mask, y_newton, y_approx);
0304 }
0305 #else
0306 template <>
0307 EIGEN_STRONG_INLINE Packet8d prsqrt<Packet8d>(const Packet8d& x) {
0308   _EIGEN_DECLARE_CONST_Packet8d(one, 1.0f);
0309   return _mm512_div_pd(p8d_one, _mm512_sqrt_pd(x));
0310 }
0311 #endif
0312
0313 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0314 Packet16f plog1p<Packet16f>(const Packet16f& _x) {
0315   return generic_plog1p(_x);
0316 }
0317
0318 F16_PACKET_FUNCTION(Packet16f, Packet16h, plog1p)
0319 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, plog1p)
0320
0321 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
0322 Packet16f pexpm1<Packet16f>(const Packet16f& _x) {
0323   return generic_expm1(_x);
0324 }
0325
0326 F16_PACKET_FUNCTION(Packet16f, Packet16h, pexpm1)
0327 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pexpm1)
0328
0329 #endif
0330
0331
0332 template <>
0333 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0334 psin<Packet16f>(const Packet16f& _x) {
0335   return psin_float(_x);
0336 }
0337
0338 template <>
0339 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0340 pcos<Packet16f>(const Packet16f& _x) {
0341   return pcos_float(_x);
0342 }
0343
0344 template <>
0345 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
0346 ptanh<Packet16f>(const Packet16f& _x) {
0347   return internal::generic_fast_tanh_float(_x);
0348 }
0349
0350 F16_PACKET_FUNCTION(Packet16f, Packet16h, psin)
0351 F16_PACKET_FUNCTION(Packet16f, Packet16h, pcos)
0352 F16_PACKET_FUNCTION(Packet16f, Packet16h, ptanh)
0353
0354 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, psin)
0355 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pcos)
0356 BF16_PACKET_FUNCTION(Packet16f, Packet16bf, ptanh)
0357
0358 }  // end namespace internal
0359
0360 }  // end namespace Eigen
0361
0362 #endif  // THIRD_PARTY_EIGEN3_EIGEN_SRC_CORE_ARCH_AVX512_MATHFUNCTIONS_H_