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

问：

在浮点运算中计算两个数字的平均值的最准确方法是什么？让我们考虑最常见的双精度 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并且似乎提供了更准确的结果，因为它们比真实值 2.4 更接近。似乎总是情况不那么精确。直观地说，这是有道理的，因为 mean3 的计算需要 3 次运算，因此失去精度的机会更多，而 mean1 的计算只需要 2 次运算。mean2mean3mean3

更新：1.2 和 3.6 存储为

1.1999999999999999555910790149937383830547332763671875
3.600000000000000088817841970012523233890533447265625

所以他们真正的意思是

2.40000000000000002220446049250313080847263336181640625

但是，此值不能用双倍数表示。最接近的双精度值恰好是 ，其中 和 aremean1mean1mean3

2.399999999999999911182158029987476766109466552734375 (difference 1.11e-16)
2.4000000000000003552713678800500929355621337890625   (difference 3.33e-16)

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

stpf = stpc + (stpq - stpc)/2

stp = stx + p5*(sty - stx)

其中为 0.5p5

我不确定这种计算均值方式的动机是什么。可能是运行平均值的 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 754 的不同硬件上运行？那么，有什么理由计算均值而不是？或者它只是轻微的不准确，应该使用简单的方法来计算平均值？a + (b - a) / 2(a + b) / 2