Kahan算法：如何获得产品的确切差异

我最近返回浮点错误分析，以完善下一版《 基于物理的渲染》一书中的一些细节。 浮点数是一个有趣的计算领域，充满了惊喜（好坏），以及摆脱不愉快惊喜的棘手技巧。

在此过程中，我遇到了StackOverflow上的这篇文章 ，从中我了解了一种用于精确计算的优雅算法。 $$$a \倍b-c \倍d$。

但是在继续算法之前，您需要了解表达式中的狡猾之处 $$$a \倍b-c \倍d$？ 拿 $$$a = 33962.035$， $$$b = -30438.8$， $$$c = 41563.4$和 $$$d = -24871.969$。 （这些是我在pbrt启动期间获得的实际值。）对于32位浮点值，我们得到： $$$a \乘以b = -1.03376365 \乘以10 ^ 9$和 $$$c \倍d = -1.03376352 \倍10 ^ 9$。 我们进行减法，我们得到 $$$-128$。 但是，如果您以双精度执行计算，最后将它们转换为浮点数，则会得到 $$$-75.1656$。 发生什么事了

问题在于每件作品的价值都远远超出了底线 $$$-1 \乘以10 ^ 9$，其中可表示的浮点值之间的距离非常大-64。即，在四舍五入时 $$$a \乘以b$和 $$$c \倍d$单独地，对于最接近的可表示浮点数，它们将变成64的倍数的数字。反过来，它们的差将是64的倍数，并且不希望它将变为 $$$-75.1656$比 $$$-64$。 在我们的案例中，由于对两个作品进行了四舍五入，结果甚至更进一步 $$$-1 \乘以10 ^ 9$。 我们将直接遇到旧的灾难性灾难减少¹ 。

这是比²更好的解决方案：

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts()计算 $$$a \倍b-c \倍d$以一种避免灾难性收缩的方式。 传奇人物William Kahan在文章“没有超精确算法的浮点计算的成本”中首次描述了这种技术。 值得注意的是，Kahan的著作整体上看起来很有趣，他们对浮点世界的当前状态以及数学和技术上的考虑有很多评论。 这是他的发现之一：

我们中那些为浮点数算法的变迁而苦苦挣扎，而对编译器的“优化”却没有深思熟虑的人，可以为这场斗争中的胜利感到自豪。 但是，如果我们将这场战斗的延续继续传递给子孙后代，这将与整个文明的本质相矛盾。 我们的经验表明，编程语言和开发系统是我们必须处理的太多混乱的根源。 可以省去太多的错误，以及一些吸引人的“优化”方法，这些方法对于整数是安全的，但有时对浮点数却是致命的。

向他的智慧表示敬意之后，让我们回到DifferenceOfProducts() ：精通此功能的基础是使用FMA ³融合乘加法。 从数学的角度来看， FMA(a,b,c)为 $$$a \乘以b + c$因此，起初看来，该操作仅对微优化有用：一条指令而不是两条指令。 但是fma
它具有特殊的属性-仅舍入一次。

通常 $$$a \乘以b + c$首先计算 $$$a \乘以b$，然后将此值（通常不能以浮点格式表示）四舍五入到最接近的浮点数。 然后向此舍入值添加 $$$c$，并且此结果再次四舍五入到最接近的浮点数。 FMA的实现方式是仅在末尾进行舍入-中间值 $$$a \乘以b$保持足够的精度，因此添加后 $$$c$最终结果将最接近真实值 $$$a \乘以b + c$浮动值。

处理完FMA之后，我们将返回DifferenceOfProducts() 。 再次，我将显示它的前两行：

  float cd = c * d; float err = std::fma(-c, d, cd); 

第一个计算舍入值 $$$c \倍d$第二个...减去 $$$c \倍d$从他们的工作？ 如果您不知道FMA的工作原理，那么您可能会认为err永远为零。 但是在使用FMA时，第二行实际上是在计算值中提取舍入误差的值 $$$c \倍d$并将其保存到err 。 之后，输出非常简单：

  float dop = std::fma(a, b, -cd); return dop + err; 

第二个FMA使用FMA计算作品之间的差异，仅在最后进行舍入。 因此，它可以防止灾难性减少，但需要以四舍五入的值工作。 $$$c \倍d$。 return “修补”此问题，并在第二行中添加了突出显示的错误。 在Jeannenrod等人的文章中。 结果表明 ，结果高达1.5 ulps（单位精度度量值）时是真实的，这非常好：FMA和简单的浮点运算在高达0.5 ulps的范围内都是准确的，因此该算法几乎是完美的。

我们用新的锤子

当您开始寻找使用DifferenceOfProducts() ，它确实非常有用。 计算二次方程的判别式？ 调用DifferenceOfProducts(b, b, 4 * a, c) ⁴ 。 计算2x2矩阵的行列式？ 该算法将解决此问题。 在下一版本的pbrt中，我发现它有大约80种用途。 在所有这些中，向量乘积的功能是最受欢迎的。 它一直是问题的根源，因此，您必须举手并在实现中使用double以避免灾难性的减少：

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

现在，我们可以继续使用float并使用DifferenceOfProducts() 。

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

从帖子开头的那个狡猾的例子实际上是矢量工作的一部分。 在某个阶段， pbrt代码需要计算向量的向量积 $$$（33962.035、41563.4、7706.415）$和 $$$（-24871.969，-30438.8，-5643.727）$。 使用float进行计算时，我们将得到一个向量 $$$（1552，-1248，-128）$。 （一般规则：如果在涉及大数的浮点计算中，得到的整数值不是那么大，那么这几乎可以肯定地表明发生了灾难性的减少。）

矢量乘积是双精度的 $$$（1556.0276，-1257.5151，-75.1656）$。 我们看到了浮动 $$$x$看起来很正常 $$$y$已经很糟糕了 $$$z$成为灾难，成为寻求解决方案的动力。 而使用DifferenceOfProducts()和float值会得到什么结果呢？ $$$（1556.0276，-1257.5153，-75.1656）$。 价值观 $$$x$和 $$$z$对应于双精度，并且 $$$y$稍微偏移-因此多余的ulp来自。

速度呢？ DifferenceOfProducts()执行两个FMA以及乘法和加法。 天真的算法可以用一个FMA和一个乘法来实现，这似乎需要花费一半的时间。 但是实际上，事实证明，从寄存器中获取值之后，这并不重要：在我笔记本电脑上的综合基准测试中， DifferenceOfProducts()仅比朴素算法贵1.09倍。 双精度运算速度慢了2.98倍。

一旦您了解了减少灾难性灾难的可能性，代码中所有看起来无辜的表达式就开始显得可疑。 DifferenceOfProducts()对于其中大多数似乎是一个很好的解决方法。 它易于使用，我们没有特殊的理由不使用它。

注意事项

1. 当减去具有不同符号的数量或添加具有相同符号的值时，灾难性的减少不是问题。 但是，在添加具有不同符号的值时可能会成为问题。 也就是说，应该以与差异相同的怀疑来考虑金额。
2. 作为读者的练习，我建议编写SumOfProducts()函数，该函数可防止灾难性收缩。 如果要使任务复杂化，请解释为什么在DifferenceOfProducts() ， dop + err == dop ，因为符号a*b和c*d不同。
3. FMA指令在GPU上已经可用了十多年，在大多数CPU上已经使用了至少五年。 使用std::fma()时，可能有必要将编译器标志添加到CPU以直接生成它们。 在gcc和clang中， -march=native作品。
4. 在IEEE浮点格式中，精确执行乘以2的幂，因此4 * a不会引起任何舍入误差，除非发生溢出。