ARM64 上浮点精度的奇怪问题
Strange issue with floating point accuracy on ARM64
我 运行 遇到了一个非常奇怪的问题,即 ARM64 上的浮点精度。我有一段非常简单的 C++ 代码,看起来像这样:
float sx = some_float_number_1;
float sy = some_float_number_2;
float ex = some_float_number_3;
float ey = some_float_number_4;
float px = ex;
float py = ey;
float d1 = (ex - sx) * (py - sy);
float d2 = (px - sx) * (ey - sy);
float d = d1 - d2;
float t = (ex - sx) * (py - sy) - (px - sx) * (ey - sy);
//32-bit output: d == t == 0
//64-bit output: d == 0, t != 0
理论上,d 应该等于 t 且等于 0,这正是 32 位 ARM 上发生的情况。但由于某些奇怪的原因,t 的输出在 64 位 ARM 上不等于 0,而 d 仍然正确。我从来没有见过这样的错误,所以我不知道是什么导致了这种问题。
编辑:更多信息
- 如果你没有注意到,d 和 t 的输出应该都是 0,因为 (ex - sx) * (py - sy) 等于 (px - sx) * (ey - sy )
- 此问题仅在输入的小数部分不等于 0 时发生。
- 我使用的编译器是 Android NDK r15c 包中包含的 Clang。
EDIT2:这是反汇编
4c: 52933348 mov w8, #0x999a // #39322
50: 72a82828 movk w8, #0x4141, lsl #16
54: b90683e8 str w8, [sp,#1664]
58: 52933348 mov w8, #0x999a // #39322
5c: 72a82728 movk w8, #0x4139, lsl #16
60: b9067fe8 str w8, [sp,#1660]
64: 52933348 mov w8, #0x999a // #39322
68: 72a838a8 movk w8, #0x41c5, lsl #16
6c: b9067be8 str w8, [sp,#1656]
70: 529999a8 mov w8, #0xcccd // #52429
74: 72a855e8 movk w8, #0x42af, lsl #16
78: b90677e8 str w8, [sp,#1652]
7c: bd467be0 ldr s0, [sp,#1656]
80: bd0673e0 str s0, [sp,#1648]
84: bd4677e0 ldr s0, [sp,#1652]
88: bd066fe0 str s0, [sp,#1644]
8c: bd467be0 ldr s0, [sp,#1656]
90: bd4683e1 ldr s1, [sp,#1664]
94: 1e213800 fsub s0, s0, s1
98: bd466fe1 ldr s1, [sp,#1644]
9c: bd467fe2 ldr s2, [sp,#1660]
a0: 1e223821 fsub s1, s1, s2
a4: 1e210800 fmul s0, s0, s1
a8: bd066be0 str s0, [sp,#1640]
ac: bd4673e0 ldr s0, [sp,#1648]
b0: bd4683e1 ldr s1, [sp,#1664]
b4: 1e213800 fsub s0, s0, s1
b8: bd4677e1 ldr s1, [sp,#1652]
bc: bd467fe2 ldr s2, [sp,#1660]
c0: 1e223821 fsub s1, s1, s2
c4: 1e210800 fmul s0, s0, s1
c8: bd0667e0 str s0, [sp,#1636]
cc: bd466be0 ldr s0, [sp,#1640]
d0: bd4667e1 ldr s1, [sp,#1636]
d4: 1e213800 fsub s0, s0, s1
d8: bd0663e0 str s0, [sp,#1632]
dc: bd467be0 ldr s0, [sp,#1656]
e0: bd4683e1 ldr s1, [sp,#1664]
e4: 1e213800 fsub s0, s0, s1
e8: bd466fe2 ldr s2, [sp,#1644]
ec: bd467fe3 ldr s3, [sp,#1660]
f0: 1e233842 fsub s2, s2, s3
f4: bd4673e4 ldr s4, [sp,#1648]
f8: 1e243821 fsub s1, s1, s4
fc: bd4677e4 ldr s4, [sp,#1652]
100: 1e233883 fsub s3, s4, s3
104: 1e230821 fmul s1, s1, s3
108: 1f020400 fmadd s0, s0, s2, s1
10c: bd065fe0 str s0, [sp,#1628]
C++ 标准允许实现以比标称类型更精确的方式计算浮点表达式。它要求实现在将值分配给对象时丢弃多余的精度。
因此,在分配给 d1
和 d2
时,多余的精度将被丢弃,不会对 d1 - d2
产生影响,但是在 (ex - sx) * (py - sy) - (px - sx) * (ey - sy)
中,多余的精度参与评价。请注意,C++ 不仅允许计算中的超精度,而且允许它用于表达式的某些部分而不是其他部分。
特别是,评估像 a*b - c*d
这样的表达式的常用方法是使用乘法指令(不使用超精度)计算 -c*d
,然后使用计算 a*b - c*d
融合乘加指令,有效地使用无限精度进行乘法运算。
您的编译器可能有一个开关来禁用此行为并始终使用标称精度。
分析你给出的反汇编,我想我找到了原因:fused-multiply-add(它是fmadd
指令)。我没有完全分析反汇编,但我认为 t
是这样计算的:
t = fma((ex - sx) * (py - sy), sx - px, ey - sy);
其中fma
的意思是:
fma(a, b, c) = a + b*c;
所以,计算有点不同,因为fma
在评估a + b*c
时只舍入一次(没有fma
,在x=b*c
处舍入,然后a+x
)
您有一些打开此功能的编译器开关,例如 -ffast-math
或 -ffp-contract=on
。
关闭FMA,我想你的问题就解决了。
(现在我看到Eric回答了同样的问题,但我在这里留下我的答案,也许对理解这个问题有一点帮助)
我 运行 遇到了一个非常奇怪的问题,即 ARM64 上的浮点精度。我有一段非常简单的 C++ 代码,看起来像这样:
float sx = some_float_number_1;
float sy = some_float_number_2;
float ex = some_float_number_3;
float ey = some_float_number_4;
float px = ex;
float py = ey;
float d1 = (ex - sx) * (py - sy);
float d2 = (px - sx) * (ey - sy);
float d = d1 - d2;
float t = (ex - sx) * (py - sy) - (px - sx) * (ey - sy);
//32-bit output: d == t == 0
//64-bit output: d == 0, t != 0
理论上,d 应该等于 t 且等于 0,这正是 32 位 ARM 上发生的情况。但由于某些奇怪的原因,t 的输出在 64 位 ARM 上不等于 0,而 d 仍然正确。我从来没有见过这样的错误,所以我不知道是什么导致了这种问题。
编辑:更多信息
- 如果你没有注意到,d 和 t 的输出应该都是 0,因为 (ex - sx) * (py - sy) 等于 (px - sx) * (ey - sy )
- 此问题仅在输入的小数部分不等于 0 时发生。
- 我使用的编译器是 Android NDK r15c 包中包含的 Clang。
EDIT2:这是反汇编
4c: 52933348 mov w8, #0x999a // #39322
50: 72a82828 movk w8, #0x4141, lsl #16
54: b90683e8 str w8, [sp,#1664]
58: 52933348 mov w8, #0x999a // #39322
5c: 72a82728 movk w8, #0x4139, lsl #16
60: b9067fe8 str w8, [sp,#1660]
64: 52933348 mov w8, #0x999a // #39322
68: 72a838a8 movk w8, #0x41c5, lsl #16
6c: b9067be8 str w8, [sp,#1656]
70: 529999a8 mov w8, #0xcccd // #52429
74: 72a855e8 movk w8, #0x42af, lsl #16
78: b90677e8 str w8, [sp,#1652]
7c: bd467be0 ldr s0, [sp,#1656]
80: bd0673e0 str s0, [sp,#1648]
84: bd4677e0 ldr s0, [sp,#1652]
88: bd066fe0 str s0, [sp,#1644]
8c: bd467be0 ldr s0, [sp,#1656]
90: bd4683e1 ldr s1, [sp,#1664]
94: 1e213800 fsub s0, s0, s1
98: bd466fe1 ldr s1, [sp,#1644]
9c: bd467fe2 ldr s2, [sp,#1660]
a0: 1e223821 fsub s1, s1, s2
a4: 1e210800 fmul s0, s0, s1
a8: bd066be0 str s0, [sp,#1640]
ac: bd4673e0 ldr s0, [sp,#1648]
b0: bd4683e1 ldr s1, [sp,#1664]
b4: 1e213800 fsub s0, s0, s1
b8: bd4677e1 ldr s1, [sp,#1652]
bc: bd467fe2 ldr s2, [sp,#1660]
c0: 1e223821 fsub s1, s1, s2
c4: 1e210800 fmul s0, s0, s1
c8: bd0667e0 str s0, [sp,#1636]
cc: bd466be0 ldr s0, [sp,#1640]
d0: bd4667e1 ldr s1, [sp,#1636]
d4: 1e213800 fsub s0, s0, s1
d8: bd0663e0 str s0, [sp,#1632]
dc: bd467be0 ldr s0, [sp,#1656]
e0: bd4683e1 ldr s1, [sp,#1664]
e4: 1e213800 fsub s0, s0, s1
e8: bd466fe2 ldr s2, [sp,#1644]
ec: bd467fe3 ldr s3, [sp,#1660]
f0: 1e233842 fsub s2, s2, s3
f4: bd4673e4 ldr s4, [sp,#1648]
f8: 1e243821 fsub s1, s1, s4
fc: bd4677e4 ldr s4, [sp,#1652]
100: 1e233883 fsub s3, s4, s3
104: 1e230821 fmul s1, s1, s3
108: 1f020400 fmadd s0, s0, s2, s1
10c: bd065fe0 str s0, [sp,#1628]
C++ 标准允许实现以比标称类型更精确的方式计算浮点表达式。它要求实现在将值分配给对象时丢弃多余的精度。
因此,在分配给 d1
和 d2
时,多余的精度将被丢弃,不会对 d1 - d2
产生影响,但是在 (ex - sx) * (py - sy) - (px - sx) * (ey - sy)
中,多余的精度参与评价。请注意,C++ 不仅允许计算中的超精度,而且允许它用于表达式的某些部分而不是其他部分。
特别是,评估像 a*b - c*d
这样的表达式的常用方法是使用乘法指令(不使用超精度)计算 -c*d
,然后使用计算 a*b - c*d
融合乘加指令,有效地使用无限精度进行乘法运算。
您的编译器可能有一个开关来禁用此行为并始终使用标称精度。
分析你给出的反汇编,我想我找到了原因:fused-multiply-add(它是fmadd
指令)。我没有完全分析反汇编,但我认为 t
是这样计算的:
t = fma((ex - sx) * (py - sy), sx - px, ey - sy);
其中fma
的意思是:
fma(a, b, c) = a + b*c;
所以,计算有点不同,因为fma
在评估a + b*c
时只舍入一次(没有fma
,在x=b*c
处舍入,然后a+x
)
您有一些打开此功能的编译器开关,例如 -ffast-math
或 -ffp-contract=on
。
关闭FMA,我想你的问题就解决了。
(现在我看到Eric回答了同样的问题,但我在这里留下我的答案,也许对理解这个问题有一点帮助)