exp10 与 pow(10) 不同

Question

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

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

我的函数中令我惊讶的一个组件是 exp10 的输出；据我所知，只要两个双精度之间的间距小于 2，那么存储为双精度的整数值应该是精确的 - 尽管 10^14 正在推动它，这应该是一个精确的整数 em>（因为 10^14 =~ 2^46.507

我的调试工作的摘录（没有什么是显而易见的），输出如下：

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: %.16lf\n", exp);
printf("powTen: %.16lf\n", powTen);
printf("powTen2: %.16lf\n", powTen2);
printf("powTen3: %.16lf\n", powTen3);
printf("powTen4: %.16lf\n", powTen4);

//these two are exact
printf("10^14: %.16lf\n", exp10(14));
printf("powTen5: %.16lf\n", powTen5);
printf("exp == 14.0: %d\n", 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 仍然行为不端。

Answer 1

很可能您的exp10 实现行为不端。请注意，它返回的结果有时会偏离 ulp（相对于 10^14 为 0.125。）

这是一个相当令人发指的错误；您有一个案例，正确答案可以表示为 double 但 exp10 没有这样做。

我会回应 Ben Voigt 的评论，即编译器有时可能会自行评估事物，而不是将它们传递给数学库。它可能做得更好，因为它可能链接到任意精度的数学库。您可以尝试使用 -fno-builtin 选项，看看它是否有任何改变。

不幸的是，我认为crlibm 没有实现exp10。否则，我建议您使用它并停止担心。

编辑：我拥有的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);
}

不要指望这会很好。