arch/SSE/MathFunctions.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2007 Julien Pommier
5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 /* The sin, cos, exp, and log functions of this file come from
27  * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
28  */
29 
30 #ifndef EIGEN_MATH_FUNCTIONS_SSE_H
31 #define EIGEN_MATH_FUNCTIONS_SSE_H
32 
33 namespace Eigen {
34 
35 namespace internal {
36 
39 {
40  Packet4f x = _x;
44 
45  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(inv_mant_mask, ~0x7f800000);
46 
47  /* the smallest non denormalized float number */
48  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos, 0x00800000);
49 
50  /* natural logarithm computed for 4 simultaneous float
51  return NaN for x <= 0
52  */
53  _EIGEN_DECLARE_CONST_Packet4f(cephes_SQRTHF, 0.707106781186547524f);
54  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p0, 7.0376836292E-2f);
55  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p1, - 1.1514610310E-1f);
56  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p2, 1.1676998740E-1f);
57  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p3, - 1.2420140846E-1f);
58  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p4, + 1.4249322787E-1f);
59  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p5, - 1.6668057665E-1f);
60  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p6, + 2.0000714765E-1f);
61  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p7, - 2.4999993993E-1f);
62  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p8, + 3.3333331174E-1f);
63  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q1, -2.12194440e-4f);
64  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q2, 0.693359375f);
65 
66 
67  Packet4i emm0;
68 
69  Packet4f invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
70 
71  x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */
72  emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
73 
74  /* keep only the fractional part */
75  x = _mm_and_ps(x, p4f_inv_mant_mask);
76  x = _mm_or_ps(x, p4f_half);
77 
78  emm0 = _mm_sub_epi32(emm0, p4i_0x7f);
79  Packet4f e = padd(_mm_cvtepi32_ps(emm0), p4f_1);
80 
81  /* part2:
82  if( x < SQRTHF ) {
83  e -= 1;
84  x = x + x - 1.0;
85  } else { x = x - 1.0; }
86  */
87  Packet4f mask = _mm_cmplt_ps(x, p4f_cephes_SQRTHF);
88  Packet4f tmp = _mm_and_ps(x, mask);
89  x = psub(x, p4f_1);
90  e = psub(e, _mm_and_ps(p4f_1, mask));
91  x = padd(x, tmp);
92 
93  Packet4f x2 = pmul(x,x);
94  Packet4f x3 = pmul(x2,x);
95 
96  Packet4f y, y1, y2;
97  y = pmadd(p4f_cephes_log_p0, x, p4f_cephes_log_p1);
98  y1 = pmadd(p4f_cephes_log_p3, x, p4f_cephes_log_p4);
99  y2 = pmadd(p4f_cephes_log_p6, x, p4f_cephes_log_p7);
100  y = pmadd(y , x, p4f_cephes_log_p2);
101  y1 = pmadd(y1, x, p4f_cephes_log_p5);
102  y2 = pmadd(y2, x, p4f_cephes_log_p8);
103  y = pmadd(y, x3, y1);
104  y = pmadd(y, x3, y2);
105  y = pmul(y, x3);
106 
107  y1 = pmul(e, p4f_cephes_log_q1);
108  tmp = pmul(x2, p4f_half);
109  y = padd(y, y1);
110  x = psub(x, tmp);
111  y2 = pmul(e, p4f_cephes_log_q2);
112  x = padd(x, y);
113  x = padd(x, y2);
114  return _mm_or_ps(x, invalid_mask); // negative arg will be NAN
115 }
116 
119 {
120  Packet4f x = _x;
122  _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
123  _EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f);
124 
125 
126  _EIGEN_DECLARE_CONST_Packet4f(exp_hi, 88.3762626647950f);
127  _EIGEN_DECLARE_CONST_Packet4f(exp_lo, -88.3762626647949f);
128 
129  _EIGEN_DECLARE_CONST_Packet4f(cephes_LOG2EF, 1.44269504088896341f);
130  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C1, 0.693359375f);
131  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C2, -2.12194440e-4f);
132 
133  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p0, 1.9875691500E-4f);
134  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p1, 1.3981999507E-3f);
135  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p2, 8.3334519073E-3f);
136  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p3, 4.1665795894E-2f);
137  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p4, 1.6666665459E-1f);
138  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p5, 5.0000001201E-1f);
139 
140  Packet4f tmp = _mm_setzero_ps(), fx;
141  Packet4i emm0;
142 
143  // clamp x
144  x = pmax(pmin(x, p4f_exp_hi), p4f_exp_lo);
145 
146  /* express exp(x) as exp(g + n*log(2)) */
147  fx = pmadd(x, p4f_cephes_LOG2EF, p4f_half);
148 
149  /* how to perform a floorf with SSE: just below */
150  emm0 = _mm_cvttps_epi32(fx);
151  tmp = _mm_cvtepi32_ps(emm0);
152  /* if greater, substract 1 */
153  Packet4f mask = _mm_cmpgt_ps(tmp, fx);
154  mask = _mm_and_ps(mask, p4f_1);
155  fx = psub(tmp, mask);
156 
157  tmp = pmul(fx, p4f_cephes_exp_C1);
158  Packet4f z = pmul(fx, p4f_cephes_exp_C2);
159  x = psub(x, tmp);
160  x = psub(x, z);
161 
162  z = pmul(x,x);
163 
164  Packet4f y = p4f_cephes_exp_p0;
165  y = pmadd(y, x, p4f_cephes_exp_p1);
166  y = pmadd(y, x, p4f_cephes_exp_p2);
167  y = pmadd(y, x, p4f_cephes_exp_p3);
168  y = pmadd(y, x, p4f_cephes_exp_p4);
169  y = pmadd(y, x, p4f_cephes_exp_p5);
170  y = pmadd(y, z, x);
171  y = padd(y, p4f_1);
172 
173  // build 2^n
174  emm0 = _mm_cvttps_epi32(fx);
175  emm0 = _mm_add_epi32(emm0, p4i_0x7f);
176  emm0 = _mm_slli_epi32(emm0, 23);
177  return pmul(y, _mm_castsi128_ps(emm0));
178 }
179 
180 /* evaluation of 4 sines at onces, using SSE2 intrinsics.
181 
182  The code is the exact rewriting of the cephes sinf function.
183  Precision is excellent as long as x < 8192 (I did not bother to
184  take into account the special handling they have for greater values
185  -- it does not return garbage for arguments over 8192, though, but
186  the extra precision is missing).
187 
188  Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
189  surprising but correct result.
190 */
191 
194 {
195  Packet4f x = _x;
197  _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
198 
203 
204  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(sign_mask, 0x80000000);
205 
206  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP1,-0.78515625f);
207  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP2, -2.4187564849853515625e-4f);
208  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP3, -3.77489497744594108e-8f);
209  _EIGEN_DECLARE_CONST_Packet4f(sincof_p0, -1.9515295891E-4f);
210  _EIGEN_DECLARE_CONST_Packet4f(sincof_p1, 8.3321608736E-3f);
211  _EIGEN_DECLARE_CONST_Packet4f(sincof_p2, -1.6666654611E-1f);
212  _EIGEN_DECLARE_CONST_Packet4f(coscof_p0, 2.443315711809948E-005f);
213  _EIGEN_DECLARE_CONST_Packet4f(coscof_p1, -1.388731625493765E-003f);
214  _EIGEN_DECLARE_CONST_Packet4f(coscof_p2, 4.166664568298827E-002f);
215  _EIGEN_DECLARE_CONST_Packet4f(cephes_FOPI, 1.27323954473516f); // 4 / M_PI
216 
217  Packet4f xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;
218 
219  Packet4i emm0, emm2;
220  sign_bit = x;
221  /* take the absolute value */
222  x = pabs(x);
223 
224  /* take the modulo */
225 
226  /* extract the sign bit (upper one) */
227  sign_bit = _mm_and_ps(sign_bit, p4f_sign_mask);
228 
229  /* scale by 4/Pi */
230  y = pmul(x, p4f_cephes_FOPI);
231 
232  /* store the integer part of y in mm0 */
233  emm2 = _mm_cvttps_epi32(y);
234  /* j=(j+1) & (~1) (see the cephes sources) */
235  emm2 = _mm_add_epi32(emm2, p4i_1);
236  emm2 = _mm_and_si128(emm2, p4i_not1);
237  y = _mm_cvtepi32_ps(emm2);
238  /* get the swap sign flag */
239  emm0 = _mm_and_si128(emm2, p4i_4);
240  emm0 = _mm_slli_epi32(emm0, 29);
241  /* get the polynom selection mask
242  there is one polynom for 0 <= x <= Pi/4
243  and another one for Pi/4<x<=Pi/2
244 
245  Both branches will be computed.
246  */
247  emm2 = _mm_and_si128(emm2, p4i_2);
248  emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
249 
250  Packet4f swap_sign_bit = _mm_castsi128_ps(emm0);
251  Packet4f poly_mask = _mm_castsi128_ps(emm2);
252  sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
253 
254  /* The magic pass: "Extended precision modular arithmetic"
255  x = ((x - y * DP1) - y * DP2) - y * DP3; */
256  xmm1 = pmul(y, p4f_minus_cephes_DP1);
257  xmm2 = pmul(y, p4f_minus_cephes_DP2);
258  xmm3 = pmul(y, p4f_minus_cephes_DP3);
259  x = padd(x, xmm1);
260  x = padd(x, xmm2);
261  x = padd(x, xmm3);
262 
263  /* Evaluate the first polynom (0 <= x <= Pi/4) */
264  y = p4f_coscof_p0;
265  Packet4f z = _mm_mul_ps(x,x);
266 
267  y = pmadd(y, z, p4f_coscof_p1);
268  y = pmadd(y, z, p4f_coscof_p2);
269  y = pmul(y, z);
270  y = pmul(y, z);
271  Packet4f tmp = pmul(z, p4f_half);
272  y = psub(y, tmp);
273  y = padd(y, p4f_1);
274 
275  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
276 
277  Packet4f y2 = p4f_sincof_p0;
278  y2 = pmadd(y2, z, p4f_sincof_p1);
279  y2 = pmadd(y2, z, p4f_sincof_p2);
280  y2 = pmul(y2, z);
281  y2 = pmul(y2, x);
282  y2 = padd(y2, x);
283 
284  /* select the correct result from the two polynoms */
285  y2 = _mm_and_ps(poly_mask, y2);
286  y = _mm_andnot_ps(poly_mask, y);
287  y = _mm_or_ps(y,y2);
288  /* update the sign */
289  return _mm_xor_ps(y, sign_bit);
290 }
291 
292 /* almost the same as psin */
295 {
296  Packet4f x = _x;
298  _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
299 
304 
305  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP1,-0.78515625f);
306  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP2, -2.4187564849853515625e-4f);
307  _EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP3, -3.77489497744594108e-8f);
308  _EIGEN_DECLARE_CONST_Packet4f(sincof_p0, -1.9515295891E-4f);
309  _EIGEN_DECLARE_CONST_Packet4f(sincof_p1, 8.3321608736E-3f);
310  _EIGEN_DECLARE_CONST_Packet4f(sincof_p2, -1.6666654611E-1f);
311  _EIGEN_DECLARE_CONST_Packet4f(coscof_p0, 2.443315711809948E-005f);
312  _EIGEN_DECLARE_CONST_Packet4f(coscof_p1, -1.388731625493765E-003f);
313  _EIGEN_DECLARE_CONST_Packet4f(coscof_p2, 4.166664568298827E-002f);
314  _EIGEN_DECLARE_CONST_Packet4f(cephes_FOPI, 1.27323954473516f); // 4 / M_PI
315 
316  Packet4f xmm1, xmm2 = _mm_setzero_ps(), xmm3, y;
317  Packet4i emm0, emm2;
318 
319  x = pabs(x);
320 
321  /* scale by 4/Pi */
322  y = pmul(x, p4f_cephes_FOPI);
323 
324  /* get the integer part of y */
325  emm2 = _mm_cvttps_epi32(y);
326  /* j=(j+1) & (~1) (see the cephes sources) */
327  emm2 = _mm_add_epi32(emm2, p4i_1);
328  emm2 = _mm_and_si128(emm2, p4i_not1);
329  y = _mm_cvtepi32_ps(emm2);
330 
331  emm2 = _mm_sub_epi32(emm2, p4i_2);
332 
333  /* get the swap sign flag */
334  emm0 = _mm_andnot_si128(emm2, p4i_4);
335  emm0 = _mm_slli_epi32(emm0, 29);
336  /* get the polynom selection mask */
337  emm2 = _mm_and_si128(emm2, p4i_2);
338  emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
339 
340  Packet4f sign_bit = _mm_castsi128_ps(emm0);
341  Packet4f poly_mask = _mm_castsi128_ps(emm2);
342 
343  /* The magic pass: "Extended precision modular arithmetic"
344  x = ((x - y * DP1) - y * DP2) - y * DP3; */
345  xmm1 = pmul(y, p4f_minus_cephes_DP1);
346  xmm2 = pmul(y, p4f_minus_cephes_DP2);
347  xmm3 = pmul(y, p4f_minus_cephes_DP3);
348  x = padd(x, xmm1);
349  x = padd(x, xmm2);
350  x = padd(x, xmm3);
351 
352  /* Evaluate the first polynom (0 <= x <= Pi/4) */
353  y = p4f_coscof_p0;
354  Packet4f z = pmul(x,x);
355 
356  y = pmadd(y,z,p4f_coscof_p1);
357  y = pmadd(y,z,p4f_coscof_p2);
358  y = pmul(y, z);
359  y = pmul(y, z);
360  Packet4f tmp = _mm_mul_ps(z, p4f_half);
361  y = psub(y, tmp);
362  y = padd(y, p4f_1);
363 
364  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
365  Packet4f y2 = p4f_sincof_p0;
366  y2 = pmadd(y2, z, p4f_sincof_p1);
367  y2 = pmadd(y2, z, p4f_sincof_p2);
368  y2 = pmul(y2, z);
369  y2 = pmadd(y2, x, x);
370 
371  /* select the correct result from the two polynoms */
372  y2 = _mm_and_ps(poly_mask, y2);
373  y = _mm_andnot_ps(poly_mask, y);
374  y = _mm_or_ps(y,y2);
375 
376  /* update the sign */
377  return _mm_xor_ps(y, sign_bit);
378 }
379 
380 // This is based on Quake3's fast inverse square root.
381 // For detail see here: http://www.beyond3d.com/content/articles/8/
384 {
385  Packet4f half = pmul(_x, pset1<Packet4f>(.5f));
386 
387  /* select only the inverse sqrt of non-zero inputs */
388  Packet4f non_zero_mask = _mm_cmpgt_ps(_x, pset1<Packet4f>(std::numeric_limits<float>::epsilon()));
389  Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x));
390 
391  x = pmul(x, psub(pset1<Packet4f>(1.5f), pmul(half, pmul(x,x))));
392  return pmul(_x,x);
393 }
394 
395 } // end namespace internal
396 
397 } // end namespace Eigen
398 
399 #endif // EIGEN_MATH_FUNCTIONS_SSE_H