Eigen  3.2.92
arch/AltiVec/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_ALTIVEC_H
16 #define EIGEN_MATH_FUNCTIONS_ALTIVEC_H
17 
18 namespace Eigen {
19 
20 namespace internal {
21 
22 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
23 Packet4f 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  _EIGEN_DECLARE_CONST_Packet4i(23, 23);
30 
31  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(inv_mant_mask, ~0x7f800000);
32 
33  /* the smallest non denormalized float number */
34  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos, 0x00800000);
35  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf, 0xff800000); // -1.f/0.f
36  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_nan, 0xffffffff);
37 
38  /* natural logarithm computed for 4 simultaneous float
39  return NaN for x <= 0
40  */
41  _EIGEN_DECLARE_CONST_Packet4f(cephes_SQRTHF, 0.707106781186547524f);
42  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p0, 7.0376836292E-2f);
43  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p1, - 1.1514610310E-1f);
44  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p2, 1.1676998740E-1f);
45  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p3, - 1.2420140846E-1f);
46  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p4, + 1.4249322787E-1f);
47  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p5, - 1.6668057665E-1f);
48  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p6, + 2.0000714765E-1f);
49  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p7, - 2.4999993993E-1f);
50  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p8, + 3.3333331174E-1f);
51  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q1, -2.12194440e-4f);
52  _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q2, 0.693359375f);
53 
54 
55  Packet4i emm0;
56 
57  /* isvalid_mask is 0 if x < 0 or x is NaN. */
58  Packet4ui isvalid_mask = reinterpret_cast<Packet4ui>(vec_cmpge(x, p4f_ZERO));
59  Packet4ui iszero_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(x, p4f_ZERO));
60 
61  x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */
62  emm0 = vec_sr(reinterpret_cast<Packet4i>(x),
63  reinterpret_cast<Packet4ui>(p4i_23));
64 
65  /* keep only the fractional part */
66  x = pand(x, p4f_inv_mant_mask);
67  x = por(x, p4f_half);
68 
69  emm0 = psub(emm0, p4i_0x7f);
70  Packet4f e = padd(vec_ctf(emm0, 0), p4f_1);
71 
72  /* part2:
73  if( x < SQRTHF ) {
74  e -= 1;
75  x = x + x - 1.0;
76  } else { x = x - 1.0; }
77  */
78  Packet4f mask = reinterpret_cast<Packet4f>(vec_cmplt(x, p4f_cephes_SQRTHF));
79  Packet4f tmp = pand(x, mask);
80  x = psub(x, p4f_1);
81  e = psub(e, pand(p4f_1, mask));
82  x = padd(x, tmp);
83 
84  Packet4f x2 = pmul(x,x);
85  Packet4f x3 = pmul(x2,x);
86 
87  Packet4f y, y1, y2;
88  y = pmadd(p4f_cephes_log_p0, x, p4f_cephes_log_p1);
89  y1 = pmadd(p4f_cephes_log_p3, x, p4f_cephes_log_p4);
90  y2 = pmadd(p4f_cephes_log_p6, x, p4f_cephes_log_p7);
91  y = pmadd(y , x, p4f_cephes_log_p2);
92  y1 = pmadd(y1, x, p4f_cephes_log_p5);
93  y2 = pmadd(y2, x, p4f_cephes_log_p8);
94  y = pmadd(y, x3, y1);
95  y = pmadd(y, x3, y2);
96  y = pmul(y, x3);
97 
98  y1 = pmul(e, p4f_cephes_log_q1);
99  tmp = pmul(x2, p4f_half);
100  y = padd(y, y1);
101  x = psub(x, tmp);
102  y2 = pmul(e, p4f_cephes_log_q2);
103  x = padd(x, y);
104  x = padd(x, y2);
105  // negative arg will be NAN, 0 will be -INF
106  x = vec_sel(x, p4f_minus_inf, iszero_mask);
107  x = vec_sel(p4f_minus_nan, x, isvalid_mask);
108  return x;
109 }
110 
111 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
112 Packet4f pexp<Packet4f>(const Packet4f& _x)
113 {
114  Packet4f x = _x;
115  _EIGEN_DECLARE_CONST_Packet4f(1 , 1.0f);
116  _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
117  _EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f);
118  _EIGEN_DECLARE_CONST_Packet4i(23, 23);
119 
120 
121  _EIGEN_DECLARE_CONST_Packet4f(exp_hi, 88.3762626647950f);
122  _EIGEN_DECLARE_CONST_Packet4f(exp_lo, -88.3762626647949f);
123 
124  _EIGEN_DECLARE_CONST_Packet4f(cephes_LOG2EF, 1.44269504088896341f);
125  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C1, 0.693359375f);
126  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C2, -2.12194440e-4f);
127 
128  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p0, 1.9875691500E-4f);
129  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p1, 1.3981999507E-3f);
130  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p2, 8.3334519073E-3f);
131  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p3, 4.1665795894E-2f);
132  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p4, 1.6666665459E-1f);
133  _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p5, 5.0000001201E-1f);
134 
135  Packet4f tmp, fx;
136  Packet4i emm0;
137 
138  // clamp x
139  x = vec_max(vec_min(x, p4f_exp_hi), p4f_exp_lo);
140 
141  /* express exp(x) as exp(g + n*log(2)) */
142  fx = pmadd(x, p4f_cephes_LOG2EF, p4f_half);
143 
144  fx = vec_floor(fx);
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 = vec_cts(fx, 0);
164  emm0 = vec_add(emm0, p4i_0x7f);
165  emm0 = vec_sl(emm0, reinterpret_cast<Packet4ui>(p4i_23));
166 
167  // Altivec's max & min operators just drop silent NaNs. Check NaNs in
168  // inputs and return them unmodified.
169  Packet4ui isnumber_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(_x, _x));
170  return vec_sel(_x, pmax(pmul(y, reinterpret_cast<Packet4f>(emm0)), _x),
171  isnumber_mask);
172 }
173 
174 #ifdef __VSX__
175 // VSX support varies between different compilers and even different
176 // versions of the same compiler. For gcc version >= 4.9.3, we can use
177 // vec_cts to efficiently convert Packet2d to Packet2l. Otherwise, use
178 // a slow version that works with older compilers.
179 static inline Packet2l ConvertToPacket2l(const Packet2d& x) {
180 #if EIGEN_GNUC_AT_LEAST(5, 0) || \
181  (EIGEN_GNUC_AT(4, 9) && __GNUC_PATCHLEVEL__ >= 3)
182  return vec_cts(x, 0); // TODO: check clang version.
183 #else
184  double tmp[2];
185  memcpy(tmp, &x, sizeof(tmp));
186  Packet2l l = { static_cast<long long>(tmp[0]),
187  static_cast<long long>(tmp[1]) };
188  return l;
189 #endif
190 }
191 
192 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
193 Packet2d pexp<Packet2d>(const Packet2d& _x)
194 {
195  Packet2d x = _x;
196 
197  _EIGEN_DECLARE_CONST_Packet2d(1 , 1.0);
198  _EIGEN_DECLARE_CONST_Packet2d(2 , 2.0);
199  _EIGEN_DECLARE_CONST_Packet2d(half, 0.5);
200 
201  _EIGEN_DECLARE_CONST_Packet2d(exp_hi, 709.437);
202  _EIGEN_DECLARE_CONST_Packet2d(exp_lo, -709.436139303);
203 
204  _EIGEN_DECLARE_CONST_Packet2d(cephes_LOG2EF, 1.4426950408889634073599);
205 
206  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p0, 1.26177193074810590878e-4);
207  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p1, 3.02994407707441961300e-2);
208  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p2, 9.99999999999999999910e-1);
209 
210  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q0, 3.00198505138664455042e-6);
211  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q1, 2.52448340349684104192e-3);
212  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q2, 2.27265548208155028766e-1);
213  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q3, 2.00000000000000000009e0);
214 
215  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C1, 0.693145751953125);
216  _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C2, 1.42860682030941723212e-6);
217 
218  Packet2d tmp, fx;
219  Packet2l emm0;
220 
221  // clamp x
222  x = pmax(pmin(x, p2d_exp_hi), p2d_exp_lo);
223  /* express exp(x) as exp(g + n*log(2)) */
224  fx = pmadd(p2d_cephes_LOG2EF, x, p2d_half);
225 
226  fx = vec_floor(fx);
227 
228  tmp = pmul(fx, p2d_cephes_exp_C1);
229  Packet2d z = pmul(fx, p2d_cephes_exp_C2);
230  x = psub(x, tmp);
231  x = psub(x, z);
232 
233  Packet2d x2 = pmul(x,x);
234 
235  Packet2d px = p2d_cephes_exp_p0;
236  px = pmadd(px, x2, p2d_cephes_exp_p1);
237  px = pmadd(px, x2, p2d_cephes_exp_p2);
238  px = pmul (px, x);
239 
240  Packet2d qx = p2d_cephes_exp_q0;
241  qx = pmadd(qx, x2, p2d_cephes_exp_q1);
242  qx = pmadd(qx, x2, p2d_cephes_exp_q2);
243  qx = pmadd(qx, x2, p2d_cephes_exp_q3);
244 
245  x = pdiv(px,psub(qx,px));
246  x = pmadd(p2d_2,x,p2d_1);
247 
248  // build 2^n
249  emm0 = ConvertToPacket2l(fx);
250 
251 #ifdef __POWER8_VECTOR__
252  static const Packet2l p2l_1023 = { 1023, 1023 };
253  static const Packet2ul p2ul_52 = { 52, 52 };
254 
255  emm0 = vec_add(emm0, p2l_1023);
256  emm0 = vec_sl(emm0, p2ul_52);
257 #else
258  // Code is a bit complex for POWER7. There is actually a
259  // vec_xxsldi intrinsic but it is not supported by some gcc versions.
260  // So we shift (52-32) bits and do a word swap with zeros.
261  _EIGEN_DECLARE_CONST_Packet4i(1023, 1023);
262  _EIGEN_DECLARE_CONST_Packet4i(20, 20); // 52 - 32
263 
264  Packet4i emm04i = reinterpret_cast<Packet4i>(emm0);
265  emm04i = vec_add(emm04i, p4i_1023);
266  emm04i = vec_sl(emm04i, reinterpret_cast<Packet4ui>(p4i_20));
267  static const Packet16uc perm = {
268  0x14, 0x15, 0x16, 0x17, 0x00, 0x01, 0x02, 0x03,
269  0x1c, 0x1d, 0x1e, 0x1f, 0x08, 0x09, 0x0a, 0x0b };
270 #ifdef _BIG_ENDIAN
271  emm0 = reinterpret_cast<Packet2l>(vec_perm(p4i_ZERO, emm04i, perm));
272 #else
273  emm0 = reinterpret_cast<Packet2l>(vec_perm(emm04i, p4i_ZERO, perm));
274 #endif
275 
276 #endif
277 
278  // Altivec's max & min operators just drop silent NaNs. Check NaNs in
279  // inputs and return them unmodified.
280  Packet2ul isnumber_mask = reinterpret_cast<Packet2ul>(vec_cmpeq(_x, _x));
281  return vec_sel(_x, pmax(pmul(x, reinterpret_cast<Packet2d>(emm0)), _x),
282  isnumber_mask);
283 }
284 #endif
285 
286 } // end namespace internal
287 
288 } // end namespace Eigen
289 
290 #endif // EIGEN_MATH_FUNCTIONS_ALTIVEC_H
Definition: LDLT.h:16
Definition: Eigen_Colamd.h:54