![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2007 Julien Pommier 00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // Copyright (C) 2016 Konstantinos Margaritis <markos@freevec.org> 00007 // 00008 // This Source Code Form is subject to the terms of the Mozilla 00009 // Public License v. 2.0. If a copy of the MPL was not distributed 00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00011 00012 /* The sin, cos, exp, and log functions of this file come from 00013 * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/ 00014 */ 00015 00016 #ifndef EIGEN_MATH_FUNCTIONS_ALTIVEC_H 00017 #define EIGEN_MATH_FUNCTIONS_ALTIVEC_H 00018 00019 namespace Eigen { 00020 00021 namespace internal { 00022 00023 static _EIGEN_DECLARE_CONST_Packet4f(1 , 1.0f); 00024 static _EIGEN_DECLARE_CONST_Packet4f(half, 0.5f); 00025 static _EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f); 00026 static _EIGEN_DECLARE_CONST_Packet4i(23, 23); 00027 00028 static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(inv_mant_mask, ~0x7f800000); 00029 00030 /* the smallest non denormalized float number */ 00031 static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos, 0x00800000); 00032 static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf, 0xff800000); // -1.f/0.f 00033 static _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_nan, 0xffffffff); 00034 00035 /* natural logarithm computed for 4 simultaneous float 00036 return NaN for x <= 0 00037 */ 00038 static _EIGEN_DECLARE_CONST_Packet4f(cephes_SQRTHF, 0.707106781186547524f); 00039 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p0, 7.0376836292E-2f); 00040 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p1, - 1.1514610310E-1f); 00041 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p2, 1.1676998740E-1f); 00042 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p3, - 1.2420140846E-1f); 00043 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p4, + 1.4249322787E-1f); 00044 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p5, - 1.6668057665E-1f); 00045 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p6, + 2.0000714765E-1f); 00046 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p7, - 2.4999993993E-1f); 00047 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_p8, + 3.3333331174E-1f); 00048 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q1, -2.12194440e-4f); 00049 static _EIGEN_DECLARE_CONST_Packet4f(cephes_log_q2, 0.693359375f); 00050 00051 static _EIGEN_DECLARE_CONST_Packet4f(exp_hi, 88.3762626647950f); 00052 static _EIGEN_DECLARE_CONST_Packet4f(exp_lo, -88.3762626647949f); 00053 00054 static _EIGEN_DECLARE_CONST_Packet4f(cephes_LOG2EF, 1.44269504088896341f); 00055 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C1, 0.693359375f); 00056 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_C2, -2.12194440e-4f); 00057 00058 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p0, 1.9875691500E-4f); 00059 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p1, 1.3981999507E-3f); 00060 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p2, 8.3334519073E-3f); 00061 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p3, 4.1665795894E-2f); 00062 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p4, 1.6666665459E-1f); 00063 static _EIGEN_DECLARE_CONST_Packet4f(cephes_exp_p5, 5.0000001201E-1f); 00064 00065 #ifdef __VSX__ 00066 static _EIGEN_DECLARE_CONST_Packet2d(1 , 1.0); 00067 static _EIGEN_DECLARE_CONST_Packet2d(2 , 2.0); 00068 static _EIGEN_DECLARE_CONST_Packet2d(half, 0.5); 00069 00070 static _EIGEN_DECLARE_CONST_Packet2d(exp_hi, 709.437); 00071 static _EIGEN_DECLARE_CONST_Packet2d(exp_lo, -709.436139303); 00072 00073 static _EIGEN_DECLARE_CONST_Packet2d(cephes_LOG2EF, 1.4426950408889634073599); 00074 00075 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p0, 1.26177193074810590878e-4); 00076 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p1, 3.02994407707441961300e-2); 00077 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_p2, 9.99999999999999999910e-1); 00078 00079 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q0, 3.00198505138664455042e-6); 00080 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q1, 2.52448340349684104192e-3); 00081 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q2, 2.27265548208155028766e-1); 00082 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q3, 2.00000000000000000009e0); 00083 00084 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C1, 0.693145751953125); 00085 static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C2, 1.42860682030941723212e-6); 00086 00087 #ifdef __POWER8_VECTOR__ 00088 static Packet2l p2l_1023 = { 1023, 1023 }; 00089 static Packet2ul p2ul_52 = { 52, 52 }; 00090 #endif 00091 00092 #endif 00093 00094 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00095 Packet4f plog<Packet4f>(const Packet4f& _x) 00096 { 00097 Packet4f x = _x; 00098 00099 Packet4i emm0; 00100 00101 /* isvalid_mask is 0 if x < 0 or x is NaN. */ 00102 Packet4ui isvalid_mask = reinterpret_cast<Packet4ui>(vec_cmpge(x, p4f_ZERO)); 00103 Packet4ui iszero_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(x, p4f_ZERO)); 00104 00105 x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */ 00106 emm0 = vec_sr(reinterpret_cast<Packet4i>(x), 00107 reinterpret_cast<Packet4ui>(p4i_23)); 00108 00109 /* keep only the fractional part */ 00110 x = pand(x, p4f_inv_mant_mask); 00111 x = por(x, p4f_half); 00112 00113 emm0 = psub(emm0, p4i_0x7f); 00114 Packet4f e = padd(vec_ctf(emm0, 0), p4f_1); 00115 00116 /* part2: 00117 if( x < SQRTHF ) { 00118 e -= 1; 00119 x = x + x - 1.0; 00120 } else { x = x - 1.0; } 00121 */ 00122 Packet4f mask = reinterpret_cast<Packet4f>(vec_cmplt(x, p4f_cephes_SQRTHF)); 00123 Packet4f tmp = pand(x, mask); 00124 x = psub(x, p4f_1); 00125 e = psub(e, pand(p4f_1, mask)); 00126 x = padd(x, tmp); 00127 00128 Packet4f x2 = pmul(x,x); 00129 Packet4f x3 = pmul(x2,x); 00130 00131 Packet4f y, y1, y2; 00132 y = pmadd(p4f_cephes_log_p0, x, p4f_cephes_log_p1); 00133 y1 = pmadd(p4f_cephes_log_p3, x, p4f_cephes_log_p4); 00134 y2 = pmadd(p4f_cephes_log_p6, x, p4f_cephes_log_p7); 00135 y = pmadd(y , x, p4f_cephes_log_p2); 00136 y1 = pmadd(y1, x, p4f_cephes_log_p5); 00137 y2 = pmadd(y2, x, p4f_cephes_log_p8); 00138 y = pmadd(y, x3, y1); 00139 y = pmadd(y, x3, y2); 00140 y = pmul(y, x3); 00141 00142 y1 = pmul(e, p4f_cephes_log_q1); 00143 tmp = pmul(x2, p4f_half); 00144 y = padd(y, y1); 00145 x = psub(x, tmp); 00146 y2 = pmul(e, p4f_cephes_log_q2); 00147 x = padd(x, y); 00148 x = padd(x, y2); 00149 // negative arg will be NAN, 0 will be -INF 00150 x = vec_sel(x, p4f_minus_inf, iszero_mask); 00151 x = vec_sel(p4f_minus_nan, x, isvalid_mask); 00152 return x; 00153 } 00154 00155 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00156 Packet4f pexp<Packet4f>(const Packet4f& _x) 00157 { 00158 Packet4f x = _x; 00159 00160 Packet4f tmp, fx; 00161 Packet4i emm0; 00162 00163 // clamp x 00164 x = pmax(pmin(x, p4f_exp_hi), p4f_exp_lo); 00165 00166 // express exp(x) as exp(g + n*log(2)) 00167 fx = pmadd(x, p4f_cephes_LOG2EF, p4f_half); 00168 00169 fx = pfloor(fx); 00170 00171 tmp = pmul(fx, p4f_cephes_exp_C1); 00172 Packet4f z = pmul(fx, p4f_cephes_exp_C2); 00173 x = psub(x, tmp); 00174 x = psub(x, z); 00175 00176 z = pmul(x,x); 00177 00178 Packet4f y = p4f_cephes_exp_p0; 00179 y = pmadd(y, x, p4f_cephes_exp_p1); 00180 y = pmadd(y, x, p4f_cephes_exp_p2); 00181 y = pmadd(y, x, p4f_cephes_exp_p3); 00182 y = pmadd(y, x, p4f_cephes_exp_p4); 00183 y = pmadd(y, x, p4f_cephes_exp_p5); 00184 y = pmadd(y, z, x); 00185 y = padd(y, p4f_1); 00186 00187 // build 2^n 00188 emm0 = vec_cts(fx, 0); 00189 emm0 = vec_add(emm0, p4i_0x7f); 00190 emm0 = vec_sl(emm0, reinterpret_cast<Packet4ui>(p4i_23)); 00191 00192 // Altivec's max & min operators just drop silent NaNs. Check NaNs in 00193 // inputs and return them unmodified. 00194 Packet4ui isnumber_mask = reinterpret_cast<Packet4ui>(vec_cmpeq(_x, _x)); 00195 return vec_sel(_x, pmax(pmul(y, reinterpret_cast<Packet4f>(emm0)), _x), 00196 isnumber_mask); 00197 } 00198 00199 #ifndef EIGEN_COMP_CLANG 00200 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00201 Packet4f prsqrt<Packet4f>(const Packet4f& x) 00202 { 00203 return vec_rsqrt(x); 00204 } 00205 #endif 00206 00207 #ifdef __VSX__ 00208 #ifndef EIGEN_COMP_CLANG 00209 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00210 Packet2d prsqrt<Packet2d>(const Packet2d& x) 00211 { 00212 return vec_rsqrt(x); 00213 } 00214 #endif 00215 00216 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00217 Packet4f psqrt<Packet4f>(const Packet4f& x) 00218 { 00219 return vec_sqrt(x); 00220 } 00221 00222 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00223 Packet2d psqrt<Packet2d>(const Packet2d& x) 00224 { 00225 return vec_sqrt(x); 00226 } 00227 00228 // VSX support varies between different compilers and even different 00229 // versions of the same compiler. For gcc version >= 4.9.3, we can use 00230 // vec_cts to efficiently convert Packet2d to Packet2l. Otherwise, use 00231 // a slow version that works with older compilers. 00232 // Update: apparently vec_cts/vec_ctf intrinsics for 64-bit doubles 00233 // are buggy, https://gcc.gnu.org/bugzilla/show_bug.cgi?id=70963 00234 static inline Packet2l ConvertToPacket2l(const Packet2d& x) { 00235 #if EIGEN_GNUC_AT_LEAST(5, 4) || \ 00236 (EIGEN_GNUC_AT(6, 1) && __GNUC_PATCHLEVEL__ >= 1) 00237 return vec_cts(x, 0); // TODO: check clang version. 00238 #else 00239 double tmp[2]; 00240 memcpy(tmp, &x, sizeof(tmp)); 00241 Packet2l l = { static_cast<long long>(tmp[0]), 00242 static_cast<long long>(tmp[1]) }; 00243 return l; 00244 #endif 00245 } 00246 00247 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 00248 Packet2d pexp<Packet2d>(const Packet2d& _x) 00249 { 00250 Packet2d x = _x; 00251 00252 Packet2d tmp, fx; 00253 Packet2l emm0; 00254 00255 // clamp x 00256 x = pmax(pmin(x, p2d_exp_hi), p2d_exp_lo); 00257 00258 /* express exp(x) as exp(g + n*log(2)) */ 00259 fx = pmadd(x, p2d_cephes_LOG2EF, p2d_half); 00260 00261 fx = pfloor(fx); 00262 00263 tmp = pmul(fx, p2d_cephes_exp_C1); 00264 Packet2d z = pmul(fx, p2d_cephes_exp_C2); 00265 x = psub(x, tmp); 00266 x = psub(x, z); 00267 00268 Packet2d x2 = pmul(x,x); 00269 00270 Packet2d px = p2d_cephes_exp_p0; 00271 px = pmadd(px, x2, p2d_cephes_exp_p1); 00272 px = pmadd(px, x2, p2d_cephes_exp_p2); 00273 px = pmul (px, x); 00274 00275 Packet2d qx = p2d_cephes_exp_q0; 00276 qx = pmadd(qx, x2, p2d_cephes_exp_q1); 00277 qx = pmadd(qx, x2, p2d_cephes_exp_q2); 00278 qx = pmadd(qx, x2, p2d_cephes_exp_q3); 00279 00280 x = pdiv(px,psub(qx,px)); 00281 x = pmadd(p2d_2,x,p2d_1); 00282 00283 // build 2^n 00284 emm0 = ConvertToPacket2l(fx); 00285 00286 #ifdef __POWER8_VECTOR__ 00287 emm0 = vec_add(emm0, p2l_1023); 00288 emm0 = vec_sl(emm0, p2ul_52); 00289 #else 00290 // Code is a bit complex for POWER7. There is actually a 00291 // vec_xxsldi intrinsic but it is not supported by some gcc versions. 00292 // So we shift (52-32) bits and do a word swap with zeros. 00293 _EIGEN_DECLARE_CONST_Packet4i(1023, 1023); 00294 _EIGEN_DECLARE_CONST_Packet4i(20, 20); // 52 - 32 00295 00296 Packet4i emm04i = reinterpret_cast<Packet4i>(emm0); 00297 emm04i = vec_add(emm04i, p4i_1023); 00298 emm04i = vec_sl(emm04i, reinterpret_cast<Packet4ui>(p4i_20)); 00299 static const Packet16uc perm = { 00300 0x14, 0x15, 0x16, 0x17, 0x00, 0x01, 0x02, 0x03, 00301 0x1c, 0x1d, 0x1e, 0x1f, 0x08, 0x09, 0x0a, 0x0b }; 00302 #ifdef _BIG_ENDIAN 00303 emm0 = reinterpret_cast<Packet2l>(vec_perm(p4i_ZERO, emm04i, perm)); 00304 #else 00305 emm0 = reinterpret_cast<Packet2l>(vec_perm(emm04i, p4i_ZERO, perm)); 00306 #endif 00307 00308 #endif 00309 00310 // Altivec's max & min operators just drop silent NaNs. Check NaNs in 00311 // inputs and return them unmodified. 00312 Packet2ul isnumber_mask = reinterpret_cast<Packet2ul>(vec_cmpeq(_x, _x)); 00313 return vec_sel(_x, pmax(pmul(x, reinterpret_cast<Packet2d>(emm0)), _x), 00314 isnumber_mask); 00315 } 00316 #endif 00317 00318 } // end namespace internal 00319 00320 } // end namespace Eigen 00321 00322 #endif // EIGEN_MATH_FUNCTIONS_ALTIVEC_H