FEXP is a very important function that is used by many others. The current function is alright. If is very precise when the fraction of the exponent is near zero, but close to 1 it mounts to 1e-8. Which is good. If we would leave Bill Gates with such a small fraction of his money we'd toss him a dime when we head for work. For ZenFP it is more than sufficient (in 32bit), because that is about all the precision it can handle. But things are different for ANS float.
ANS float can handle a much higher precision and there the error can become annoying. So I set out to find a good Remez variant. And I found it in the source code of Go.
But wait a minute: this one can only handle values between 0 and 0.34-something. That won't do.. or will it.
Well, first of all e^-x equals 1/e^x, so we now expanded our range to -0.34 to 0.34. It's something, but not all.
Integer exponents are easy, since e^x equals e * e * e * (x times). So that takes care of all natural numbers. When we combine it with the previous rule: all integer negative numbers as well. Still, it's not the continuous range we need.
Fortunately, this is true as well:
e^(i+f) = e^(i) * e^(f)
Now we're getting somewhere:
e^(integer + fraction) = e^integer * e^fraction
Still, there is a significant gap between 0.34 and 1. How do we fill that one? Simply by applying the same rule! Let's say we gotta calculate e^-3.70:
e^-3.70 = 1 / ( e^3 * e^0.34 * e^0.34 * e^0.02)
So, let's break that one down:
-: check! Simply apply the inverse 1/x;
3: check! That's simply an integer exponent;
0.34: check! That's a constant, so we can calculate that one right now.
And what's left can be channeled through the Remez constant. Now let's throw in Anton Ertl's FPOW function that does e^30 in just eight multiplications and we got something going.
Now we got all the components and we can assemble our function:
Ok, that's done and the results look very good, especially in the "near 1" range.