Странная проблема с точностью с плавающей точкой на 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 не равен 0 на 64-битной ARM, в то время как d все еще корректен. Я никогда не видел подобной ошибки, поэтому понятия не имею, что могло вызвать такую проблему.
РЕДАКТИРОВАТЬ: немного больше информации
- В случае, если вы не заметили, выходные данные как d, так и t должны быть равны 0, поскольку (ex-sx) * (py-sy) равно (px-sx) * (ey-sy)
- Эта проблема возникает только тогда, когда дробная часть ввода не равна 0.
- Я использую компилятор Clang, включенный в комплект Android NDK r15c.
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, и я думаю, что ваша проблема будет решена.
(Теперь я вижу, что Эрик ответил на то же самое, но я оставляю свой ответ здесь, возможно, это поможет немного понять проблему)