3 ms·
tl;dr: expm1 does work. Something funny is going on if it doesn't work for you. ~~~ How this works out in practice is that once you try your hands on some ser
by mehrdadn 7y ago
tl;dr: expm1 does work. Something funny is going on if it doesn't work for you.
~~~
How this works out in practice is that once you try your hands on some serious floating-point, you quickly learn that an expression like x == 0 (or even like x <= 0) is a strong code smell. It's an instant red flag that you need to drop what you're doing and think things through extremely carefully first. So you'd stop right there and try to get some intuition for what's going on.
Now the intuition you want here is that floating point is logarithmically spaced around 0 rather than around 1, so if you do exp(x) for small x, it won't have adequate resolution to return anything but 1, because the next-highest jump is too large to be accurate. So when you subtract 1, you get zero.
But... that's precisely why expm1 exists: it's to evaluate exp(x) - 1 accurately. And in fact, it should work. And it does work for me. Why doesn't it work here? It sounds like either a bad implementation or the wrong compiler flags. We really need to know the compiler flags that are being used.
For what it's worth: floating-point can be a lot worse than this, where there's no simple red flag or library function like this that you can recognize. In my experience, a rule of thumb (somewhat of a yellow flag) is that, if you start doing nonlinear operations (even division), you gotta roll up your sleeves and get ready for some fun...
- acqq 7y ago> once you try your hands on some serious floating-point, you quickly learn that an expression like x == 0 (or even like x <= 0) is a strong code smell. Yes, specifically here is that comparison that surely should not be written as such: template <typename T> T e(T x) { return x == 0. ? 1. : (exp(x) - 1.)/x; } and also: https://en.cppreference.com/w/cpp/numeric/math/expm1 https://en.cppreference.com/w/cpp/numeric/math/expm1
- gumby 7y ago> x == 0 (or even like x <= 0) is a strong code smell. Careful: 0.0f is defined by the standard to be 0x0. So such comparisons can be quite reasonable. However trying for 0 as the result of a computation, as is used in this example is indeed problematic. But there is plenty of reasonable code that initializes a float or double variable to 0.0 and can expect to be able to inspect it and find its bit pattern to be 0x0 before that variable is manipulated. After that, all bets are off.
- gjm11 7y agoThe point isn't just that you can't tell exactly when something is zero, although indeed you often can't. It's that if you have some calculation that needs to be special-cased at zero (e.g., because you're dividing by something and it might explode) then it probably actually needs to be special-cased near zero (because near zero it might overflow or run into other numerical problems), and it's probably worth seeing if there's another way of organizing the calculation that doesn't misbehave for small values in the first place.
- deleted 7y ago[deleted]