Eigen  3.2.93
arch/AVX/MathFunctions.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com)
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_MATH_FUNCTIONS_AVX_H
11 #define EIGEN_MATH_FUNCTIONS_AVX_H
12 
13 /* The sin, cos, exp, and log functions of this file are loosely derived from
14  * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
15  */
16 
17 namespace Eigen {
18 
19 namespace internal {
20 
21 inline Packet8i pshiftleft(Packet8i v, int n)
22 {
23 #ifdef EIGEN_VECTORIZE_AVX2
24  return _mm256_slli_epi32(v, n);
25 #else
26  __m128i lo = _mm_slli_epi32(_mm256_extractf128_si256(v, 0), n);
27  __m128i hi = _mm_slli_epi32(_mm256_extractf128_si256(v, 1), n);
28  return _mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1);
29 #endif
30 }
31 
32 inline Packet8f pshiftright(Packet8f v, int n)
33 {
34 #ifdef EIGEN_VECTORIZE_AVX2
35  return _mm256_cvtepi32_ps(_mm256_srli_epi32(_mm256_castps_si256(v), n));
36 #else
37  __m128i lo = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 0), n);
38  __m128i hi = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 1), n);
39  return _mm256_cvtepi32_ps(_mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1));
40 #endif
41 }
42 
43 // Sine function
44 // Computes sin(x) by wrapping x to the interval [-Pi/4,3*Pi/4] and
45 // evaluating interpolants in [-Pi/4,Pi/4] or [Pi/4,3*Pi/4]. The interpolants
46 // are (anti-)symmetric and thus have only odd/even coefficients
47 template <>
48 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
49 psin<Packet8f>(const Packet8f& _x) {
50  Packet8f x = _x;
51 
52  // Some useful values.
53  _EIGEN_DECLARE_CONST_Packet8i(one, 1);
54  _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f);
55  _EIGEN_DECLARE_CONST_Packet8f(two, 2.0f);
56  _EIGEN_DECLARE_CONST_Packet8f(one_over_four, 0.25f);
57  _EIGEN_DECLARE_CONST_Packet8f(one_over_pi, 3.183098861837907e-01f);
58  _EIGEN_DECLARE_CONST_Packet8f(neg_pi_first, -3.140625000000000e+00f);
59  _EIGEN_DECLARE_CONST_Packet8f(neg_pi_second, -9.670257568359375e-04f);
60  _EIGEN_DECLARE_CONST_Packet8f(neg_pi_third, -6.278329571784980e-07f);
61  _EIGEN_DECLARE_CONST_Packet8f(four_over_pi, 1.273239544735163e+00f);
62 
63  // Map x from [-Pi/4,3*Pi/4] to z in [-1,3] and subtract the shifted period.
64  Packet8f z = pmul(x, p8f_one_over_pi);
65  Packet8f shift = _mm256_floor_ps(padd(z, p8f_one_over_four));
66  x = pmadd(shift, p8f_neg_pi_first, x);
67  x = pmadd(shift, p8f_neg_pi_second, x);
68  x = pmadd(shift, p8f_neg_pi_third, x);
69  z = pmul(x, p8f_four_over_pi);
70 
71  // Make a mask for the entries that need flipping, i.e. wherever the shift
72  // is odd.
73  Packet8i shift_ints = _mm256_cvtps_epi32(shift);
74  Packet8i shift_isodd = _mm256_castps_si256(_mm256_and_ps(_mm256_castsi256_ps(shift_ints), _mm256_castsi256_ps(p8i_one)));
75  Packet8i sign_flip_mask = pshiftleft(shift_isodd, 31);
76 
77  // Create a mask for which interpolant to use, i.e. if z > 1, then the mask
78  // is set to ones for that entry.
79  Packet8f ival_mask = _mm256_cmp_ps(z, p8f_one, _CMP_GT_OQ);
80 
81  // Evaluate the polynomial for the interval [1,3] in z.
82  _EIGEN_DECLARE_CONST_Packet8f(coeff_right_0, 9.999999724233232e-01f);
83  _EIGEN_DECLARE_CONST_Packet8f(coeff_right_2, -3.084242535619928e-01f);
84  _EIGEN_DECLARE_CONST_Packet8f(coeff_right_4, 1.584991525700324e-02f);
85  _EIGEN_DECLARE_CONST_Packet8f(coeff_right_6, -3.188805084631342e-04f);
86  Packet8f z_minus_two = psub(z, p8f_two);
87  Packet8f z_minus_two2 = pmul(z_minus_two, z_minus_two);
88  Packet8f right = pmadd(p8f_coeff_right_6, z_minus_two2, p8f_coeff_right_4);
89  right = pmadd(right, z_minus_two2, p8f_coeff_right_2);
90  right = pmadd(right, z_minus_two2, p8f_coeff_right_0);
91 
92  // Evaluate the polynomial for the interval [-1,1] in z.
93  _EIGEN_DECLARE_CONST_Packet8f(coeff_left_1, 7.853981525427295e-01f);
94  _EIGEN_DECLARE_CONST_Packet8f(coeff_left_3, -8.074536727092352e-02f);
95  _EIGEN_DECLARE_CONST_Packet8f(coeff_left_5, 2.489871967827018e-03f);
96  _EIGEN_DECLARE_CONST_Packet8f(coeff_left_7, -3.587725841214251e-05f);
97  Packet8f z2 = pmul(z, z);
98  Packet8f left = pmadd(p8f_coeff_left_7, z2, p8f_coeff_left_5);
99  left = pmadd(left, z2, p8f_coeff_left_3);
100  left = pmadd(left, z2, p8f_coeff_left_1);
101  left = pmul(left, z);
102 
103  // Assemble the results, i.e. select the left and right polynomials.
104  left = _mm256_andnot_ps(ival_mask, left);
105  right = _mm256_and_ps(ival_mask, right);
106  Packet8f res = _mm256_or_ps(left, right);
107 
108  // Flip the sign on the odd intervals and return the result.
109  res = _mm256_xor_ps(res, _mm256_castsi256_ps(sign_flip_mask));
110  return res;
111 }
112 
113 // Natural logarithm
114 // Computes log(x) as log(2^e * m) = C*e + log(m), where the constant C =log(2)
115 // and m is in the range [sqrt(1/2),sqrt(2)). In this range, the logarithm can
116 // be easily approximated by a polynomial centered on m=1 for stability.
117 // TODO(gonnet): Further reduce the interval allowing for lower-degree
118 // polynomial interpolants -> ... -> profit!
119 template <>
120 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
121 plog<Packet8f>(const Packet8f& _x) {
122  Packet8f x = _x;
123  _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f);
124  _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f);
125  _EIGEN_DECLARE_CONST_Packet8f(126f, 126.0f);
126 
127  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inv_mant_mask, ~0x7f800000);
128 
129  // The smallest non denormalized float number.
130  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(min_norm_pos, 0x00800000);
131  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(minus_inf, 0xff800000);
132 
133  // Polynomial coefficients.
134  _EIGEN_DECLARE_CONST_Packet8f(cephes_SQRTHF, 0.707106781186547524f);
135  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p0, 7.0376836292E-2f);
136  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p1, -1.1514610310E-1f);
137  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p2, 1.1676998740E-1f);
138  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p3, -1.2420140846E-1f);
139  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p4, +1.4249322787E-1f);
140  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p5, -1.6668057665E-1f);
141  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p6, +2.0000714765E-1f);
142  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p7, -2.4999993993E-1f);
143  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p8, +3.3333331174E-1f);
144  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q1, -2.12194440e-4f);
145  _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q2, 0.693359375f);
146 
147  Packet8f invalid_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_NGE_UQ); // not greater equal is true if x is NaN
148  Packet8f iszero_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_EQ_OQ);
149 
150  // Truncate input values to the minimum positive normal.
151  x = pmax(x, p8f_min_norm_pos);
152 
153  Packet8f emm0 = pshiftright(x,23);
154  Packet8f e = _mm256_sub_ps(emm0, p8f_126f);
155 
156  // Set the exponents to -1, i.e. x are in the range [0.5,1).
157  x = _mm256_and_ps(x, p8f_inv_mant_mask);
158  x = _mm256_or_ps(x, p8f_half);
159 
160  // part2: Shift the inputs from the range [0.5,1) to [sqrt(1/2),sqrt(2))
161  // and shift by -1. The values are then centered around 0, which improves
162  // the stability of the polynomial evaluation.
163  // if( x < SQRTHF ) {
164  // e -= 1;
165  // x = x + x - 1.0;
166  // } else { x = x - 1.0; }
167  Packet8f mask = _mm256_cmp_ps(x, p8f_cephes_SQRTHF, _CMP_LT_OQ);
168  Packet8f tmp = _mm256_and_ps(x, mask);
169  x = psub(x, p8f_1);
170  e = psub(e, _mm256_and_ps(p8f_1, mask));
171  x = padd(x, tmp);
172 
173  Packet8f x2 = pmul(x, x);
174  Packet8f x3 = pmul(x2, x);
175 
176  // Evaluate the polynomial approximant of degree 8 in three parts, probably
177  // to improve instruction-level parallelism.
178  Packet8f y, y1, y2;
179  y = pmadd(p8f_cephes_log_p0, x, p8f_cephes_log_p1);
180  y1 = pmadd(p8f_cephes_log_p3, x, p8f_cephes_log_p4);
181  y2 = pmadd(p8f_cephes_log_p6, x, p8f_cephes_log_p7);
182  y = pmadd(y, x, p8f_cephes_log_p2);
183  y1 = pmadd(y1, x, p8f_cephes_log_p5);
184  y2 = pmadd(y2, x, p8f_cephes_log_p8);
185  y = pmadd(y, x3, y1);
186  y = pmadd(y, x3, y2);
187  y = pmul(y, x3);
188 
189  // Add the logarithm of the exponent back to the result of the interpolation.
190  y1 = pmul(e, p8f_cephes_log_q1);
191  tmp = pmul(x2, p8f_half);
192  y = padd(y, y1);
193  x = psub(x, tmp);
194  y2 = pmul(e, p8f_cephes_log_q2);
195  x = padd(x, y);
196  x = padd(x, y2);
197 
198  // Filter out invalid inputs, i.e. negative arg will be NAN, 0 will be -INF.
199  return _mm256_or_ps(
200  _mm256_andnot_ps(iszero_mask, _mm256_or_ps(x, invalid_mask)),
201  _mm256_and_ps(iszero_mask, p8f_minus_inf));
202 }
203 
204 // Exponential function. Works by writing "x = m*log(2) + r" where
205 // "m = floor(x/log(2)+1/2)" and "r" is the remainder. The result is then
206 // "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1).
207 template <>
208 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
209 pexp<Packet8f>(const Packet8f& _x) {
210  _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f);
211  _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f);
212  _EIGEN_DECLARE_CONST_Packet8f(127, 127.0f);
213 
214  _EIGEN_DECLARE_CONST_Packet8f(exp_hi, 88.3762626647950f);
215  _EIGEN_DECLARE_CONST_Packet8f(exp_lo, -88.3762626647949f);
216 
217  _EIGEN_DECLARE_CONST_Packet8f(cephes_LOG2EF, 1.44269504088896341f);
218 
219  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p0, 1.9875691500E-4f);
220  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p1, 1.3981999507E-3f);
221  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p2, 8.3334519073E-3f);
222  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p3, 4.1665795894E-2f);
223  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p4, 1.6666665459E-1f);
224  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p5, 5.0000001201E-1f);
225 
226  // Clamp x.
227  Packet8f x = pmax(pmin(_x, p8f_exp_hi), p8f_exp_lo);
228 
229  // Express exp(x) as exp(m*ln(2) + r), start by extracting
230  // m = floor(x/ln(2) + 0.5).
231  Packet8f m = _mm256_floor_ps(pmadd(x, p8f_cephes_LOG2EF, p8f_half));
232 
233 // Get r = x - m*ln(2). If no FMA instructions are available, m*ln(2) is
234 // subtracted out in two parts, m*C1+m*C2 = m*ln(2), to avoid accumulating
235 // truncation errors. Note that we don't use the "pmadd" function here to
236 // ensure that a precision-preserving FMA instruction is used.
237 #ifdef EIGEN_VECTORIZE_FMA
238  _EIGEN_DECLARE_CONST_Packet8f(nln2, -0.6931471805599453f);
239  Packet8f r = _mm256_fmadd_ps(m, p8f_nln2, x);
240 #else
241  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C1, 0.693359375f);
242  _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C2, -2.12194440e-4f);
243  Packet8f r = psub(x, pmul(m, p8f_cephes_exp_C1));
244  r = psub(r, pmul(m, p8f_cephes_exp_C2));
245 #endif
246 
247  Packet8f r2 = pmul(r, r);
248 
249  // TODO(gonnet): Split into odd/even polynomials and try to exploit
250  // instruction-level parallelism.
251  Packet8f y = p8f_cephes_exp_p0;
252  y = pmadd(y, r, p8f_cephes_exp_p1);
253  y = pmadd(y, r, p8f_cephes_exp_p2);
254  y = pmadd(y, r, p8f_cephes_exp_p3);
255  y = pmadd(y, r, p8f_cephes_exp_p4);
256  y = pmadd(y, r, p8f_cephes_exp_p5);
257  y = pmadd(y, r2, r);
258  y = padd(y, p8f_1);
259 
260  // Build emm0 = 2^m.
261  Packet8i emm0 = _mm256_cvttps_epi32(padd(m, p8f_127));
262  emm0 = pshiftleft(emm0, 23);
263 
264  // Return 2^m * exp(r).
265  return pmax(pmul(y, _mm256_castsi256_ps(emm0)), _x);
266 }
267 
268 // Hyperbolic Tangent function.
269 // Doesn't do anything fancy, just a 13/6-degree rational interpolant which
270 // is accurate up to a couple of ulp in the range [-9, 9], outside of which the
271 // fl(tanh(x)) = +/-1.
272 template <>
273 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
274 ptanh<Packet8f>(const Packet8f& _x) {
275  // Clamp the inputs to the range [-9, 9] since anything outside
276  // this range is +/-1.0f in single-precision.
277  _EIGEN_DECLARE_CONST_Packet8f(plus_9, 9.0f);
278  _EIGEN_DECLARE_CONST_Packet8f(minus_9, -9.0f);
279  const Packet8f x = pmax(p8f_minus_9, pmin(p8f_plus_9, _x));
280 
281  // The monomial coefficients of the numerator polynomial (odd).
282  _EIGEN_DECLARE_CONST_Packet8f(alpha_1, 4.89352455891786e-03f);
283  _EIGEN_DECLARE_CONST_Packet8f(alpha_3, 6.37261928875436e-04f);
284  _EIGEN_DECLARE_CONST_Packet8f(alpha_5, 1.48572235717979e-05f);
285  _EIGEN_DECLARE_CONST_Packet8f(alpha_7, 5.12229709037114e-08f);
286  _EIGEN_DECLARE_CONST_Packet8f(alpha_9, -8.60467152213735e-11f);
287  _EIGEN_DECLARE_CONST_Packet8f(alpha_11, 2.00018790482477e-13f);
288  _EIGEN_DECLARE_CONST_Packet8f(alpha_13, -2.76076847742355e-16f);
289 
290  // The monomial coefficients of the denominator polynomial (even).
291  _EIGEN_DECLARE_CONST_Packet8f(beta_0, 4.89352518554385e-03f);
292  _EIGEN_DECLARE_CONST_Packet8f(beta_2, 2.26843463243900e-03f);
293  _EIGEN_DECLARE_CONST_Packet8f(beta_4, 1.18534705686654e-04f);
294  _EIGEN_DECLARE_CONST_Packet8f(beta_6, 1.19825839466702e-06f);
295 
296  // Since the polynomials are odd/even, we need x^2.
297  const Packet8f x2 = pmul(x, x);
298 
299  // Evaluate the numerator polynomial p.
300  Packet8f p = pmadd(x2, p8f_alpha_13, p8f_alpha_11);
301  p = pmadd(x2, p, p8f_alpha_9);
302  p = pmadd(x2, p, p8f_alpha_7);
303  p = pmadd(x2, p, p8f_alpha_5);
304  p = pmadd(x2, p, p8f_alpha_3);
305  p = pmadd(x2, p, p8f_alpha_1);
306  p = pmul(x, p);
307 
308  // Evaluate the denominator polynomial p.
309  Packet8f q = pmadd(x2, p8f_beta_6, p8f_beta_4);
310  q = pmadd(x2, q, p8f_beta_2);
311  q = pmadd(x2, q, p8f_beta_0);
312 
313  // Divide the numerator by the denominator.
314  return pdiv(p, q);
315 }
316 
317 template <>
318 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d
319 pexp<Packet4d>(const Packet4d& _x) {
320  Packet4d x = _x;
321 
322  _EIGEN_DECLARE_CONST_Packet4d(1, 1.0);
323  _EIGEN_DECLARE_CONST_Packet4d(2, 2.0);
324  _EIGEN_DECLARE_CONST_Packet4d(half, 0.5);
325 
326  _EIGEN_DECLARE_CONST_Packet4d(exp_hi, 709.437);
327  _EIGEN_DECLARE_CONST_Packet4d(exp_lo, -709.436139303);
328 
329  _EIGEN_DECLARE_CONST_Packet4d(cephes_LOG2EF, 1.4426950408889634073599);
330 
331  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p0, 1.26177193074810590878e-4);
332  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p1, 3.02994407707441961300e-2);
333  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p2, 9.99999999999999999910e-1);
334 
335  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q0, 3.00198505138664455042e-6);
336  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q1, 2.52448340349684104192e-3);
337  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q2, 2.27265548208155028766e-1);
338  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q3, 2.00000000000000000009e0);
339 
340  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C1, 0.693145751953125);
341  _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C2, 1.42860682030941723212e-6);
342  _EIGEN_DECLARE_CONST_Packet4i(1023, 1023);
343 
344  Packet4d tmp, fx;
345 
346  // clamp x
347  x = pmax(pmin(x, p4d_exp_hi), p4d_exp_lo);
348  // Express exp(x) as exp(g + n*log(2)).
349  fx = pmadd(p4d_cephes_LOG2EF, x, p4d_half);
350 
351  // Get the integer modulus of log(2), i.e. the "n" described above.
352  fx = _mm256_floor_pd(fx);
353 
354  // Get the remainder modulo log(2), i.e. the "g" described above. Subtract
355  // n*log(2) out in two steps, i.e. n*C1 + n*C2, C1+C2=log2 to get the last
356  // digits right.
357  tmp = pmul(fx, p4d_cephes_exp_C1);
358  Packet4d z = pmul(fx, p4d_cephes_exp_C2);
359  x = psub(x, tmp);
360  x = psub(x, z);
361 
362  Packet4d x2 = pmul(x, x);
363 
364  // Evaluate the numerator polynomial of the rational interpolant.
365  Packet4d px = p4d_cephes_exp_p0;
366  px = pmadd(px, x2, p4d_cephes_exp_p1);
367  px = pmadd(px, x2, p4d_cephes_exp_p2);
368  px = pmul(px, x);
369 
370  // Evaluate the denominator polynomial of the rational interpolant.
371  Packet4d qx = p4d_cephes_exp_q0;
372  qx = pmadd(qx, x2, p4d_cephes_exp_q1);
373  qx = pmadd(qx, x2, p4d_cephes_exp_q2);
374  qx = pmadd(qx, x2, p4d_cephes_exp_q3);
375 
376  // I don't really get this bit, copied from the SSE2 routines, so...
377  // TODO(gonnet): Figure out what is going on here, perhaps find a better
378  // rational interpolant?
379  x = _mm256_div_pd(px, psub(qx, px));
380  x = pmadd(p4d_2, x, p4d_1);
381 
382  // Build e=2^n by constructing the exponents in a 128-bit vector and
383  // shifting them to where they belong in double-precision values.
384  __m128i emm0 = _mm256_cvtpd_epi32(fx);
385  emm0 = _mm_add_epi32(emm0, p4i_1023);
386  emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(3, 1, 2, 0));
387  __m128i lo = _mm_slli_epi64(emm0, 52);
388  __m128i hi = _mm_slli_epi64(_mm_srli_epi64(emm0, 32), 52);
389  __m256i e = _mm256_insertf128_si256(_mm256_setzero_si256(), lo, 0);
390  e = _mm256_insertf128_si256(e, hi, 1);
391 
392  // Construct the result 2^n * exp(g) = e * x. The max is used to catch
393  // non-finite values in the input.
394  return pmax(pmul(x, _mm256_castsi256_pd(e)), _x);
395 }
396 
397 // Functions for sqrt.
398 // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step
399 // of Newton's method, at a cost of 1-2 bits of precision as opposed to the
400 // exact solution. The main advantage of this approach is not just speed, but
401 // also the fact that it can be inlined and pipelined with other computations,
402 // further reducing its effective latency.
403 #if EIGEN_FAST_MATH
404 template <>
405 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
406 psqrt<Packet8f>(const Packet8f& _x) {
407  _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f);
408  _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f);
409  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000);
410 
411  Packet8f neg_half = pmul(_x, p8f_minus_half);
412 
413  // select only the inverse sqrt of positive normal inputs (denormals are
414  // flushed to zero and cause infs as well).
415  Packet8f non_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_GE_OQ);
416  Packet8f x = _mm256_and_ps(non_zero_mask, _mm256_rsqrt_ps(_x));
417 
418  // Do a single step of Newton's iteration.
419  x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five));
420 
421  // Multiply the original _x by it's reciprocal square root to extract the
422  // square root.
423  return pmul(_x, x);
424 }
425 #else
426 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
427 Packet8f psqrt<Packet8f>(const Packet8f& x) {
428  return _mm256_sqrt_ps(x);
429 }
430 #endif
431 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
432 Packet4d psqrt<Packet4d>(const Packet4d& x) {
433  return _mm256_sqrt_pd(x);
434 }
435 #if EIGEN_FAST_MATH
436 
437 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
438 Packet8f prsqrt<Packet8f>(const Packet8f& _x) {
439  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inf, 0x7f800000);
440  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(nan, 0x7fc00000);
441  _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f);
442  _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f);
443  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000);
444 
445  Packet8f neg_half = pmul(_x, p8f_minus_half);
446 
447  // select only the inverse sqrt of positive normal inputs (denormals are
448  // flushed to zero and cause infs as well).
449  Packet8f le_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_LT_OQ);
450  Packet8f x = _mm256_andnot_ps(le_zero_mask, _mm256_rsqrt_ps(_x));
451 
452  // Fill in NaNs and Infs for the negative/zero entries.
453  Packet8f neg_mask = _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_LT_OQ);
454  Packet8f zero_mask = _mm256_andnot_ps(neg_mask, le_zero_mask);
455  Packet8f infs_and_nans = _mm256_or_ps(_mm256_and_ps(neg_mask, p8f_nan),
456  _mm256_and_ps(zero_mask, p8f_inf));
457 
458  // Do a single step of Newton's iteration.
459  x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five));
460 
461  // Insert NaNs and Infs in all the right places.
462  return _mm256_or_ps(x, infs_and_nans);
463 }
464 
465 #else
466 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
467 Packet8f prsqrt<Packet8f>(const Packet8f& x) {
468  _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f);
469  return _mm256_div_ps(p8f_one, _mm256_sqrt_ps(x));
470 }
471 #endif
472 
473 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
474 Packet4d prsqrt<Packet4d>(const Packet4d& x) {
475  _EIGEN_DECLARE_CONST_Packet4d(one, 1.0);
476  return _mm256_div_pd(p4d_one, _mm256_sqrt_pd(x));
477 }
478 
479 
480 } // end namespace internal
481 
482 } // end namespace Eigen
483 
484 #endif // EIGEN_MATH_FUNCTIONS_AVX_H
Namespace containing all symbols from the Eigen library.
Definition: Core:271
Definition: Eigen_Colamd.h:50