使用 glibc/libm 和 float32 的 atan2 的错误结果

Question

我目前正在为医疗设备开发固件，其中涉及很多困难的数学运算。目标处理器支持硬件中的浮点运算，但仅支持 float32（又名 single）。

为了模拟行为并证明我的公式和代码的正确性，我已将固件的相关/数学部分移植到 Linux 中的 GCC 工具链（gcc 6.3.0，libc6 2.24），仔细检查 float32 在任何地方都被使用，并且没有使用编译器开关，这可能会降低数学运算的精度或标准兼容性；值得注意的是，有 none 个 -ffast-math 或其朋友。

现在，对于一小组输入参数，我得到了意想不到的结果。我已经追踪到这个问题并得出结论，libm 为非常小的一组输入参数计算了 arctan（准确地说：atan2）的错误结果。

例如，如果我有

#include <math.h>

#define C_RAD2DEG (57.29577951308f)

int main(void)
{
  float f_Temp = C_RAD2DEG * atan2f(0.713114202f, 0.665558934f);
}

f_Temp 被计算为 46.9755516f，其中正确的结果将是 46.975548972f.

请注意，我通常了解不同浮点数据类型、舍入错误等问题。

然而，我的感觉是，即使考虑到 float32 的低精度，上面显示的误差也高了一个数量级，不幸的是，对于随后的计算，该误差太大了。

此外，atan2 函数的可能输入参数中只有很小一部分受到此问题的影响。

任何人都可以快速解释一下这是 libm 中的错误，还是仅仅是由于 float32 的不精确以及计算 [=17= 所需的大量顺序操作]?

Answer 1

您作为观察结果报告的数字 46.9755516f 对应于 float 值 46.975551605224609375。

您作为预期结果报告的数字 46.975548972f 对应于 float 值 46.97554779052734375。

这些是相邻的 float 值，这意味着它们相差 1 个最小精度单位 (ULP)。（它们的差异是 3.814697265625e-06，这是当最高有效位的值为 32 时 float 有效位中最低有效位的值，就像 47 左右的数字一样。）这是最小的可能数量float 可以按该比例变化。

通常，数学库例程很难实现，并且没有人用正确的舍入（舍入到最接近精确数学值的可表示数字）和已知的有界运行时间来实现所有这些例程。一些 ULP 错误在三角函数例程中并不罕见。

即使您使用的 libc 代码提供了正确的舍入结果，将其从弧度转换为度数也会引入两个舍入错误（将 180/π 转换为可表示的值并乘以它）。期望最终结果是最接近理想数学结果的 float 是不合理的；你应该期待几个 ULP 错误。

使用 glibc/libm 和 float32 的 atan2 的错误结果

Wrong result for atan2 with glibc / libm and float32

c

precision

libc

libm

gcc6