exp10 different to pow(10)

exp10 different to pow(10)

本文关键字：pow different to exp10      更新时间：2023-10-16

首先，我意识到大多数以10为基数的数字不能精确地以2为基数表示，所以我的问题并不是关于浮点算术的缺陷。

我正试图写一个函数，该函数将尝试通过检查最后6个有意义的数字是否在一定范围内，并将其更改为下一个可表示的一些假定的精确值(仅用于显示目的-除非它是整数或2的幂)来纠正被累积舍入误差污染的双精度。

我的函数的一个组件让我惊讶的是exp10的输出;据我所知，只要两个双精度之间的间距小于2，那么存储为双精度的整数值应该是精确的-尽管10^14正在推动它，这应该是一个精确的整数(因为10^14 =~ 2^46.507 <2 ^ 53)。然而，这不是我的测试显示。

我调试工作的摘录(没有什么明显的)和输出如下:

double test = 0.000699;
double tmp = fabs(test);
double exp = 10.0 - floor(log10(tmp));
double powTen = exp10(10.0 - floor(log10(tmp)));
double powTen2 = exp10(exp);
double powTen3 = exp10((int)exp);
double powTen4 = exp10(exp);
double powTen5 = pow(10, exp);
printf("exp: %.16lfn", exp);
printf("powTen: %.16lfn", powTen);
printf("powTen2: %.16lfn", powTen2);
printf("powTen3: %.16lfn", powTen3);
printf("powTen4: %.16lfn", powTen4);
//these two are exact
printf("10^14: %.16lfn", exp10(14));
printf("powTen5: %.16lfn", powTen5);
printf("exp == 14.0: %dn", exp == 14.0);
输出:

exp: 14.0000000000000000
powTen: 100000000000000.1250000000000000
powTen2: 100000000000000.1250000000000000
powTen3: 100000000000000.1250000000000000
powTen4: 100000000000000.1250000000000000
10^14: 100000000000000.0000000000000000
powTen5: 100000000000000.0000000000000000
exp == 14.0: 1

pow得到确切的答案，exp10与硬编码的int一样。对于所有其他情况，我加1/8(10^14和10^14 +下一个可表示的间隔是1/64)。文档中说exp10应该等于pow。有人能看出我漏掉了什么吗?

编辑-与O3, O2, O1优化我得到预期的输出- 除非数据不能知道，直到运行时。此时，exp10仍然不正常

很可能您的exp10实现行为不正常。请注意，它返回的结果有时会偏离一个ulp(相对于您的10^14，这是0.125)

这是一个相当可恶的错误;你有一个情况，正确的答案是可以表示为double，但exp10不这样做。

我赞同Ben Voigt的评论，编译器有时可能自己求值，而不是把它们传递给数学库。它可能做得更好，因为它可能链接到任意精度的数学库。您可以尝试-fno-builtin选项，看看它是否改变了什么。

不幸的是，我认为clibm没有实现exp10。否则，我建议你使用它，不要担心。

EDIT:我所拥有的eglibc源的副本似乎实现了exp10:

double
__ieee754_exp10 (double arg)
{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_exp (M_LN10 * arg);
}

相关文章：