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