I use a C library which uses the pow
function on two double
values
double a = pow(b, c)
At the moment, I have b = 0.62
and c = 1504
, which means that a
should be nearly 0 (3.6e-312).
But I have a floating point exception. How to avoid it and directly return 0? Can we anticipate this case?
I use Debian 9, 64-bit, and I compile with gcc 6.3. The library is ccmaes and here is the problematic line:
https://github.com/CMA-ES/c-cmaes/blob/eda8268ee4c8c9fbe4d2489555ae08f8a8c949b5/src/cmaes.c#L893
I have used gdb so the floating point exception does not come from the division (t->chiN = 2.74)
If I try to reproduce it, with the values when the FPE occurs, I have no problem (compilation option : -fopenmp -O3 -DNDEBUG -fPIC -Wall -Wextra -Wno-long-long -Wconversion -o, like the library)
#include <math.h>
#include <stdio.h>
int main() {
double psxps = 5.6107247793270769;
double cs = 0.37564049253818982;
int gen = 752;
double chiN = 2.7421432615656891;
int N =8;
double foo = sqrt(psxps) / sqrt(1. - pow(1.-cs, 2*gen)) / chiN
< 1.4 + 2./(N+1);
printf("%lf\n",foo);
}
Result: 1.00000000000
Conceptually at least, pow(b, c)
can be considered as being implemented as exp(c * ln(b))
.
So one way of intercepting a potential numerical problem would be to compute the first part of this conceptual pow
yourself, using
double i = ln(b);
If c * i
is sufficiently small (the threshold value will be a function of the floating point scheme used on your platform), then you can proceed with the evaluation of pow
using the C standard library function.
It's probably unwise to finish the job yourself with exp(c * i)
, since the standard pow
function may well have various tricks up its sleeve to attain a result superior in accuracy to exp(c * ln(b))
.