Eigen  3.3.9
 
Loading...
Searching...
No Matches
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// Copyright (C) 2016 Konstantinos Margaritis <markos@freevec.org>
7//
8// This Source Code Form is subject to the terms of the Mozilla
9// Public License v. 2.0. If a copy of the MPL was not distributed
10// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12/* The sin, cos, exp, and log functions of this file come from
13 * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
14 */
15
16#ifndef EIGEN_MATH_FUNCTIONS_ALTIVEC_H
17#define EIGEN_MATH_FUNCTIONS_ALTIVEC_H
18
19namespace Eigen {
20
21namespace internal {
22
23static _EIGEN_DECLARE_CONST_Packet4f(1 , 1.0f);
24static _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
25static _EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f);
26static _EIGEN_DECLARE_CONST_Packet4i(23, 23);
27
28static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(inv_mant_mask, ~0x7f800000);
29
30/* the smallest non denormalized float number */
31static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos, 0x00800000);
32static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf, 0xff800000); // -1.f/0.f
33static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_nan, 0xffffffff);
34
35/* natural logarithm computed for 4 simultaneous float
36 return NaN for x <= 0
37*/
38static _EIGEN_DECLARE_CONST_Packet4f(cephes_SQRTHF, 0.707106781186547524f);
39static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p0, 7.0376836292E-2f);
40static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p1, - 1.1514610310E-1f);
41static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p2, 1.1676998740E-1f);
42static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p3, - 1.2420140846E-1f);
43static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p4, + 1.4249322787E-1f);
44static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p5, - 1.6668057665E-1f);
45static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p6, + 2.0000714765E-1f);
46static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p7, - 2.4999993993E-1f);
47static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p8, + 3.3333331174E-1f);
48static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q1, -2.12194440e-4f);
49static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q2, 0.693359375f);
50
51static _EIGEN_DECLARE_CONST_Packet4f(exp_hi, 88.3762626647950f);
52static _EIGEN_DECLARE_CONST_Packet4f(exp_lo, -88.3762626647949f);
53
54static _EIGEN_DECLARE_CONST_Packet4f(cephes_LOG2EF, 1.44269504088896341f);
55static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C1, 0.693359375f);
56static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C2, -2.12194440e-4f);
57
58static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p0, 1.9875691500E-4f);
59static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p1, 1.3981999507E-3f);
60static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p2, 8.3334519073E-3f);
61static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p3, 4.1665795894E-2f);
62static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p4, 1.6666665459E-1f);
63static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p5, 5.0000001201E-1f);
64
65#ifdef __VSX__
66static _EIGEN_DECLARE_CONST_Packet2d(1 , 1.0);
67static _EIGEN_DECLARE_CONST_Packet2d(2 , 2.0);
68static _EIGEN_DECLARE_CONST_Packet2d(half, 0.5);
69
70static _EIGEN_DECLARE_CONST_Packet2d(exp_hi, 709.437);
71static _EIGEN_DECLARE_CONST_Packet2d(exp_lo, -709.436139303);
72
73static _EIGEN_DECLARE_CONST_Packet2d(cephes_LOG2EF, 1.4426950408889634073599);
74
75static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p0, 1.26177193074810590878e-4);
76static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p1, 3.02994407707441961300e-2);
77static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p2, 9.99999999999999999910e-1);
78
79static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q0, 3.00198505138664455042e-6);
80static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q1, 2.52448340349684104192e-3);
81static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q2, 2.27265548208155028766e-1);
82static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q3, 2.00000000000000000009e0);
83
84static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C1, 0.693145751953125);
85static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C2, 1.42860682030941723212e-6);
86
87#ifdef __POWER8_VECTOR__
88static Packet2l p2l_1023 = { 1023, 1023 };
89static Packet2ul p2ul_52 = { 52, 52 };
90#endif
91
92#endif
93
94template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
95Packet4f plog<Packet4f>(const Packet4f& _x)
96{
97 Packet4f x = _x;
98
99 Packet4i emm0;
100
101 /* isvalid_mask is 0 if x < 0 or x is NaN. */
102 Packet4ui isvalid_mask = reinterpret_cast<Packet4ui>(vec_cmpge(x, p4f_ZERO));
103 Packet4ui iszero_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(x, p4f_ZERO));
104
105 x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */
106 emm0 = vec_sr(reinterpret_cast<Packet4i>(x),
107 reinterpret_cast<Packet4ui>(p4i_23));
108
109 /* keep only the fractional part */
110 x = pand(x, p4f_inv_mant_mask);
111 x = por(x, p4f_half);
112
113 emm0 = psub(emm0, p4i_0x7f);
114 Packet4f e = padd(vec_ctf(emm0, 0), p4f_1);
115
116 /* part2:
117 if( x < SQRTHF ) {
118 e -= 1;
119 x = x + x - 1.0;
120 } else { x = x - 1.0; }
121 */
122 Packet4f mask = reinterpret_cast<Packet4f>(vec_cmplt(x, p4f_cephes_SQRTHF));
123 Packet4f tmp = pand(x, mask);
124 x = psub(x, p4f_1);
125 e = psub(e, pand(p4f_1, mask));
126 x = padd(x, tmp);
127
128 Packet4f x2 = pmul(x,x);
129 Packet4f x3 = pmul(x2,x);
130
131 Packet4f y, y1, y2;
132 y = pmadd(p4f_cephes_log_p0, x, p4f_cephes_log_p1);
133 y1 = pmadd(p4f_cephes_log_p3, x, p4f_cephes_log_p4);
134 y2 = pmadd(p4f_cephes_log_p6, x, p4f_cephes_log_p7);
135 y = pmadd(y , x, p4f_cephes_log_p2);
136 y1 = pmadd(y1, x, p4f_cephes_log_p5);
137 y2 = pmadd(y2, x, p4f_cephes_log_p8);
138 y = pmadd(y, x3, y1);
139 y = pmadd(y, x3, y2);
140 y = pmul(y, x3);
141
142 y1 = pmul(e, p4f_cephes_log_q1);
143 tmp = pmul(x2, p4f_half);
144 y = padd(y, y1);
145 x = psub(x, tmp);
146 y2 = pmul(e, p4f_cephes_log_q2);
147 x = padd(x, y);
148 x = padd(x, y2);
149 // negative arg will be NAN, 0 will be -INF
150 x = vec_sel(x, p4f_minus_inf, iszero_mask);
151 x = vec_sel(p4f_minus_nan, x, isvalid_mask);
152 return x;
153}
154
155template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
156Packet4f pexp<Packet4f>(const Packet4f& _x)
157{
158 Packet4f x = _x;
159
160 Packet4f tmp, fx;
161 Packet4i emm0;
162
163 // clamp x
164 x = pmax(pmin(x, p4f_exp_hi), p4f_exp_lo);
165
166 // express exp(x) as exp(g + n*log(2))
167 fx = pmadd(x, p4f_cephes_LOG2EF, p4f_half);
168
169 fx = pfloor(fx);
170
171 tmp = pmul(fx, p4f_cephes_exp_C1);
172 Packet4f z = pmul(fx, p4f_cephes_exp_C2);
173 x = psub(x, tmp);
174 x = psub(x, z);
175
176 z = pmul(x,x);
177
178 Packet4f y = p4f_cephes_exp_p0;
179 y = pmadd(y, x, p4f_cephes_exp_p1);
180 y = pmadd(y, x, p4f_cephes_exp_p2);
181 y = pmadd(y, x, p4f_cephes_exp_p3);
182 y = pmadd(y, x, p4f_cephes_exp_p4);
183 y = pmadd(y, x, p4f_cephes_exp_p5);
184 y = pmadd(y, z, x);
185 y = padd(y, p4f_1);
186
187 // build 2^n
188 emm0 = vec_cts(fx, 0);
189 emm0 = vec_add(emm0, p4i_0x7f);
190 emm0 = vec_sl(emm0, reinterpret_cast<Packet4ui>(p4i_23));
191
192 // Altivec's max & min operators just drop silent NaNs. Check NaNs in
193 // inputs and return them unmodified.
194 Packet4ui isnumber_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(_x, _x));
195 return vec_sel(_x, pmax(pmul(y, reinterpret_cast<Packet4f>(emm0)), _x),
196 isnumber_mask);
197}
198
199#ifndef EIGEN_COMP_CLANG
200template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
201Packet4f prsqrt<Packet4f>(const Packet4f& x)
202{
203 return vec_rsqrt(x);
204}
205#endif
206
207#ifdef __VSX__
208#ifndef EIGEN_COMP_CLANG
209template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
210Packet2d prsqrt<Packet2d>(const Packet2d& x)
211{
212 return vec_rsqrt(x);
213}
214#endif
215
216template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
217Packet4f psqrt<Packet4f>(const Packet4f& x)
218{
219 return vec_sqrt(x);
220}
221
222template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
223Packet2d psqrt<Packet2d>(const Packet2d& x)
224{
225 return vec_sqrt(x);
226}
227
228// VSX support varies between different compilers and even different
229// versions of the same compiler. For gcc version >= 4.9.3, we can use
230// vec_cts to efficiently convert Packet2d to Packet2l. Otherwise, use
231// a slow version that works with older compilers.
232// Update: apparently vec_cts/vec_ctf intrinsics for 64-bit doubles
233// are buggy, https://gcc.gnu.org/bugzilla/show_bug.cgi?id=70963
234static inline Packet2l ConvertToPacket2l(const Packet2d& x) {
235#if EIGEN_GNUC_AT_LEAST(5, 4) || \
236 (EIGEN_GNUC_AT(6, 1) && __GNUC_PATCHLEVEL__ >= 1)
237 return vec_cts(x, 0); // TODO: check clang version.
238#else
239 double tmp[2];
240 memcpy(tmp, &x, sizeof(tmp));
241 Packet2l l = { static_cast<long long>(tmp[0]),
242 static_cast<long long>(tmp[1]) };
243 return l;
244#endif
245}
246
247template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
248Packet2d pexp<Packet2d>(const Packet2d& _x)
249{
250 Packet2d x = _x;
251
252 Packet2d tmp, fx;
253 Packet2l emm0;
254
255 // clamp x
256 x = pmax(pmin(x, p2d_exp_hi), p2d_exp_lo);
257
258 /* express exp(x) as exp(g + n*log(2)) */
259 fx = pmadd(x, p2d_cephes_LOG2EF, p2d_half);
260
261 fx = pfloor(fx);
262
263 tmp = pmul(fx, p2d_cephes_exp_C1);
264 Packet2d z = pmul(fx, p2d_cephes_exp_C2);
265 x = psub(x, tmp);
266 x = psub(x, z);
267
268 Packet2d x2 = pmul(x,x);
269
270 Packet2d px = p2d_cephes_exp_p0;
271 px = pmadd(px, x2, p2d_cephes_exp_p1);
272 px = pmadd(px, x2, p2d_cephes_exp_p2);
273 px = pmul (px, x);
274
275 Packet2d qx = p2d_cephes_exp_q0;
276 qx = pmadd(qx, x2, p2d_cephes_exp_q1);
277 qx = pmadd(qx, x2, p2d_cephes_exp_q2);
278 qx = pmadd(qx, x2, p2d_cephes_exp_q3);
279
280 x = pdiv(px,psub(qx,px));
281 x = pmadd(p2d_2,x,p2d_1);
282
283 // build 2^n
284 emm0 = ConvertToPacket2l(fx);
285
286#ifdef __POWER8_VECTOR__
287 emm0 = vec_add(emm0, p2l_1023);
288 emm0 = vec_sl(emm0, p2ul_52);
289#else
290 // Code is a bit complex for POWER7. There is actually a
291 // vec_xxsldi intrinsic but it is not supported by some gcc versions.
292 // So we shift (52-32) bits and do a word swap with zeros.
293 _EIGEN_DECLARE_CONST_Packet4i(1023, 1023);
294 _EIGEN_DECLARE_CONST_Packet4i(20, 20); // 52 - 32
295
296 Packet4i emm04i = reinterpret_cast<Packet4i>(emm0);
297 emm04i = vec_add(emm04i, p4i_1023);
298 emm04i = vec_sl(emm04i, reinterpret_cast<Packet4ui>(p4i_20));
299 static const Packet16uc perm = {
300 0x14, 0x15, 0x16, 0x17, 0x00, 0x01, 0x02, 0x03,
301 0x1c, 0x1d, 0x1e, 0x1f, 0x08, 0x09, 0x0a, 0x0b };
302#ifdef _BIG_ENDIAN
303 emm0 = reinterpret_cast<Packet2l>(vec_perm(p4i_ZERO, emm04i, perm));
304#else
305 emm0 = reinterpret_cast<Packet2l>(vec_perm(emm04i, p4i_ZERO, perm));
306#endif
307
308#endif
309
310 // Altivec's max & min operators just drop silent NaNs. Check NaNs in
311 // inputs and return them unmodified.
312 Packet2ul isnumber_mask = reinterpret_cast<Packet2ul>(vec_cmpeq(_x, _x));
313 return vec_sel(_x, pmax(pmul(x, reinterpret_cast<Packet2d>(emm0)), _x),
314 isnumber_mask);
315}
316#endif
317
318} // end namespace internal
319
320} // end namespace Eigen
321
322#endif // EIGEN_MATH_FUNCTIONS_ALTIVEC_H
Eigen
Namespace containing all symbols from the Eigen library.
Definition A05_PortingFrom2To3.dox:1