In adherence with IEEE-754, the following
#include <iostream>
#include <cmath>
int main() {
std::cout << std::exp(-1/0.) << "\n";
return 0;
}
yields 0
as the result. The reason being, -1/0.0
gives -INF
, and exp
of that gives 0
.
Now, I need my code to run with -ffast-math
. I understand that this option makes, in particular, infinities non-IEEE-conformant. So, I don't think I can rely on the intermediate -INF
result.
However, the exponential function is (at least, mathematically) very stable for negative arguments: any finite value < -1e+3
already yields 0
as the result. So, I would be completely ok if the compiler just replaced the infinity with a large finite value.
At least in g++-5.4
, this seems to work reliably, but is it prudent to rely on this?
I'm not yet convinced by the comments and answer saying that the original code is unreliable, but here's what I've decided to go with to be sure:
#include <iostream>
#include <cmath>
#include <limits>
int main() {
double a = 0.; // In practice, this is actually a small, nonnegative,
// usually-but-not-always >0 number
double tiny = 1e-50;
double b = std::exp(-1/(a+tiny));
std::cout << b << "\n";
return 0;
}
This is ugly, but it ensures that the divisor is always greater than zero whilst avoiding any conditional.
I considered using epsilon
instead of tiny
, but it doesn't really express the intent any better. (In fact it would skew the results; what I actually want to express is a + tiny == a
for any a
such that exp(-1/a)
is defined.)