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

Question

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

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

现在，事实证明，对于一小组输入参数，我得到了意想不到的结果。我已经追踪到问题并得出结论，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。

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

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

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

谁能简短地解释一下这是 libm 中的错误，还是仅仅是由于 float32 的不精确以及计算所需的大量顺序操作atan2?

Answer 1

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

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

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

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

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