[dietrichepp]

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.

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.

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.

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.

[dvt]

This is way cool! Mad respect to you audio engineers out there. I tried implementing FFT last year in one of my open-source projects[1] (key word: tried).

[1] https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de...

[nayuki]

Glancing at your code diff, I think you're missing the bit-reversed permutation. Either way, feel free to steal from: https://www.nayuki.io/page/free-small-fft-in-multiple-langua...

[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.

[exmadscientist]

Remez is the most common implementation of Minimax approximation algorithms, so it's probably coming when the author makes it to section 6.

[fractionalhare]

Yes, this is the plan :)

[nitrogen]

I used something similar to implement an approximate square root function in a graphical DSP platform that only has addition and multiplication. A 4th-order polynomial gives satisfactory accuracy from -20dBFS to 0 dBFS, which was the region I cared about. It's hard to do better near the origin (-inf dB) though, because the slope of the square root approaches infinity.