1  /*


2  * exp.h


3  * The basic idea is to exploit Pade polynomials.


4  * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier


5  * moshier@nanet.ornl.gov) as well as actual code.


6  * The Cephes library can be found here: http://www.netlib.org/cephes/


7  *


8  * Created on: Jun 23, 2012


9  * Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente


10  */


11 


12  /*


13  * VDT is free software: you can redistribute it and/or modify


14  * it under the terms of the GNU Lesser Public License as published by


15  * the Free Software Foundation, either version 3 of the License, or


16  * (at your option) any later version.


17  *


18  * This program is distributed in the hope that it will be useful,


19  * but WITHOUT ANY WARRANTY; without even the implied warranty of


20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


21  * GNU Lesser Public License for more details.


22  *


23  * You should have received a copy of the GNU Lesser Public License


24  * along with this program. If not, see <http://www.gnu.org/licenses/>.


25  */


26 


27  #ifndef _VDT_EXP_


28  #define _VDT_EXP_


29 


30  #include "vdtcore_common.h"


31  #include <limits>


32 


33  namespace vdt{


34 


35  namespace details{


36 


37  const double EXP_LIMIT = 708;


38 


39  const double PX1exp = 1.26177193074810590878E4;


40  const double PX2exp = 3.02994407707441961300E2;


41  const double PX3exp = 9.99999999999999999910E1;


42  const double QX1exp = 3.00198505138664455042E6;


43  const double QX2exp = 2.52448340349684104192E3;


44  const double QX3exp = 2.27265548208155028766E1;


45  const double QX4exp = 2.00000000000000000009E0;


46 


47  const double LOG2E = 1.4426950408889634073599; // 1/log(2)


48 


49  const float MAXLOGF = 88.72283905206835f;


50  const float MINLOGF = 88.f;


51 


52  const float C1F = 0.693359375f;


53  const float C2F = 2.12194440e4f;


54 


55  const float PX1expf = 1.9875691500E4f;


56  const float PX2expf =1.3981999507E3f;


57  const float PX3expf =8.3334519073E3f;


58  const float PX4expf =4.1665795894E2f;


59  const float PX5expf =1.6666665459E1f;


60  const float PX6expf =5.0000001201E1f;


61 


62  const float LOG2EF = 1.44269504088896341f;


63 


64  }


65 


66  // Exp double precision 


67 


68 


69  /// Exponential Function double precision


70  inline double fast_exp(double initial_x){


71 


72  double x = initial_x;


73  double px=details::fpfloor(details::LOG2E * x +0.5);


74 


75  const int32_t n = int32_t(px);


76 


77  x = px * 6.93145751953125E1;


78  x = px * 1.42860682030941723212E6;


79 


80  const double xx = x * x;


81 


82  // px = x * P(x**2).


83  px = details::PX1exp;


84  px *= xx;


85  px += details::PX2exp;


86  px *= xx;


87  px += details::PX3exp;


88  px *= x;


89 


90  // Evaluate Q(x**2).


91  double qx = details::QX1exp;


92  qx *= xx;


93  qx += details::QX2exp;


94  qx *= xx;


95  qx += details::QX3exp;


96  qx *= xx;


97  qx += details::QX4exp;


98 


99  // e**x = 1 + 2x P(x**2)/( Q(x**2)  P(x**2) )


100  x = px / (qx  px);


101  x = 1.0 + 2.0 * x;


102 


103  // Build 2^n in double.


104  x *= details::uint642dp(( ((uint64_t)n) +1023)<<52);


105 


106  if (initial_x > details::EXP_LIMIT)


107  x = std::numeric_limits<double>::infinity();


108  if (initial_x < details::EXP_LIMIT)


109  x = 0.;


110 


111  return x;


112 


113  }


114 


115  // Exp single precision 


116 


117  /// Exponential Function single precision


118  inline float fast_expf(float initial_x) {


119 


120  float x = initial_x;


121 


122  float z = details::fpfloor( details::LOG2EF * x +0.5f ); /* floor() truncates toward infinity. */


123 


124  x = z * details::C1F;


125  x = z * details::C2F;


126  const int32_t n = int32_t ( z );


127 


128  const float x2 = x * x;


129 


130  z = x*details::PX1expf;


131  z += details::PX2expf;


132  z *= x;


133  z += details::PX3expf;


134  z *= x;


135  z += details::PX4expf;


136  z *= x;


137  z += details::PX5expf;


138  z *= x;


139  z += details::PX6expf;


140  z *= x2;


141  z += x + 1.0f;


142 


143  /* multiply by power of 2 */


144  z *= details::uint322sp((n+0x7f)<<23);


145 


146  if (initial_x > details::MAXLOGF) z=std::numeric_limits<float>::infinity();


147  if (initial_x < details::MINLOGF) z=0.f;


148 


149  return z;


150 


151  }


152 


153  //


154 


155  } // end namespace vdt


156 


157  #endif

