![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com) 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_MATH_FUNCTIONS_AVX_H 00011 #define EIGEN_MATH_FUNCTIONS_AVX_H 00012 00013 /* The sin, cos, exp, and log functions of this file are loosely derived from 00014 * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/ 00015 */ 00016 00017 namespace Eigen { 00018 00019 namespace internal { 00020 00021 inline Packet8i pshiftleft(Packet8i v, int n) 00022 { 00023 #ifdef EIGEN_VECTORIZE_AVX2 00024 return _mm256_slli_epi32(v, n); 00025 #else 00026 __m128i lo = _mm_slli_epi32(_mm256_extractf128_si256(v, 0), n); 00027 __m128i hi = _mm_slli_epi32(_mm256_extractf128_si256(v, 1), n); 00028 return _mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1); 00029 #endif 00030 } 00031 00032 inline Packet8f pshiftright(Packet8f v, int n) 00033 { 00034 #ifdef EIGEN_VECTORIZE_AVX2 00035 return _mm256_cvtepi32_ps(_mm256_srli_epi32(_mm256_castps_si256(v), n)); 00036 #else 00037 __m128i lo = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 0), n); 00038 __m128i hi = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 1), n); 00039 return _mm256_cvtepi32_ps(_mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1)); 00040 #endif 00041 } 00042 00043 // Sine function 00044 // Computes sin(x) by wrapping x to the interval [-Pi/4,3*Pi/4] and 00045 // evaluating interpolants in [-Pi/4,Pi/4] or [Pi/4,3*Pi/4]. The interpolants 00046 // are (anti-)symmetric and thus have only odd/even coefficients 00047 template <> 00048 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 00049 psin<Packet8f>(const Packet8f& _x) { 00050 Packet8f x = _x; 00051 00052 // Some useful values. 00053 _EIGEN_DECLARE_CONST_Packet8i(one, 1); 00054 _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f); 00055 _EIGEN_DECLARE_CONST_Packet8f(two, 2.0f); 00056 _EIGEN_DECLARE_CONST_Packet8f(one_over_four, 0.25f); 00057 _EIGEN_DECLARE_CONST_Packet8f(one_over_pi, 3.183098861837907e-01f); 00058 _EIGEN_DECLARE_CONST_Packet8f(neg_pi_first, -3.140625000000000e+00f); 00059 _EIGEN_DECLARE_CONST_Packet8f(neg_pi_second, -9.670257568359375e-04f); 00060 _EIGEN_DECLARE_CONST_Packet8f(neg_pi_third, -6.278329571784980e-07f); 00061 _EIGEN_DECLARE_CONST_Packet8f(four_over_pi, 1.273239544735163e+00f); 00062 00063 // Map x from [-Pi/4,3*Pi/4] to z in [-1,3] and subtract the shifted period. 00064 Packet8f z = pmul(x, p8f_one_over_pi); 00065 Packet8f shift = _mm256_floor_ps(padd(z, p8f_one_over_four)); 00066 x = pmadd(shift, p8f_neg_pi_first, x); 00067 x = pmadd(shift, p8f_neg_pi_second, x); 00068 x = pmadd(shift, p8f_neg_pi_third, x); 00069 z = pmul(x, p8f_four_over_pi); 00070 00071 // Make a mask for the entries that need flipping, i.e. wherever the shift 00072 // is odd. 00073 Packet8i shift_ints = _mm256_cvtps_epi32(shift); 00074 Packet8i shift_isodd = _mm256_castps_si256(_mm256_and_ps(_mm256_castsi256_ps(shift_ints), _mm256_castsi256_ps(p8i_one))); 00075 Packet8i sign_flip_mask = pshiftleft(shift_isodd, 31); 00076 00077 // Create a mask for which interpolant to use, i.e. if z > 1, then the mask 00078 // is set to ones for that entry. 00079 Packet8f ival_mask = _mm256_cmp_ps(z, p8f_one, _CMP_GT_OQ); 00080 00081 // Evaluate the polynomial for the interval [1,3] in z. 00082 _EIGEN_DECLARE_CONST_Packet8f(coeff_right_0, 9.999999724233232e-01f); 00083 _EIGEN_DECLARE_CONST_Packet8f(coeff_right_2, -3.084242535619928e-01f); 00084 _EIGEN_DECLARE_CONST_Packet8f(coeff_right_4, 1.584991525700324e-02f); 00085 _EIGEN_DECLARE_CONST_Packet8f(coeff_right_6, -3.188805084631342e-04f); 00086 Packet8f z_minus_two = psub(z, p8f_two); 00087 Packet8f z_minus_two2 = pmul(z_minus_two, z_minus_two); 00088 Packet8f right = pmadd(p8f_coeff_right_6, z_minus_two2, p8f_coeff_right_4); 00089 right = pmadd(right, z_minus_two2, p8f_coeff_right_2); 00090 right = pmadd(right, z_minus_two2, p8f_coeff_right_0); 00091 00092 // Evaluate the polynomial for the interval [-1,1] in z. 00093 _EIGEN_DECLARE_CONST_Packet8f(coeff_left_1, 7.853981525427295e-01f); 00094 _EIGEN_DECLARE_CONST_Packet8f(coeff_left_3, -8.074536727092352e-02f); 00095 _EIGEN_DECLARE_CONST_Packet8f(coeff_left_5, 2.489871967827018e-03f); 00096 _EIGEN_DECLARE_CONST_Packet8f(coeff_left_7, -3.587725841214251e-05f); 00097 Packet8f z2 = pmul(z, z); 00098 Packet8f left = pmadd(p8f_coeff_left_7, z2, p8f_coeff_left_5); 00099 left = pmadd(left, z2, p8f_coeff_left_3); 00100 left = pmadd(left, z2, p8f_coeff_left_1); 00101 left = pmul(left, z); 00102 00103 // Assemble the results, i.e. select the left and right polynomials. 00104 left = _mm256_andnot_ps(ival_mask, left); 00105 right = _mm256_and_ps(ival_mask, right); 00106 Packet8f res = _mm256_or_ps(left, right); 00107 00108 // Flip the sign on the odd intervals and return the result. 00109 res = _mm256_xor_ps(res, _mm256_castsi256_ps(sign_flip_mask)); 00110 return res; 00111 } 00112 00113 // Natural logarithm 00114 // Computes log(x) as log(2^e * m) = C*e + log(m), where the constant C =log(2) 00115 // and m is in the range [sqrt(1/2),sqrt(2)). In this range, the logarithm can 00116 // be easily approximated by a polynomial centered on m=1 for stability. 00117 // TODO(gonnet): Further reduce the interval allowing for lower-degree 00118 // polynomial interpolants -> ... -> profit! 00119 template <> 00120 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 00121 plog<Packet8f>(const Packet8f& _x) { 00122 Packet8f x = _x; 00123 _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f); 00124 _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f); 00125 _EIGEN_DECLARE_CONST_Packet8f(126f, 126.0f); 00126 00127 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inv_mant_mask, ~0x7f800000); 00128 00129 // The smallest non denormalized float number. 00130 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(min_norm_pos, 0x00800000); 00131 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(minus_inf, 0xff800000); 00132 00133 // Polynomial coefficients. 00134 _EIGEN_DECLARE_CONST_Packet8f(cephes_SQRTHF, 0.707106781186547524f); 00135 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p0, 7.0376836292E-2f); 00136 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p1, -1.1514610310E-1f); 00137 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p2, 1.1676998740E-1f); 00138 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p3, -1.2420140846E-1f); 00139 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p4, +1.4249322787E-1f); 00140 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p5, -1.6668057665E-1f); 00141 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p6, +2.0000714765E-1f); 00142 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p7, -2.4999993993E-1f); 00143 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p8, +3.3333331174E-1f); 00144 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q1, -2.12194440e-4f); 00145 _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q2, 0.693359375f); 00146 00147 Packet8f invalid_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_NGE_UQ); // not greater equal is true if x is NaN 00148 Packet8f iszero_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_EQ_OQ); 00149 00150 // Truncate input values to the minimum positive normal. 00151 x = pmax(x, p8f_min_norm_pos); 00152 00153 Packet8f emm0 = pshiftright(x,23); 00154 Packet8f e = _mm256_sub_ps(emm0, p8f_126f); 00155 00156 // Set the exponents to -1, i.e. x are in the range [0.5,1). 00157 x = _mm256_and_ps(x, p8f_inv_mant_mask); 00158 x = _mm256_or_ps(x, p8f_half); 00159 00160 // part2: Shift the inputs from the range [0.5,1) to [sqrt(1/2),sqrt(2)) 00161 // and shift by -1. The values are then centered around 0, which improves 00162 // the stability of the polynomial evaluation. 00163 // if( x < SQRTHF ) { 00164 // e -= 1; 00165 // x = x + x - 1.0; 00166 // } else { x = x - 1.0; } 00167 Packet8f mask = _mm256_cmp_ps(x, p8f_cephes_SQRTHF, _CMP_LT_OQ); 00168 Packet8f tmp = _mm256_and_ps(x, mask); 00169 x = psub(x, p8f_1); 00170 e = psub(e, _mm256_and_ps(p8f_1, mask)); 00171 x = padd(x, tmp); 00172 00173 Packet8f x2 = pmul(x, x); 00174 Packet8f x3 = pmul(x2, x); 00175 00176 // Evaluate the polynomial approximant of degree 8 in three parts, probably 00177 // to improve instruction-level parallelism. 00178 Packet8f y, y1, y2; 00179 y = pmadd(p8f_cephes_log_p0, x, p8f_cephes_log_p1); 00180 y1 = pmadd(p8f_cephes_log_p3, x, p8f_cephes_log_p4); 00181 y2 = pmadd(p8f_cephes_log_p6, x, p8f_cephes_log_p7); 00182 y = pmadd(y, x, p8f_cephes_log_p2); 00183 y1 = pmadd(y1, x, p8f_cephes_log_p5); 00184 y2 = pmadd(y2, x, p8f_cephes_log_p8); 00185 y = pmadd(y, x3, y1); 00186 y = pmadd(y, x3, y2); 00187 y = pmul(y, x3); 00188 00189 // Add the logarithm of the exponent back to the result of the interpolation. 00190 y1 = pmul(e, p8f_cephes_log_q1); 00191 tmp = pmul(x2, p8f_half); 00192 y = padd(y, y1); 00193 x = psub(x, tmp); 00194 y2 = pmul(e, p8f_cephes_log_q2); 00195 x = padd(x, y); 00196 x = padd(x, y2); 00197 00198 // Filter out invalid inputs, i.e. negative arg will be NAN, 0 will be -INF. 00199 return _mm256_or_ps( 00200 _mm256_andnot_ps(iszero_mask, _mm256_or_ps(x, invalid_mask)), 00201 _mm256_and_ps(iszero_mask, p8f_minus_inf)); 00202 } 00203 00204 // Exponential function. Works by writing "x = m*log(2) + r" where 00205 // "m = floor(x/log(2)+1/2)" and "r" is the remainder. The result is then 00206 // "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1). 00207 template <> 00208 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 00209 pexp<Packet8f>(const Packet8f& _x) { 00210 _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f); 00211 _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f); 00212 _EIGEN_DECLARE_CONST_Packet8f(127, 127.0f); 00213 00214 _EIGEN_DECLARE_CONST_Packet8f(exp_hi, 88.3762626647950f); 00215 _EIGEN_DECLARE_CONST_Packet8f(exp_lo, -88.3762626647949f); 00216 00217 _EIGEN_DECLARE_CONST_Packet8f(cephes_LOG2EF, 1.44269504088896341f); 00218 00219 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p0, 1.9875691500E-4f); 00220 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p1, 1.3981999507E-3f); 00221 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p2, 8.3334519073E-3f); 00222 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p3, 4.1665795894E-2f); 00223 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p4, 1.6666665459E-1f); 00224 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p5, 5.0000001201E-1f); 00225 00226 // Clamp x. 00227 Packet8f x = pmax(pmin(_x, p8f_exp_hi), p8f_exp_lo); 00228 00229 // Express exp(x) as exp(m*ln(2) + r), start by extracting 00230 // m = floor(x/ln(2) + 0.5). 00231 Packet8f m = _mm256_floor_ps(pmadd(x, p8f_cephes_LOG2EF, p8f_half)); 00232 00233 // Get r = x - m*ln(2). If no FMA instructions are available, m*ln(2) is 00234 // subtracted out in two parts, m*C1+m*C2 = m*ln(2), to avoid accumulating 00235 // truncation errors. Note that we don't use the "pmadd" function here to 00236 // ensure that a precision-preserving FMA instruction is used. 00237 #ifdef EIGEN_VECTORIZE_FMA 00238 _EIGEN_DECLARE_CONST_Packet8f(nln2, -0.6931471805599453f); 00239 Packet8f r = _mm256_fmadd_ps(m, p8f_nln2, x); 00240 #else 00241 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C1, 0.693359375f); 00242 _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C2, -2.12194440e-4f); 00243 Packet8f r = psub(x, pmul(m, p8f_cephes_exp_C1)); 00244 r = psub(r, pmul(m, p8f_cephes_exp_C2)); 00245 #endif 00246 00247 Packet8f r2 = pmul(r, r); 00248 00249 // TODO(gonnet): Split into odd/even polynomials and try to exploit 00250 // instruction-level parallelism. 00251 Packet8f y = p8f_cephes_exp_p0; 00252 y = pmadd(y, r, p8f_cephes_exp_p1); 00253 y = pmadd(y, r, p8f_cephes_exp_p2); 00254 y = pmadd(y, r, p8f_cephes_exp_p3); 00255 y = pmadd(y, r, p8f_cephes_exp_p4); 00256 y = pmadd(y, r, p8f_cephes_exp_p5); 00257 y = pmadd(y, r2, r); 00258 y = padd(y, p8f_1); 00259 00260 // Build emm0 = 2^m. 00261 Packet8i emm0 = _mm256_cvttps_epi32(padd(m, p8f_127)); 00262 emm0 = pshiftleft(emm0, 23); 00263 00264 // Return 2^m * exp(r). 00265 return pmax(pmul(y, _mm256_castsi256_ps(emm0)), _x); 00266 } 00267 00268 // Hyperbolic Tangent function. 00269 template <> 00270 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 00271 ptanh<Packet8f>(const Packet8f& x) { 00272 return internal::generic_fast_tanh_float(x); 00273 } 00274 00275 template <> 00276 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d 00277 pexp<Packet4d>(const Packet4d& _x) { 00278 Packet4d x = _x; 00279 00280 _EIGEN_DECLARE_CONST_Packet4d(1, 1.0); 00281 _EIGEN_DECLARE_CONST_Packet4d(2, 2.0); 00282 _EIGEN_DECLARE_CONST_Packet4d(half, 0.5); 00283 00284 _EIGEN_DECLARE_CONST_Packet4d(exp_hi, 709.437); 00285 _EIGEN_DECLARE_CONST_Packet4d(exp_lo, -709.436139303); 00286 00287 _EIGEN_DECLARE_CONST_Packet4d(cephes_LOG2EF, 1.4426950408889634073599); 00288 00289 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p0, 1.26177193074810590878e-4); 00290 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p1, 3.02994407707441961300e-2); 00291 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p2, 9.99999999999999999910e-1); 00292 00293 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q0, 3.00198505138664455042e-6); 00294 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q1, 2.52448340349684104192e-3); 00295 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q2, 2.27265548208155028766e-1); 00296 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q3, 2.00000000000000000009e0); 00297 00298 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C1, 0.693145751953125); 00299 _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C2, 1.42860682030941723212e-6); 00300 _EIGEN_DECLARE_CONST_Packet4i(1023, 1023); 00301 00302 Packet4d tmp, fx; 00303 00304 // clamp x 00305 x = pmax(pmin(x, p4d_exp_hi), p4d_exp_lo); 00306 // Express exp(x) as exp(g + n*log(2)). 00307 fx = pmadd(p4d_cephes_LOG2EF, x, p4d_half); 00308 00309 // Get the integer modulus of log(2), i.e. the "n" described above. 00310 fx = _mm256_floor_pd(fx); 00311 00312 // Get the remainder modulo log(2), i.e. the "g" described above. Subtract 00313 // n*log(2) out in two steps, i.e. n*C1 + n*C2, C1+C2=log2 to get the last 00314 // digits right. 00315 tmp = pmul(fx, p4d_cephes_exp_C1); 00316 Packet4d z = pmul(fx, p4d_cephes_exp_C2); 00317 x = psub(x, tmp); 00318 x = psub(x, z); 00319 00320 Packet4d x2 = pmul(x, x); 00321 00322 // Evaluate the numerator polynomial of the rational interpolant. 00323 Packet4d px = p4d_cephes_exp_p0; 00324 px = pmadd(px, x2, p4d_cephes_exp_p1); 00325 px = pmadd(px, x2, p4d_cephes_exp_p2); 00326 px = pmul(px, x); 00327 00328 // Evaluate the denominator polynomial of the rational interpolant. 00329 Packet4d qx = p4d_cephes_exp_q0; 00330 qx = pmadd(qx, x2, p4d_cephes_exp_q1); 00331 qx = pmadd(qx, x2, p4d_cephes_exp_q2); 00332 qx = pmadd(qx, x2, p4d_cephes_exp_q3); 00333 00334 // I don't really get this bit, copied from the SSE2 routines, so... 00335 // TODO(gonnet): Figure out what is going on here, perhaps find a better 00336 // rational interpolant? 00337 x = _mm256_div_pd(px, psub(qx, px)); 00338 x = pmadd(p4d_2, x, p4d_1); 00339 00340 // Build e=2^n by constructing the exponents in a 128-bit vector and 00341 // shifting them to where they belong in double-precision values. 00342 __m128i emm0 = _mm256_cvtpd_epi32(fx); 00343 emm0 = _mm_add_epi32(emm0, p4i_1023); 00344 emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(3, 1, 2, 0)); 00345 __m128i lo = _mm_slli_epi64(emm0, 52); 00346 __m128i hi = _mm_slli_epi64(_mm_srli_epi64(emm0, 32), 52); 00347 __m256i e = _mm256_insertf128_si256(_mm256_setzero_si256(), lo, 0); 00348 e = _mm256_insertf128_si256(e, hi, 1); 00349 00350 // Construct the result 2^n * exp(g) = e * x. The max is used to catch 00351 // non-finite values in the input. 00352 return pmax(pmul(x, _mm256_castsi256_pd(e)), _x); 00353 } 00354 00355 // Functions for sqrt. 00356 // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step 00357 // of Newton's method, at a cost of 1-2 bits of precision as opposed to the 00358 // exact solution. It does not handle +inf, or denormalized numbers correctly. 00359 // The main advantage of this approach is not just speed, but also the fact that 00360 // it can be inlined and pipelined with other computations, further reducing its 00361 // effective latency. This is similar to Quake3's fast inverse square root. 00362 // For detail see here: http://www.beyond3d.com/content/articles/8/ 00363 #if EIGEN_FAST_MATH 00364 template <> 00365 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 00366 psqrt<Packet8f>(const Packet8f& _x) { 00367 Packet8f half = pmul(_x, pset1<Packet8f>(.5f)); 00368 Packet8f denormal_mask = _mm256_and_ps( 00369 _mm256_cmp_ps(_x, pset1<Packet8f>((std::numeric_limits<float>::min)()), 00370 _CMP_LT_OQ), 00371 _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_GE_OQ)); 00372 00373 // Compute approximate reciprocal sqrt. 00374 Packet8f x = _mm256_rsqrt_ps(_x); 00375 // Do a single step of Newton's iteration. 00376 x = pmul(x, psub(pset1<Packet8f>(1.5f), pmul(half, pmul(x,x)))); 00377 // Flush results for denormals to zero. 00378 return _mm256_andnot_ps(denormal_mask, pmul(_x,x)); 00379 } 00380 #else 00381 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00382 Packet8f psqrt<Packet8f>(const Packet8f& x) { 00383 return _mm256_sqrt_ps(x); 00384 } 00385 #endif 00386 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00387 Packet4d psqrt<Packet4d>(const Packet4d& x) { 00388 return _mm256_sqrt_pd(x); 00389 } 00390 #if EIGEN_FAST_MATH 00391 00392 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00393 Packet8f prsqrt<Packet8f>(const Packet8f& _x) { 00394 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inf, 0x7f800000); 00395 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(nan, 0x7fc00000); 00396 _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f); 00397 _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f); 00398 _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000); 00399 00400 Packet8f neg_half = pmul(_x, p8f_minus_half); 00401 00402 // select only the inverse sqrt of positive normal inputs (denormals are 00403 // flushed to zero and cause infs as well). 00404 Packet8f le_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_LT_OQ); 00405 Packet8f x = _mm256_andnot_ps(le_zero_mask, _mm256_rsqrt_ps(_x)); 00406 00407 // Fill in NaNs and Infs for the negative/zero entries. 00408 Packet8f neg_mask = _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_LT_OQ); 00409 Packet8f zero_mask = _mm256_andnot_ps(neg_mask, le_zero_mask); 00410 Packet8f infs_and_nans = _mm256_or_ps(_mm256_and_ps(neg_mask, p8f_nan), 00411 _mm256_and_ps(zero_mask, p8f_inf)); 00412 00413 // Do a single step of Newton's iteration. 00414 x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five)); 00415 00416 // Insert NaNs and Infs in all the right places. 00417 return _mm256_or_ps(x, infs_and_nans); 00418 } 00419 00420 #else 00421 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00422 Packet8f prsqrt<Packet8f>(const Packet8f& x) { 00423 _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f); 00424 return _mm256_div_ps(p8f_one, _mm256_sqrt_ps(x)); 00425 } 00426 #endif 00427 00428 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00429 Packet4d prsqrt<Packet4d>(const Packet4d& x) { 00430 _EIGEN_DECLARE_CONST_Packet4d(one, 1.0); 00431 return _mm256_div_pd(p4d_one, _mm256_sqrt_pd(x)); 00432 } 00433 00434 00435 } // end namespace internal 00436 00437 } // end namespace Eigen 00438 00439 #endif // EIGEN_MATH_FUNCTIONS_AVX_H