부동소수점 수의 비교

얼마전에 재민군이 질문했던 내용인데, 태준형의 도움으로 해결했던 문제입니다.

Question

먼저 다음과 같은 코드가 있습니다.

#include 
int main()
{
int i = 10;
double val = i*12.34;
if(val == i * 12.34){
printf("Equaln");
} else {
printf("Not Equaln");
}
}

결과는 무엇일까요?

물론 답은 이 프로그램을 컴파일하고 실행하는 환경에 따라 다르다입니다. 즉, 어셈블리 코드를 생성하는 컴파일러와 이 코드를 실행하는 아키텍쳐에 따라 다를 것이라고 예상할 수 있습니다. 그런데, 신기한 것은 같은 아키텍쳐(Pentium 3)와 같은 컴파일러(gcc 3.4)를 사용하더라도 OS에 따라(Linux와 FreeBSD) 다른 결과가 나온다는 것입니다. 자세히 얘기하면, Linux에서는 “Not Equal”이 출력되고, FreeBSD에서는 “Equal”이 출력됩니다. 어째서 그럴까요?

Answer

설마, 어셈블리 코드가 다르게 나오는 걸까요? OS가 달라서 그럴지도, 컴파일러 설정이 달라서 그럴지도 모르지 않습니까? 아니면 컴파일러 버그? 그런데, 그렇지 않습니다. 부동소수점 연산과 비교를 수행하는 부분의 코드는 다음과 같이 똑같습니다.

movl    $10, -4(%ebp)
fildl   -4(%ebp)
fldl    .LC0
fmulp   %st, %st(1)
fstpl   -16(%ebp)
fildl   -4(%ebp)
fldl    .LC0
fmulp   %st, %st(1)
fldl    -16(%ebp)
fxch    %st(1)
fucompp
fnstsw  %ax
andb    $69, %ah
cmpb    $64, %ah
je  .L3
jmp .L2

간단히 설명하면, 10과 12.34를 FPU stack에 넣고 곱셈을 수행한 후, (FPU stack에 들어있는) 결과를 다시 메모리에 저장합니다. 두번째 연산을 위해 역시 10과 12.34를 FPU stack에 넣고 곱셈을 수행한 후, 메모리에 저장했던 결과도 다시 FPU stack에 불러들여, 두 값의 비교를 수행하는 코드입니다.

두 값의 차이는 곱셈의 결과가 메모리로 잠시 다녀왔는가의 여부 뿐이죠. OS에 따른 결과의 차이도 바로 여기서 발생합니다.

The difference between Linux and FreeBSD is that in Linux, the FPU is operated by default in “extended precision” mode, where 80-bit internal registers are used for floating point numbers. In FreeBSD, the default is to use “double precision” mode, where 64-bit precision is retained.

Excerpt from Resolution of Differences Between X86 Linux/glibc Floating-Point to Integer Conversion Behavior Relative to Other Platforms

즉, 387 내에서 부동소수점 수는 80bit의 정밀도를 가질 수 있는데, 메모리 상에서는 IEEE 754대로 64bit의 정밀도를 가지므로 일종의 rounding error가 일어난다는 거죠. 그리고 FreeBSD에서는 80bit 정밀도 기능을 OS 수준에서 꺼버리지만 Linux는 그렇지 않기 때문에 차이가 나는 것이구요.

More Fun?

참고로, 제 amd64 머신에서는 “Equal”이 출력되고 다음과 같은 코드를 생성합니다.

movl    $10, -4(%rbp)
cvtsi2sd    -4(%rbp), %xmm1
movlpd  .LC0(%rip), %xmm0
mulsd   %xmm1, %xmm0
movsd   %xmm0, -16(%rbp)
cvtsi2sd    -4(%rbp), %xmm1
movlpd  .LC0(%rip), %xmm0
mulsd   %xmm0, %xmm1
movlpd  -16(%rbp), %xmm0
ucomisd %xmm0, %xmm1
jp  .L5
je  .L3

즉, SSE instruction을 사용하죠. gcc는 x86-64 (aka. amd64) 아키텍쳐에서는 기본적으로 SSE를 사용한다고 합니다.

sse Use scalar floating point instructions present in the SSE
instruction set. This instruction set is supported by Pentium3
and newer chips, in the AMD line by Athlon-4, Athlon-xp and
Athlon-mp chips. The earlier version of SSE instruction set
supports only single precision arithmetics, thus the double and
extended precision arithmetics is still done using 387. Later
version, present only in Pentium4 and the future AMD x86-64
chips supports double precision arithmetics too.
For i387 you need to use -march=cpu-type, -msse or -msse2
switches to enable SSE extensions and make this option effec-
tive. For x86-64 compiler, these extensions are enabled by
default.
The resulting code should be considerably faster in the major-
ity of cases and avoid the numerical instability problems of
387 code, but may break some existing code that expects tempo-
raries to be 80bit.
This is the default choice for the x86-64 compiler.

Pentium 4이상의 머신에서 다음과 같이 컴파일 하면 역시 “Equal”이 출력되는 것을 확인할 수 있습니다. (Pentium 3에서도 컴파일은 되지만, Illegal Instruction 예외가 발생합니다.)

gcc -msse -msse2 -mfpmath=sse test_double.c

Conclusion

문제의 원인은 x86 FPU가 80bit extension을 사용하고 있고, FreeBSD에서는 이 기능을 꺼버리고 Linux는 그대로 두는 것에 있습니다.

어떻게 보면, x86의 FPU가 사용하는 80bit extension이 표준(IEEE754)을 잘 지키지 않은 것이 문제라거나 Linux가 문제의 소지가 있는 CPU 기능을 켜놓은 것이 문제라고 볼 수도 있습니다. 하지만, 모두들 아시다시피 부동소수점 수의 비교 연산은 보장받기 어려운 것이므로, 이런 간단한 비교 연산의 정확도를 버리고 연산의 정밀도를 택한 trade-off는 정당하다고도 생각됩니다.

결국, 부동소수점 수의 비교 연산은 아무리 간단하더라도 rounding error와 같은 문제를 무시해서는 안된다는 것이죠.

위에서 인용했던 글에서 읽어보라고 하고 있는 David Goldberg의 What Every Computer Scientist Should Know About Floating-Point Arithmetic를 읽어볼 생각입니다.

부동소수점 수의 비교

Question

Answer

More Fun?

Conclusion

관련