A couple other ways to compute the exponential function: The Intel 8087 coprocessor used CORDIC. The chip contained a ROM with the values of log2(1+2^-n) for various n. These constants allowed the exponential to be rapidly computed with shifts and adds.<p>The Sinclair Scientific calculator was notable for cramming transcendental functions into a chip designed as a 4-function calculator, so they took some severe shortcuts. Exponentials were based on computing 0.99^n through repeated multiplication. Multiplying by 0.99 is super-cheap in BCD since you shift by two digits and subtract. To compute 10^x, it computes .99^(-229.15*x). Slow and inaccurate, but easy to implement.
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).<p>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!<p>[1] (p35) <a href="http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf" rel="nofollow">http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradi...</a><p>[2] <a href="http://blog.mrtz.org/2013/09/07/the-zen-of-gradient-descent.html" rel="nofollow">http://blog.mrtz.org/2013/09/07/the-zen-of-gradient-descent....</a><p>[3] <a href="https://benrbray.com/static/notes/chebychev-polynomials_dec16.pdf" rel="nofollow">https://benrbray.com/static/notes/chebychev-polynomials_dec1...</a><p>[4] <a href="https://benrbray.com/static/notes/minimax-approx_nov16.pdf" rel="nofollow">https://benrbray.com/static/notes/minimax-approx_nov16.pdf</a>
This is an awesome article! One of my favorite SO answers I remember researching dealt with the Padé approximation of tanh[1] (which I found was significantly better than the Taylor expansion), but the caveat was that it only worked within a very narrow neighborhood.<p>I will say that the article didn't really touch on techniques that minimize IEEE 754 subnormal[2] performance impact, which is a very interesting problem in itself. A lot of real-life implementations of something like e^x will have various clever ways to avoid subnormals.<p>[1] <a href="https://stackoverflow.com/a/6118100/243613" rel="nofollow">https://stackoverflow.com/a/6118100/243613</a><p>[2] <a href="https://mathworld.wolfram.com/SubnormalNumber.html" rel="nofollow">https://mathworld.wolfram.com/SubnormalNumber.html</a>
I recently implemented the exponential function for a sound synthesizer toolkit I’m working on. The method I used was not mentioned in the article, so I’ll explain it here.<p>I used the Remez algorithm, modified to minimize equivalent input error. This algorithm lets you find a polynomial with the smallest maximum error. You start with a set of X coordinates, and create a polynomial which oscillates up and down around the target function, so p(x0) = f(x0) + E, p(x1) = f(x1) - E, p(x2) = f(x2) + E, etc.<p>You then iterate by replacing the set of x values with the x values which have maximum error. This will converge to a polynomial with smallest maximum error.<p>From my experiments, you can use a fairly low order approximation for musical applications. Used to calculate frequencies of musical notes, 2nd order is within 3.0 cents, and 3rd order is within 0.13 cents.
Need to be careful about confusing smooth (infinitely differentiable) and analytic (= to Taylor series) functions.<p>> Any function with this property can be uniquely represented as a Taylor series expansion “centered” at a<p>It's fairly easy to construct tame looking functions where this isn't true.
[Aside on Taylor Series] The section on recurrence relations is the basis of a little-known ODE solving technique sometimes known as the Taylor Series Method (occasionally Differential Transform). It allows an "extension" of Euler's Method to arbitrary order (using arbitrary precision arithmetic), by _calculating_ the higher derivatives rather than approximating them as in RK4 and friends.<p>Here is an open access link to just one of many papers on the subject: <a href="https://projecteuclid.org/euclid.em/1120145574" rel="nofollow">https://projecteuclid.org/euclid.em/1120145574</a><p>It includes a link to the actual software they used.<p>[EDIT] Obviously applying the technique above to dx/dt = x reproduces the Taylor Series in the article!
Really interesting article!<p>Some other resources you may find interesting<p><a href="https://fpbench.org/index.html" rel="nofollow">https://fpbench.org/index.html</a> - I think I saw a talk at FPTalks that felt related to your topic "Creating correctly rounded math libraries for real number approximations"<p><a href="https://fpbench.org/community.html" rel="nofollow">https://fpbench.org/community.html</a><p><a href="https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804" rel="nofollow">https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804</a> - Correctly Rounded libm<p><a href="https://www.mpfr.org/" rel="nofollow">https://www.mpfr.org/</a> - GNU Multi precision floating point<p><a href="http://sollya.gforge.inria.fr/" rel="nofollow">http://sollya.gforge.inria.fr/</a> - A tool for floating point code development
A collaborator and I recently had fun code golfing integer implementations of floor(106*log2(x)) and floor(log2(x!)):<p><a href="https://github.com/sipa/minisketch/blob/master/src/false_positives.h#L20" rel="nofollow">https://github.com/sipa/minisketch/blob/master/src/false_pos...</a>
Excellent article; it's very detailed in the mathematical explanations and code.<p>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. <a href="https://www.nayuki.io/page/approximating-eulers-number-correctly" rel="nofollow">https://www.nayuki.io/page/approximating-eulers-number-corre...</a>
"We can see the values are relatively densely clustered near 0 and increasingly sparse the further you move away from the origin. As another example, half of all 32-bit floating point numbers reside in the real interval [-1,1]. "<p>I knew this in general but this specific statistic was eye-opening to me. Systems like OpenGL normalize coordinates so they only use values in this [-1,1] range. That means they effectively lose half the expressive power of the 32 bits. Are there other machine representations of real numbers which would be better for that use case? i.e. restricted to that range and evenly distributed?
Ugh I had to implement all the transcendental functions a while ago. It was miserable, especially dealing with intermediate rounding, etc<p>I don’t recommend it.<p>That said I liked that there’s a specific exp-1 function (otherwise you drop most precision).
Enjoying the article. Looks like there's a math typo at the tail of the "Precision Bound" section:<p><pre><code> $$ \log\left(\frac{x e}{n}\right) \leq \frac{1}{n}\log(\eta) $$
</code></pre>
is not equivalent to<p><pre><code> $$ \frac{xe}{n} \leq \eta^{-n} $$
</code></pre>
The latter should be<p><pre><code> $$ \frac{xe}{n} \leq \eta^{1/n} $$
</code></pre>
and the "$x^{-n}$" in the following paragraph should be changed accordingly.<p>Certainly understandable flub, though. Working with logarithms, it's easy to mix up minus sign vs. 1/n conversions.
There's 2 cool things related to floating point error I learned recently.<p>First, `x + x + x + x + x == 5<i>x`. 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.<p>Second, SMT solvers can actually prove stuff like this relatively easy! In Z3py, one can prove this in 4 lines of code.<p>from z3 import </i>
x = FP('x', FPSort(8, 24))
set_default_rounding_mode(RNE())
solve(5*x != x + x+ x+ x + x)
IMHO it's a bit cheeky to use exp2 from the standard library. It's likely got all the same math exp does. Fortunately, if you know that your exponent is an integer, there's a constant time version. It's 5 instructions: add, and, shift, mov, ret.<p><a href="https://godbolt.org/z/KNyoVd" rel="nofollow">https://godbolt.org/z/KNyoVd</a>
Nice article, thanks. For anyone interested in the topic, I highly recommend Trefethen's Approximation Theory and Approximation Practice. It's approachable, intuitive, and fun while still covering a lot of technical detail.
Few miscellaneous meta nitpicks:<p>* The article looks beautiful, but am I the only one that hates left-centered content? Wide monitors are the norm these days, and reading left centered content literally hurts my neck (need to resize manually to center the content, approximately.<p>* Why would someone make the effort to write such detailed analysis and conceal their identity? I just don't get what's the objective of not providing any sort of provenance... just a gripe of mine.