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