Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.ExtLibs/HeuristicLab.NativeInterpreter/0.1/NativeInterpreter-0.1/lib/vdt/exp.h @ 17071

Last change on this file since 17071 was 17071, checked in by mkommend, 5 years ago

#2958: Merged 16266, 16269, 16274, 16276, 16277, 16285, 16286, 16287, 16289, 16293, 16296, 16297, 16298, 16333, 16334 into stable.

File size: 4.2 KB
Line 
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@na-net.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
33namespace vdt{
34
35namespace details{
36 
37  const double EXP_LIMIT = 708;
38
39  const double PX1exp = 1.26177193074810590878E-4;
40  const double PX2exp = 3.02994407707441961300E-2;
41  const double PX3exp = 9.99999999999999999910E-1;
42  const double QX1exp = 3.00198505138664455042E-6;
43  const double QX2exp = 2.52448340349684104192E-3;
44  const double QX3exp = 2.27265548208155028766E-1;
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.12194440e-4f;
54
55  const float PX1expf = 1.9875691500E-4f;
56  const float PX2expf =1.3981999507E-3f;
57  const float PX3expf =8.3334519073E-3f;
58  const float PX4expf =4.1665795894E-2f;
59  const float PX5expf =1.6666665459E-1f;
60  const float PX6expf =5.0000001201E-1f;
61
62  const float LOG2EF = 1.44269504088896341f;
63
64}
65
66// Exp double precision --------------------------------------------------------
67
68
69/// Exponential Function double precision
70inline 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.93145751953125E-1;
78    x -= px * 1.42860682030941723212E-6;
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
118inline 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
Note: See TracBrowser for help on using the repository browser.