我认为 Higham 的数值算法的准确性和稳定性解决了如何分析这些类型的问题。参见第 2 章，尤其是练习 2.8。

在这个答案中，我想指出海厄姆的书中没有真正解决的问题（就此而言，它似乎并不广为人知）。如果您对证明诸如此类的简单数值算法的属性感兴趣，您可以使用现代 SMT 求解器（可满足性模理论）的强大功能，例如z3，使用Haskell中的sbv等包。这比使用铅笔和纸要容易一些。

假设给定 $0\leq x\leq y$，我想知道 $z=(x+y)/2$ 是否满足 $x\leq z\leq y$。下面的 Haskell 代码

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

会让我自动执行此操作。这test1 fun是命题$x \leq \mathit{fun}(x,y) \leq y$ 对于所有有限浮点 $x,y$ 和 $0\leq x\leq y$。

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

它溢出了。假设我现在采用您的另一个公式：$z=x/2+y/2$

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

不起作用（由于逐渐下溢：$(x/2)\times2\neq x$，由于所有算术都是以 2 为底的，这可能不直观）。

现在试试 $z=x + (yx)/2$：

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

作品！这Q.E.D.是一个证明该test1属性适用于上面定义的所有浮点数。

相同的情况如何，但仅限于 $x\leq y$（而不是 $0\leq x\leq y$）？

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

好的，如果 $yx$ 溢出，那么 $z = x + (y/2-x/2)$ 怎么样？

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

因此，似乎在我在这里尝试过的公式中， $x + (y/2 - x/2)$ 似乎有效（也有证明）。在我看来，SMT 求解器方法比用铅笔和纸进行浮点错误分析更快地回答了对简单浮点公式的怀疑。

最后，准确性和稳定性的目标往往与性能目标不一致。就性能而言，我真的看不出你能比 $(x+y)/2$ 做得更好，特别是因为编译器仍然会为你完成将其翻译成机器指令的繁重工作。

PS这都是单精度IEEE754浮点运算。我用双精度算术（替换SFloat为SDouble）检查了 $x \leq x + (y/2-x/2) \leq y$，它也可以工作。

PPS在代码中实现这一点时要记住的一件事是编译器标志-ffast-math（某些形式的此类标志有时在某些常见编译器中默认打开）不会导致 IEEE754 算术，这将使上述证明无效。如果您确实使用了启用例如关联加法优化的标志，那么除了$(x+y)/2$ 之外，没有任何意义。

PPPS只看没有条件的简单代数表达式，我有点得意忘形。严格来说，Don Hatch的公式更好。

首先，请注意，如果您有一种方法可以在所有情况下给出最准确的答案，那么它将满足您的要求。（请注意，我说的是最准确的答案，而不是最准确的答案，因为可能有两个获胜者。） 证明：相反，如果你有一个不满足要求条件的尽可能准确的答案，那意味着要么answer<min(x,y)<=max(x,y)（在这种情况下min(x,y)是一个更好的答案，一个矛盾），或者min(x,y)<=max(x,y)<answer（在这种情况下max(x,y)是一个更好的答案，一个矛盾）。

所以我认为这意味着你的问题归结为找到一个最准确的答案。假设始终采用 IEEE754 算法，我提出以下建议：

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

我认为这给出了最准确答案的论点是一个有点乏味的案例分析。开始：

• 案例max(abs(x),abs(y)) >= 1.：

  • 子情况 x 和 y 都不是非规范化的：在这种情况下，计算的答案操作相同的尾数，因此如果我们假设扩展指数以防止溢出，则x/2.+y/2.给出与计算结果完全相同的答案。(x+y)/2这个答案可能取决于舍入模式，但在任何情况下，IEEE754 都保证它是一个最佳答案（因为计算x+y得到的保证是数学 x+y 的最佳近似值，并且除以 2 是准确的案件）。

  • 子案例 x 是非规范化的（因此abs(y)>=1）：

    answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

  • 子案例 y 是非规范化的（因此abs(x)>=1）：类似。

• 案例max(abs(x),abs(y)) < 1.：
  • 计算x+y的子情况是非非规范化或非规范化和“偶数”：尽管计算x+y可能不准确，但 IEEE754 保证它是数学 x+y 的最佳近似值。在这种情况下，表达式中随后除以 2(x+y)/2.是精确的，因此计算出的答案(x+y)/2.是数学 (x+y)/2 的最佳近似值。
  • 子情况下计算x+y的非规范化和“奇数”：在这种情况下，x,y 中的一个也必须非规范化和“奇数”，这意味着 x,y 的另一个用相反的符号非规范化，因此计算x+y的是恰好是数学 x+y，因此(x+y)/2.IEEE754 保证计算结果是数学 (x+y)/2 的最佳近似值。

对于 IEEE-754 二进制浮点格式，以binary64（双精度）计算为例，S. Boldo 正式证明了下面显示的简单算法提供了正确的舍入平均值。

Sylvie Boldo，“计算浮点平均值的程序的形式验证”。在形式工程方法国际会议上，第 17-32 页。Springer, Cham, 2015。（在线草稿）

由于在二进制浮点算术中除以二是准确的，只要不发生下溢，通过选择$(x+y)/2$和$x/2 + y/2$这两个公式之一似乎很直观在适当的情况下（基于输入的大小）一次应该达到一个准确四舍五入的平均值。Boldo 的论文表明，对于 IEEE-754，binary64任何切换点$C \in [2^{-967}, 2^{970}]$就足够了。人们可能会选择$C$以便为特定用例提供最佳性能。

这会产生以下示例ISO-C99代码：

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

在最近的后续工作中，S. Boldo 和合著者展示了如何通过使用融合乘加 (FMA) 操作和众所周知的精度来实现 IEEE-754 十进制浮点格式的最佳结果。加倍构建块（TwoSum）：

Sylvie Boldo、Florian Faissole 和 Vincent Tourneur，“一种正式证明的算法，用于计算十进制浮点数的正确平均值”。在第 25 届 IEEE 计算机算术研讨会 (ARITH 25) 上，2018 年 6 月，第 69-75 页。（网上草稿）

虽然它在性能方面可能不是超级高效，但有一种非常简单的方法可以（1）确保没有一个数字大于x或y（没有溢出）和（2）保持浮点“准确”为可能（和(3)，作为额外的奖励，即使使用减法，也不会将任何值存储为负数。

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

事实上，如果你真的想要精确，你甚至不需要当场进行除法；只需返回 和 的值min(x, y)，difference您可以使用这些值在逻辑上简化或稍后进行操作。