使用 pow 时如何避免浮点异常?
How to avoid floating point exception when using pow?
我使用一个 C 库,它在两个 double
值上使用 pow
函数
double a = pow(b, c)
目前,我有b = 0.62
和c = 1504
,这意味着a
应该接近0(3.6e-312)。
但是我有一个浮点异常。如何避免直接return0?我们可以预见这种情况吗?
我使用 Debian 9,64 位,我用 gcc 6.3 编译。库是 ccmaes,这里是有问题的行:
https://github.com/CMA-ES/c-cmaes/blob/eda8268ee4c8c9fbe4d2489555ae08f8a8c949b5/src/cmaes.c#L893
我用过gdb所以浮点异常不是来自除法(t->chiN = 2.74)
如果我尝试重现它,使用 FPE 发生时的值,我没有问题(编译选项:-fopenmp -O3 -DNDEBUG -fPIC -Wall -Wextra -Wno-long-long -Wconversion -o , 喜欢图书馆)
#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);
}
结果:1.00000000000
至少在概念上,pow(b, c)
可以被视为 exp(c * ln(b))
的实现。
因此,拦截潜在数值问题的一种方法是自己计算这个概念的第一部分pow
,使用
double i = ln(b);
如果 c * i
足够小(阈值将是您平台上使用的浮点方案的函数),那么您可以使用 C 标准继续评估 pow
库函数。
用 exp(c * i)
自己完成这项工作可能是不明智的,因为标准 pow
函数很可能有各种技巧来获得比 exp(c * ln(b))
更准确的结果.
我使用一个 C 库,它在两个 double
值上使用 pow
函数
double a = pow(b, c)
目前,我有b = 0.62
和c = 1504
,这意味着a
应该接近0(3.6e-312)。
但是我有一个浮点异常。如何避免直接return0?我们可以预见这种情况吗?
我使用 Debian 9,64 位,我用 gcc 6.3 编译。库是 ccmaes,这里是有问题的行:
https://github.com/CMA-ES/c-cmaes/blob/eda8268ee4c8c9fbe4d2489555ae08f8a8c949b5/src/cmaes.c#L893
我用过gdb所以浮点异常不是来自除法(t->chiN = 2.74)
如果我尝试重现它,使用 FPE 发生时的值,我没有问题(编译选项:-fopenmp -O3 -DNDEBUG -fPIC -Wall -Wextra -Wno-long-long -Wconversion -o , 喜欢图书馆)
#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);
}
结果:1.00000000000
至少在概念上,pow(b, c)
可以被视为 exp(c * ln(b))
的实现。
因此,拦截潜在数值问题的一种方法是自己计算这个概念的第一部分pow
,使用
double i = ln(b);
如果 c * i
足够小(阈值将是您平台上使用的浮点方案的函数),那么您可以使用 C 标准继续评估 pow
库函数。
用 exp(c * i)
自己完成这项工作可能是不明智的,因为标准 pow
函数很可能有各种技巧来获得比 exp(c * ln(b))
更准确的结果.