I made a fast approximation of the exponential function (e^x) based on a 2nd-degree polynomial. However, the program doesn't give reasonable numbers if x<0.0 .
The code is like this:
// Compile with: gcc -std=c99 fastexp.c -o fastexp
#include <stdio.h>
#include <math.h>
#define ALPHA 0.50617024f
#define ALPHAc 0.49382976f
#define INV_LN2 1.4426950408889634f
float fast2pxm1_m(float x) { // 0.0 <= x <= 1.0
return x*(ALPHA+ALPHAc*x);
}
float fast2px(float x) {
float m, mp;
float e;
m = modff(x, &e);
mp = fast2pxm1_m(m) + 1.0f;
return ldexpf(mp, e);
}
float fastexp(float x) { // exportable
return fast2px(x*INV_LN2);
}
int main(int argc, char *argv[]) {
float x;
for (int i = -10; i <= 10; i++) {
x= (float) i/10.0*2;
printf("%.8f => %.8f\n",x,fastexp(x));
}
return 0;
}
The results of the function fastexp are not good for x<0.0: its evolution along x is erroneous, with local minima in several places (x). See the results of the program:
-2.00000000 => 0.23474069
-1.79999995 => 0.21845233
-1.60000002 => 0.22272082
-1.39999998 => 0.24754614
-1.20000005 => 0.44696173
-1.00000000 => 0.43635058
-0.80000001 => 0.46685311
-0.60000002 => 0.93187356
-0.40000001 => 0.87235498
-0.20000000 => 0.89506382
0.00000000 => 1.00000000
0.20000000 => 1.18716359
0.40000001 => 1.45655441
0.60000002 => 1.80817270
0.80000001 => 2.17952919
1.00000000 => 2.64171839
1.20000005 => 3.26836252
1.39999998 => 4.04080629
1.60000002 => 4.81200027
1.79999995 => 5.91210270
2.00000000 => 7.34111547
EDIT:
The constants below give a better approximation:
#define ALPHA 0.6563658616549711f
#define ALPHAc 0.34363413834502887f
The polynomial you're using in fast2pxm1_m approximates the function pow(2, x) - 1, but it does this well only in the range [0, 1], as evidenced by the comment // 0.0 <= x <= 1.0.
The value you're passing to that function is the return value from modff which has the same sign as the input (and therefore is in the range (-1, 1)), so the result is inaccurate for negative inputs.
One way to fix this is to modify the output of modff so that the fractional part is always non-negative:
m = modff(x, &e);
if (m < 0.0f) {
m += 1.0f; // fractional part
e -= 1.0f; // integral part
}
This gives more reasonable outputs:
x fastexp(x) exp(x)
----------------------------------------
-2.00000000 => 0.13306235 0.13533528
-1.79999995 => 0.16054048 0.16529890
-1.60000002 => 0.19829696 0.20189651
-1.39999998 => 0.24633193 0.24659697
-1.20000005 => 0.29292828 0.30119420
-1.00000000 => 0.35886729 0.36787944
-0.80000001 => 0.44536310 0.44932896
-0.60000002 => 0.53846931 0.54881162
-0.40000001 => 0.65119916 0.67032004
-0.20000000 => 0.80504274 0.81873075
0.00000000 => 1.00000000 1.00000000
0.20000000 => 1.18716359 1.22140276
0.40000001 => 1.45655441 1.49182471
0.60000002 => 1.80817270 1.82211884
0.80000001 => 2.17952919 2.22554096
1.00000000 => 2.64171839 2.71828183
1.20000005 => 3.26836252 3.32011708
1.39999998 => 4.04080629 4.05519987
1.60000002 => 4.81200027 4.95303254
1.79999995 => 5.91210270 6.04964718
2.00000000 => 7.34111547 7.38905610