Implementing the Exponential Function

I believe that they're all part of the IEEE 754 spec, but some of these functions are just recommended.

x87 has pretty much all of them, including two remainder implementations \o/

Great post! The author touches on Chebychev polynomials, which are the basis for quite a few numerical analysis tricks, including conjugate gradient [1], accelerated gradient descent [2], and Chebychev semi-iterative methods (which find the best combination of past iterates to use during an optimization procedure; sadly I can't find a good reference).

There are a number of facts/folklore about Chebychev polynomials that seem to go stated without proof in a lot of papers, so a few years ago I wrote up some notes [3,4] with the most common claims. Maybe someone will find them useful!

[1] (p35) http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradi...

[2] http://blog.mrtz.org/2013/09/07/the-zen-of-gradient-descent....

[3] https://benrbray.com/static/notes/chebychev-polynomials_dec1...

[4] https://benrbray.com/static/notes/minimax-approx_nov16.pdf

Excellent article; it's very detailed in the mathematical explanations and code.

I was interested in the exponential function too, but approached it from the standpoint of arbitrary-precision arithmetic. My explanation is an order of magnitude shorter, but of course my code is a lot slower than a typical floating-point algorithm. I offer a way to compute correctly rounded exponentials without resorting to a big math package like Mathematica/Maple/etc. https://www.nayuki.io/page/approximating-eulers-number-corre...

There's 2 cool things related to floating point error I learned recently.

First, `x + x + x + x + x == 5x`. This is true for values up to 5, but is not true for 6. Proving this is a fun exercise, but is a little bit painful and has a lot of cases.

Second, SMT solvers can actually prove stuff like this relatively easy! In Z3py, one can prove this in 4 lines of code.

from z3 import x = FP('x', FPSort(8, 24)) set_default_rounding_mode(RNE()) solve(5*x != x + x+ x+ x + x)

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.