[jbay808]

I wouldn't recommend doing this directly -- there's no good polynomial that can represent the exponential function well over a wide range.

Instead, it's better to exploit the definition of IEEE754 as an exponent and a mantissa. You calculate the exponent directly (because e^x = 2^(x/ln2)), and use the Remez algorithm to find a polynomial that fits just the mantissa.

[dietrichepp]

Honestly, reading this comment made me a bit unhappy. Don't mansplain to me the things that I've done. The Remez algorithm requires that you specify the domain. If you tried to use it to approximate the exponential function over its entire domain the algorithm would diverge.

You have to explicitly choose a range [x1,x2] for the Remez algorithm. It minimizes the maximum error over that range.

I'm ordinarily happy to share the work I'm doing but this comment reminds me of the reasons I don't share much with HN.

[jbay808]

I'm writing this because I've done it too, for a performance-critical embedded application, and I think it's an insight worth sharing, for anybody else who comes across this and is trying to do something similar. Not to put you down or one-up you or anything; I'm very sorry to have come across that way. Can I try again?

For a domain like [0,1], the Remez algorithm will do a great job -- because the exponent is barely changing, just the mantissa. For a range like [-10, 10], you won't find any polynomial of reasonable order that has acceptable error.

But, there's a really neat trick that does let the Remez algorithm work over the entire domain of floats, and get excellent performance with just a fourth or fifth order polynomial, or even third order if you're pressed for time, and minimal extra computation.

That is, with just a little simplification, you multiply x by (1/ln(2)), floor it, and stuff the resulting integer into the exponent bits. Then, you take the remainder (what is left over after taking the floor), and make that your input to the polynomial, and stuff the result into the mantissa bits.

There's a little more to it; signs and NaNs need to be handled correctly, subnormal numbers can be treated as a special case with a different Remez polynomial, and a few other corner cases and exact values might be important (e^0, e^1, e^-1). But the result works wonders and it's basically just using Remez but transforming the input to a domain that behaves more like a polynomial, leveraging the magic of IEEE-754.

[fractionalhare]

This is a good explanation of range reduction with respect to floating point. I use a similar (but not quite identical) techniue in the article; I may expand it with more information about what's happening with the exponent and mantissa bits like you did here.

That reminds me, I'm using the domain [-1, 1] for range reduction. But by symmetry of e^x, I can actually use [0, 1]. Thank you!

[dietrichepp]

Depending on the particular architecture you may prefer [-0.5, +0.5].

[ekelsen]

For those interested in more detail / explanation, the paper below does a good job explaining the basic idea of how this works, but with only a 1st order approximation for the mantissa bits. The first order approximation reduces the calculation of exp to just one multiply and one add - which are often combined into just one instruction on modern processors.

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.9.4...

[dietrichepp]

Most of this is covered in the original article.

[6nf]

> Don't mansplain to me

I honestly don't think this has anything to do with your gender or the gender of the person you're replying to.