How to avoid floating point exception when using pow?

Question

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

Are you sure your exception is due to this line? When I try at home your exact example I get 5.71612e-313. Can you give more details about the library you use for "pow", the compiler you use, on what system, etc? Perhaps you should try to upgrade something. Floating-point exceptions are mostly due to divisions by 0. — Benjamin Barrois, Jan 24 '18 at 14:54
In addition to a minimal, complete, verifiable example, what platform are you running on—hardware, compiler, operating system, including versions? `pow(.62, 1504)` is representable in the common 64-bit `double` format. It is subnormal, but even implementations that flush subnormals to zero ought to return a result for this, not trap. — Eric Postpischil, Jan 24 '18 at 14:55
Ok, I use Debian 9 64 bits and I compile with gcc 6.3. The library is c cmaes and here is the problematic line: https://github.com/CMA-ES/c-cmaes/blob/eda8268ee4c8c9fbe4d2489555ae08f8a8c949b5/src/cmaes.c#L893 — AntoineMazuyer, Jan 24 '18 at 14:59
This "problematic line" is containing division. Which can give division by zero. I have no idea how you deduced the problem is in `pow`. — Eugene Sh., Jan 24 '18 at 15:01
If you take the lines you presented in the question, give it the values you presented in the question and compile it with the compiler you state in the question - does it give you the exception? — Eugene Sh., Jan 24 '18 at 15:04
About the edit: there are *two* divisions. I don't understand how this question is getting so many upvotes as it is most likely completely misleading. — Eugene Sh., Jan 24 '18 at 15:05
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/163808/discussion-between-antoinemazuyer-and-eugene-sh). — AntoineMazuyer, Jan 24 '18 at 15:29
Is the code you posted supposed to result in a floating point exception? — Christian Gibbons, Jan 24 '18 at 15:33
No It works whereas it is the same values given by gdb when I run the library from which this line is taken. I have a FPE using the library, but it works like a charm using the code sample — AntoineMazuyer, Jan 24 '18 at 15:38
The snippet you posted is a bit different from the "offending" line in the linked code, which is basically `int hsign = /*some calc with doubles*/ < /*some double value*/ `. — Bob__, Jan 24 '18 at 15:41
i suspect FPE is overflow due to division by small number (maybe 0) and not anyhing to do with pow() troubles. break the line of code into steps to see the true culprit. — chux - Reinstate Monica, Jan 24 '18 at 20:25

Bathsheba · Accepted Answer · 2018-01-26T11:31:41.737

2

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

edited Jan 26 '18 at 11:31

answered Jan 24 '18 at 14:37

Bathsheba

231,907
34
361
483

1

Don't you mean `ln` instead of `log`? – lost_in_the_source Jan 24 '18 at 14:50
@stackptr: Oops. – Bathsheba Jan 24 '18 at 14:51
1

@Bathsheba: A `pow` implementation generally uses logarithms and exponents, but it is more likely to use base two (to match the base of the floating-point format) and specialized code to calculate the things it needs, rather than calling the actual `exp`, `exp2`, `ln`, or `log2`. Those routines would not provide the needed precision for a good `pow`. – Eric Postpischil Jan 24 '18 at 14:58
@EricPostpischil: A good point, as always. I wonder if it's feasible to still proceed as above, check c * i, and if that satisfies some criteria or other, call the standard pow function. – Bathsheba Jan 24 '18 at 15:01
in c language, log is ln (and log10 is log) https://stackoverflow.com/questions/34205208/how-can-you-use-ln-function-in-c-programing – AntoineMazuyer Jan 24 '18 at 17:19
`pow()` as `exp(c * ln(b))` when 1) `b <= 0` 2) is unnessisarially imprecise, – chux - Reinstate Monica Jan 24 '18 at 21:04
Do not implement pow in terms of exp and ln. It will introduce a loss of precision and will likely be slower. – Paul Floyd Jan 26 '18 at 11:29
@PaulFloyd: I might reword this. – Bathsheba Jan 26 '18 at 11:32
@PaulFloyd: Dunnit. Obviously this invalidates some of the previous comments, but the edit history is always viewable. – Bathsheba Jan 26 '18 at 11:35
I would suggest that before anyone tries to do better than their libc version of pow that they take a look at the source https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/e_pow.c and https://github.com/freebsd/freebsd/blob/master/lib/msun/src/e_powf.c. – Paul Floyd Jan 26 '18 at 11:40

How to avoid floating point exception when using pow?

1 Answers1