什么给出最好的精度，差异的指数或指数的商？

Question

给定一种典型的编程语言。

我有两个浮点数 a 和 b，它们彼此接近（即它们的差值绝对值比它们的平均值小得多）。

|a-b| << |a+b|/2

从数学上讲，我们有

exp(a-b) = exp(a)/exp(b) .

但是你编程的时候可以选择先计算(a-b)再求幂，或者先对a求幂，再对b求幂，然后再除。

如果 a 和 b 真的很接近，那么 (a-b) 的精度可能会很差。

示例

(1+pi*10^-20) -(1+1.1*pi*10^-20) = -pi*10^-21

但是如果你使用的是精度只有 19 位小数点的浮点数。您将得到零作为答案，这是一个很差的精度。您可以通过按如下方式重新排序操作来获得更好的精度

(1-1) + (pi*10^-20 -1.1*pi*10^-20) = -pi*10^-21

这将为您提供 -pi*10^-21，精确到小数点后 19 位。

因此我的问题是，给定一个有限的浮点精度，哪种计算 exp(a-b) 的方法精度更高？

差异的指数：

exp(a-b)

或指数的商

exp(a)/exp(b)

?

Answer 1

If a and b are really close to one another, then the precision on (a-b) might be bad.

相反，如果 a 和 b 接近，则 a - b 是精确的 (Sterbenz lemma)。因此 exp(a-b) 只涉及一个舍入步骤（假设一个正确舍入的 exp 函数以简化推理）。因此，exp(a-b) 为您尝试计算的表达式提供了正确舍入的结果。任何其他方式都不会比这更好。

exp(a-b) 明显优于替代方案的情况是 exp(a) 和 exp(b) 分别溢出，导致除法为 NaN。相比之下，exp(a-b) 仅在数学结果的最佳浮点近似值时产生 +inf，并且从不为有限 a 和 b.[=36= 产生 NaN ]

注意：您可能面临的问题是 a 和 b 的计算过于近似，以至于 a - b 的相对准确度非常糟糕。这是浮点数 (cancellation) 固有的问题，但在以任何方式计算 exp(a-b) 时，试图解决此问题已为时已晚：信息已经丢失。您只能为您拥有的 a 和 b 计算 exp(a-b)。正确的方法是使用 exp(a-b).

对于期望额外精度位为最终结果提供额外准确性的简单表达式，也可以通过计算具有额外精度的参考结果来凭经验测试这些东西。如果我以这种方式处理这个问题，我可能已经编写了 C 程序（对于我的平台，其中 FLT_EVAL_METHOD 被编译器定义为 0 并且 long double 是 IEEE 754 80 位双扩展） :

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

double best(double a, double b) {
  return exp(a-b);
}

double other(double a, double b) {
  return exp(a) / exp(b);
}

long double reference(long double a, long double b) {
  // assume the method doesn't matter so much with
  // long double computations, use any method:
  return expl(a-b);
}

int main(void) {
  for (int i = 0; i < 10; i++) {
    double a = rand() ^ ((long long)rand())<<16 ^ ((long long)rand()) << 32;
    a /= 0x1.0p64;
    double b = a * (1 + (double)rand() / (5.0 * RAND_MAX));
    printf("a=%a\nb=%a\n", a, b);

    double be = best(a, b);
    double o = other(a, b);

    long double r = reference(a, b);

    printf("error sub then exp:%La\n", fabsl(r - be));
    printf("error exp then div:%La\n", fabsl(r - o));
  }
}

上述程序的运行产生结果：

~ $ gcc ex.c && ./a.out
a=0x1.82def03cebc5p-2
b=0x1.a65bc869d479cp-2
error sub then exp:0xa.dp-60
error exp then div:0xa.dp-60
a=0x1.8164b7a7be6dep-6
b=0x1.b5b82c63fa0bcp-6
error sub then exp:0x8.1p-58
error exp then div:0x8.1p-58
a=0x1.88b779e295f18p-3
b=0x1.b1836e39a5186p-3
error sub then exp:0x8.5p-59
error exp then div:0xd.ecp-57
a=0x1.b5f437ec6fc4ap-6
b=0x1.e459d0a1d64b9p-6
error sub then exp:0xa.1p-59
error exp then div:0xd.7cp-57
a=0x1.889e1b20a7c9dp-3
b=0x1.8dddc5217d12fp-3
error sub then exp:0x9.4p-61
error exp then div:0xf.6cp-57
a=0x1.2d8f07147984cp-2
b=0x1.65acc944015e6p-2
error sub then exp:0x8.ep-58
error exp then div:0x8.ep-58
a=0x1.78b8432157465p-5
b=0x1.a9fd15e131562p-5
error sub then exp:0x9.4p-61
error exp then div:0xf.6cp-57
a=0x1.d214f566a0e1ep-2
b=0x1.0c90cb63e50c4p-1
error sub then exp:0xd.54p-58
error exp then div:0xd.54p-58
a=0x1.78dfa1c5c2ac4p-2
b=0x1.919d3723ae61p-2
error sub then exp:0x9.e8p-59
error exp then div:0x9.e8p-59
a=0x1.fb68c3c57fdd3p-2
b=0x1.103e0374df0bbp-1
error sub then exp:0xd.54p-58
error exp then div:0x9.56p-57

什么给出最好的精度，差异的指数或指数的商？

What gives best precision, exponential of the difference or the quotient of the exponentials?

floating-point

precision

exponentiation

mathematical-expressions