Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

Implementing the Exponential Function

pseudorandom.com

Post reply on HN
Print view View original

51–60 of 62 posts

Re: Implementing the Exponential Function

#51

Post by fluffything »

benrbray wrote: post #43

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 t…

And also quite useful in a wide range of numerical methods for PDEs like discontinuous Galerkin (and spectral methods in general).

Re: Implementing the Exponential Function

#52

Post by _0ffh »

nwallin wrote: post #24

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. https://godbolt.org/z/KNyoVd

I think you can spare the "and", because if you overflow the exponent field, the result will be wrong either way.

Re: Implementing the Exponential Function

#53

Post by ekelsen »

jbay808 wrote: post #49
dietrichepp wrote: post #38

Earlier quoted context omitted.

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 s…

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…

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

Re: Implementing the Exponential Function

#54

Post by fractionalhare »

jbay808 wrote: post #49
dietrichepp wrote: post #38

Earlier quoted context omitted.

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 s…

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…

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!

Re: Implementing the Exponential Function

#55

Post by fractionalhare »

xelxebar wrote: post #44

Enjoying the article. Looks like there's a math typo at the tail of the "Precision Bound" section: $$ \log\left(\frac{x e}{n}\right) \leq \frac{1}{n}\log(\eta) $$ is not equivalent to $$ \frac{xe}{n} \leq \eta^{-n} $$ The latter should be $$ \frac{xe}{n} \leq \eta^{1/n} $$ and the "$x^{-n}$" in the following paragraph should be changed accordingly. Certainly understandable flub, though. Working with logarithms, it's…

Nice catch, thank you. I'll edit that and credit your comments in the errata :)

Re: Implementing the Exponential Function

#56

Post by defaultcompany »

"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]. "

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?

Re: Implementing the Exponential Function

#57

Post by mota7 »

defaultcompany wrote: post #56

"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]. " 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…

'half' is ambiguous here: They lose 1 bit of expressiveness for the mantissa, and 1 bit in the exponent, giving an efficiency loss of 2 bits from 32 == ~3%.

So there's very little loss in using single-precision floating point, and a lot of gains in smoothly handling larger numbers that arise from addition et al.

edit: I'm half wrong. The bit in mantissa isn't wasted, the only bit that's wasted in the sign bit in the exponent, so the efficiency loss is ~1.5%

Re: Implementing the Exponential Function

#58

Post by widdma »

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.

Re: Implementing the Exponential Function

#59

Post by dietrichepp »

jbay808 wrote: post #49
dietrichepp wrote: post #38

Earlier quoted context omitted.

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 s…

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…

Most of this is covered in the original article.

Re: Implementing the Exponential Function

#60

Post by dietrichepp »

fractionalhare wrote: post #54
jbay808 wrote: post #49

Earlier quoted context omitted.

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…

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!

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