浮点乘法的错误界限答案

【问题标题】：Error bound on floating point multiplication浮点乘法的错误界限
【发布时间】：2019-12-14 04:15:41
【问题描述】：

我正在寻找一个公式来限制浮点乘法传播的错误。我已经看到以下陈述（但未证明）。

如果xa 是x 的浮点逼近，绑定ux ulps，同样ya 是y 的浮点逼近，绑定uy ulps 和p 是浮点点精度，则正确舍入乘积 xa*ya 的界限为 ux + uy + ux*uy/(2^p) + 0.5 ulps。

我已经看到了以下更简单的结果。如果ya 是y 与边界uy ulps 的近似值，则正确舍入乘积x*ya 上的边界是uy + 0.5 ulps。

所以我想我正在寻找第一个结果的证明，或者至少是确认它是正确的。

TIA

编辑

接受了 Mark Dickinson 在下面的 cmets 中所说的内容并重新阅读了 Florian Loitsch 的论文（来自上面的第二个公式），很明显，这个问题中的 ulps 与潜在的非规范化结果有关。例如，使用 Mark 给出的示例并通过非 IEEE 方法计算乘积，即乘以尾数并丢弃低 53 位（但必要时舍入），我们得到

 x   = 0x1.8000000000000p-1
 y   = 0x1.4000000000000p+0
ya   = 0x1.3fffffffffff6p+0
x*y  = 0x0.f000000000000p+0 (denormalised)
x*ya - 0x0.ffffffffffff8p+8 (denormalised)

这是 8 ulp 的误差，并且在第二个公式的范围内。

【问题讨论】：

两种结果都不正确：当您从一个 binade 的底部移动到下面一个 binade 的顶部时（同时保持相对误差常数）。对于更简单结果的反例，假设 IEEE 754 binary64 格式和语义，x = 0.75、y = 1.25 和uy = 10（比如ya = 1.25 - 10 * 2**-52）。然后x * ya 与真实结果x * y 相差15 ulps。（舍入模式不相关，因为两个乘法都是精确的。）所以我很想看看第二个语句的“证明”。
证明在 Florian Loitsch 论文 'Printing Floating-Point Numbers Quick and Accurately with Integers'，引理 3.4。这是一篇众所周知的论文，所以我可能歪曲了它所说的内容。
@MarkDickinson 你说同时保持相对误差不变是什么意思？我相信我引用的公式假设误差范围与结果相关。因此，在您引用的示例中，是针对 ya 或 x*ya 测量 15 ulps 的误差？
我的意思是 15 * ulp(x*ya) 的结果（查看 ulps 的 y 是没有意义的，因为 x*ya 的大小可能完全不同）。我承认我对那篇论文中发生的事情感到困惑。有类似“Since, for x⊗y, 1 ulp = 2^(q+ex+ey) ...”之类的陈述（在引理 3.3 的证明中）似乎不正确，除非 ulp 有一些非常规的定义继续 - 鉴于 x 的 binade 和 y 的 binade 知识，x*y 可能会出现两个可能的 binade，因此 x*y 的 ulp 有两种可能性。有什么我想念的想法吗？
啊，我想我明白了，仔细阅读：对于“diy_fp”类型，乘法结果不一定必须归一化，所以我认为有一些技巧在 ulp 的定义中进行（这不是很明确）。因此，确实在那篇论文中使用“ulp”并不完全符合标准用法（或者更确切地说，任何标准用途，因为对于 ulp 的确切含义并没有普遍共识，尤其是在基数边界附近）。这是some reading 的主题。

标签： floating-point floating-accuracy ieee-754

【解决方案1】：

如果你假设 ulps 是相对于近似值的（我相信另一种情况更困难，而且用处不大），GNU MPFR 的doc/algorithms.tex 文档有以下结果的证明，这里重写遵循您的符号并简化以获得独特的精度 p。

如果你有

|xa - x| ≤ ux · ulp(xa)
|ya - y| ≤uy·ulp(ya)

那么四舍五入到最接近的乘积的误差为 [(2 + 2^1−p)(ux + uy) + 1/2] · ulp(xa*ya)。

这通常比您给出的界限大，大约是 2 倍。但是，您给出的界限的有效性需要确认。通常很难找到严格的误差范围，因为当 ulp 以 2 的幂次变化 2 倍时，很容易失去因子 2（摆动），这可能是由于次优证明或实际精度损失极端案例。

【讨论】：