Eigen  3.3.3
MathFunctions.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends