如何计算两个浮点数的平均值？

Question

在浮点运算中计算两个数字的平均值最准确的方法是什么？让我们考虑一下最常见的双精度 64 位数字。

1. (a + b) / 2

2. a / 2 + b / 2

3. a + (b - a) / 2

这些计算平均值的方法可能会给出不同的结果，如下面的 C++ 代码所示：

double a = 1.2;
double b = 3.6;
double mean1 = (a + b) / 2.0;
double mean2 = a / 2.0 + b / 2.0;
double mean3 = a + (b - a) / 2.0;
cout << fixed << setprecision(20);
cout << "mean1: " << mean1 << endl;
cout << "mean2: " << mean2 << endl;
cout << "mean3: " << mean3 << endl;

产生结果

mean1: 2.39999999999999991118
mean2: 2.39999999999999991118
mean3: 2.40000000000000035527

问题 计算两个浮点数的平均值的数学性质？ 比较情况 1 和 2。简单的结论是精度相同，但情况 1 可以溢出，而情况 2 可以下溢。

mean1

和

mean2

似乎提供了更准确的结果，因为它们比

mean3

更接近真实值 2.4。似乎总是

mean3

不太精确。直观上这是有道理的，因为计算mean3需要3次运算，因此有更多机会丢失精度，而计算mean1只需要2次运算。

我偶然发现了 Fortran 代码中的第三种情况 https://www.cs.umd.edu/~oleary/LBFGS/FORTRAN/linesearch.f

stpf = stpc + (stpq - stpc)/2

stp = stx + p5*(sty - stx)

其中

p5

是 0.5

我不确定这种计算平均值的动机是什么。可能是 n=2 运行平均值的情况，如浮点值的数值稳定运行平均值或https://diego.assencio.com/?index=c34d06f4f4de2375658ed41f70177d59中所述，但它提高了准确性仅适用于大量元素。

如果这是为了避免溢出，那么它似乎不起作用。所有三种情况都可能产生上溢或下溢。

#include <iostream>
#include <iomanip>

using namespace std;

double mean1(double a, double b) {
    return (a + b) / 2.0;
}

double mean2(double a, double b) {
    return a / 2.0 + b / 2.0;
}

double mean3(double a, double b) {
    return a + (b - a) / 2.0;
}

void test(double a, double b) {
    cout << setprecision(20);
    cout << "a: " << a << ", b: " << b << endl;
    cout << "mean1: " << mean1(a, b) << endl;
    cout << "mean2: " << mean2(a, b) << endl;
    cout << "mean3: " << mean3(a, b) << endl;
    cout << endl;
}

int main() {
    test(1e+308, 1e+308); // case 1 overflow
    test(4e-324, 4e-324); // case 2 underflow
    test(-1e+308, 1e+308); // case 3 overflow
}

a: 1.000000000000000011e+308, b: 1.000000000000000011e+308
mean1: inf
mean2: 1.000000000000000011e+308
mean3: 1.000000000000000011e+308

a: 4.9406564584124654418e-324, b: 4.9406564584124654418e-324
mean1: 4.9406564584124654418e-324
mean2: 0
mean3: 4.9406564584124654418e-324

a: -1.000000000000000011e+308, b: 1.000000000000000011e+308
mean1: 0
mean2: 0
mean3: inf

此外，第三种情况的性能稍差一些（但是应谨慎考虑结果，因为它们可能会根据上下文而有很大差异，但通常指令数量较少的代码运行速度更快）

static void case_1(benchmark::State& state) {
  volatile double a = 1.2;
  volatile double b = 3.6;
  for (auto _ : state) {
    benchmark::DoNotOptimize((a + b) / 2.0);
  }
}
BENCHMARK(case_1);

static void case_2(benchmark::State& state) {
  volatile double a = 1.2;
  volatile double b = 3.6;
  for (auto _ : state) {
    benchmark::DoNotOptimize(a / 2.0 + b / 2.0);
  }
}
BENCHMARK(case_2);

static void case_3(benchmark::State& state) {
  volatile double a = 1.2;
  volatile double b = 3.6;
  for (auto _ : state) {
    benchmark::DoNotOptimize(a + (b - a) / 2.0);
  }
}
BENCHMARK(case_3);

https://quick-bench.com/q/V_STrIZ5CmIL0Q3KpbLX6G9-G90

我可能缺少这个选择背后的一些逻辑。上面的代码可能是要在不支持 IEEE 758 的不同硬件上运行？那么，是否有任何理由将平均值计算为

a + (b - a) / 2

而不是

(a + b) / 2

？或者只是有点不准确，应该使用简单的方法来计算平均值？

Answer 1

为了准确性，应首选

(a + b) / 2

或

a / 2 + b / 2

，因为除以 2 在浮点中是无损，除非它“下溢”。

/ 2

只是递减指数，因此，实际上，错误仅在加法中引入。

如果您愿意，可以用

* 0.5

代替

/ 2.0

— 这没有区别。

表格

a + (b - a) / 2

，但是保证：

• 如果

  a

  和

  b

  都是非负有限值，那么结果也是非负有限值。

  (a + b)/2

  表格不提供此保证；和
• 如果

  a

  和

  b

  都是非负有限值，则结果始终是

  >= a

  和

  <= b

  。其他两种形式都没有提供这种保证。