c - 如何比较两个数学库实现？

Question

如您所知，C 标准库定义了几个标准函数调用，应该由任何兼容的实现来实现，例如 Newlib、MUSL、GLIBC ...

例如，如果我以 Linux 为目标，我必须在 glibc 和 MUSL 之间进行选择，而我的标准是数学库的准确性libm。例如，我如何比较两种可能的sin()实现cos()？

一种天真的方法是使用参考输入（例如来自 Matlab）在一组随机生成的输入上测试两种实现结果的输出质量，但是还有其他更可靠/正式/结构化/引导的方法来比较/两个模型？我试图看看在这个方向上是否有任何研究，但我发现了，所以任何指针都值得赞赏。

Answer 1

一些想法：

• 您可以使用GNU 多精度算术库 (GnuMP来生成良好的参考结果。
• 可以详尽地测试大多数（如果不是全部）单参数单精度（IEEE-754 binary32）例程。（对于一些 macOS 三角函数例程，例如sinf，我们详尽地测试了一种实现，验证它是否返回了如实舍入的结果，这意味着结果是数学值 [如果可表示] 或两个相邻值之一 [如果不是]。然后，在更改实现时，我们比较了一个，如果新实现的结果与旧实现的结果相同，则通过。否则，使用GnuMP进行测试。由于新实现与旧实现基本一致，因此导致在 GnuMP 的几次调用中，如果我没记错的话，我们能够在大约三分钟内彻底测试一个新的例程实现。）
• 详尽地测试多参数或双精度例程是不可行的。
• 在比较实现时，您必须选择一个或多个指标。具有良好的最坏情况错误的库有利于证明；可以断言它的界限对任何参数都成立，并且可以用来在后续计算中推导出进一步的界限。但是，具有良好平均误差的库可能会为使用大量数据的物理模拟产生更好的结果。对于某些应用程序，只有“正常”域中的误差可能是相关的（大约 -2π 到 +2π 的角度），因此减少大参数（高达 10 ³⁰⁸左右）的误差可能是不相关的，因为这些参数从未使用过。
• 有一些共同点应该测试各种例程。例如，对于三角函数例程，测试不同的 π 分数。除了在数学上有趣之外，这些往往是实现在内部在近似值之间切换的地方。还要测试大量可表示但恰好非常接近 π 的简单分数的倍数。这些是减少参数的最坏情况，如果没有正确完成，可能会产生巨大的相对错误。他们需要数论才能找到。以任何一种散漫的方法进行测试，甚至是没有考虑这种归约问题的有序方法，都无法找到这些麻烦的论据，因此很容易将具有巨大错误的例程报告为准确的。
• 另一方面，如果没有实现的内部知识，则无法知道要测试的重要点。例如，在设计一个正弦例程时，我会使用 Remez 算法来找到一个极小极大多项式，目标是使其从 –π/2 到 +π/2（对于这类事情来说相当大，但是仅举例）。然后我会看看在参数减少期间可能发生的算术和舍入错误。有时他们会在这个区间之外产生一个结果。所以我会回到极小多项式生成并推动稍大的间隔。而且我还会寻求减少论点的改进。最后，我会得到一个保证在某个区间内产生结果的归约，以及一个已知在该区间内具有一定精度的多项式。为了测试我的例程，您需要知道该区间的端点，并且您必须能够找到一些参数，为这些参数减少在这些端点附近产生点，这意味着您必须了解我的参数减少的方式已实现——它使用了多少位，等等。就像上面提到的麻烦的论点一样，这些点不能用散漫的方法找到。但与上述不同的是，它们无法从纯数学中找到；您需要有关实施的信息。这使得实际上不可能知道您已经比较了实现的最糟糕的潜在参数。这意味着您必须了解我的参数缩减是如何实现的——它使用了多少位，等等。就像上面提到的麻烦的论点一样，这些点不能用散漫的方法找到。但与上述不同的是，它们无法从纯数学中找到；您需要有关实施的信息。这使得实际上不可能知道您已经比较了实现的最糟糕的潜在参数。这意味着您必须了解我的参数缩减是如何实现的——它使用了多少位，等等。就像上面提到的麻烦的论点一样，这些点不能用散漫的方法找到。但与上述不同的是，它们无法从纯数学中找到；您需要有关实施的信息。这使得实际上不可能知道您已经比较了实现的最糟糕的潜在参数。