Eigen  3.3.3
MathFunctions.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends