cfloating-pointnumericalfloating-point-exceptions

How to avoid floating point exception when using pow?


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


Solution

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