[완료] 찾기어려운 버그~ floating point 뺄셈.
찾기 어려웠던 버그 이야기가 나와서,
아래 문제는 현재 진행중인 문제입니다.
큰 코드인데, 코드를 바꾸면 재현이 안되는 관계로 일부만 올립니다.
double a = preac1->GetJumpRate(system_temperature);
double b = preac2->GetJumpRate(system_temperature);
double d = a - b;
printf("a = %llx\n", *(unsigned long long *)&a);
printf("b = %llx\n", *(unsigned long long *)&b);
printf("d = %llx\n", *(unsigned long long *)&d);결과는 아래와 같습니다.
a = 4248900326281980 b = 4248900326281980 d = bee9b80000000000
우선 바이너리 시퀀스가 똑같은 a,b 두 값을 뺀 결과가 exactly 0가 안나오고 있고,
floating point의 rounding off 오차를 감안하더라도
10진수로 -0.000012264 정도로 double precision에서 꽤 큰값입니다.
결국 뺄때의 a,b값과 출력할때의 a,b값이 다르다는 결론밖에 안나옵니다.
volatile double b
로 선언하면 정확히 0가 나옵니다.
intel cpu, ubuntu linux, gcc 4.3.3 환경에서 -O1으로 컴파일했고,
멀티쓰레드환경은 아닙니다.
optimize 안하면 정확히 0 나옵니다.
메모리 침범 문제는 어느정도 배제됐고 optimize관련 컴파일러 버그를 의심하고 있습니다.
volatile이 있는 경우와 없는 경우의 위 c코드에 해당하는 asm 코드는 아래와 같이 딱 두줄만 다릅니다.
%%% 를 기준으로 앞부분이 오동작하는경우, 뒤부분이 volatile을 넣어서 정상동작하는 경우입니다.
fstl이 fstpl로 바뀌고, fldl -32(%ebp)라인이 추가되는 군요..
맨앞에 첨부한 코드에서 프린트하기 직전까지의 asm코드 입니다.
.LEHE104:
.loc 1 268 0
movl 8(%ebp), %edx
fldl 220(%edx)
fstpl 4(%esp)
movl %esi, (%esp)
.LEHB105:
call _ZN8Particle11GetJumpRateEd
fstpl -24(%ebp)
.LVL559:
.loc 1 269 0
movl 8(%ebp), %ecx
fldl 220(%ecx)
fstpl 4(%esp)
movl %ebx, (%esp)
call _ZN8Particle11GetJumpRateEd
fstl -32(%ebp) %%%%%%%%%%%%%%%%%%% fstpl -32(%ebp)
.LVL560:
.loc 1 271 0
%%%%%%%%%%%%%%%%%% fldl -32(%ebp)
fsubrl -24(%ebp)
fstpl -40(%ebp)중요한건 아닌데 개인적으로 너무 궁금하군요..
===
답글 감사합니다.
궁금증을 해소하는데 많은 도움이 됐습니다.
아래 jick님이 알려주신 내용을 바탕으로
문제가 된 부분을 재현하는 코드를 구성했습니다.
컴파일 옵션에 따라 차이가 0이 되기도 하고 아니기도 하네요.
#include <cmath>
#include <stdio.h>
int main()
{
unsigned long long temp = 0xc009db657e605192LL;
double d2 = *(double*)&temp;
double nu = 5.34513e12;
double a = nu * std::exp(d2);
double b = nu * std::exp(d2);
double d1 = a - b;
printf("a = %llx\n", *(unsigned long long *)&a);
printf("b = %llx\n", *(unsigned long long *)&b);
printf("d1 = %llx\n", *(unsigned long long *)&d1);
printf("a = %12.9e\n", a);
printf("d1 = %12.9e\n", d1);
return 0;
}
추측입니다만 CPU 버그가 아닐까요?
위의 두 코드의 차이는 실제로 없다고 봐도 무방합니다.
메모리에 썼다가 다시 읽는 것 외에는 차이가 없으니까요.
코드를 바꾸면 재현되지 않는다는 점에서 더욱 CPU 버그를 의심하게 됩니다.
Intel에 정식으로 문의해 보시는게 어떨까 싶네요.
전 초보입니다만 CPU
전 초보입니다만
CPU 내부에서 floating 연산을 할 때
80bit 레지스터를 이용해서 처리를 한다는 이야기를 들었는데
그 변환에 관련된 문제는 아닐까요?
gcc버그로 대충
gcc버그로 대충 결론나는 분위기네요. auditory님이 0022??
x86 CPU의 구조적인(?) 문제입니다
재현되는 코드가 없으니 100% 그 문제라고 확신할 수는 없습니다만...
x86의 floating point register는 64bit double이 아니라 80bit extended double입니다. 그런데 C에서 double이라고 하면 64bit double을 의미하게 되죠.
따라서, 변수값을 floating point register에 담았다가, 다시 메모리에 (double 형태로) 저장했다가, 다시 레지스터로 읽어들이면 원래와 값이 달라질 수 있습니다.
gcc manual을 뒤져보니 -ffloat-store를 사용하면 변수값을 무조건 메모리에 썼다가 다시 가져오기 때문에 (물론 성능이 나빠지겠지만) 이런 문제를 예방할 수 있다고 하는군요. 혹은 387 floating point instruction set 대신에 SSE 명령어를 쓰도록 -msse -mfpmath=sse를 써도 된다는 것 같은데 테스트를 해보지 않아서 100% 장담은 못하겠습니다.
* 이 문제에 대해 좀 더 자세한 설명을 원하시면, Computer Systems: A Programmer's Perspective (Bryant & O'Hallaron) 이라는 책에서 section 2.4.6 Floating Point in C를 찾아보세요. 문제를 일으키는 예제 코드도 나옵니다.
* 키즈에서 메모리 문제라는데 손모가지 걸겠다는 분 이제 손모가지 자르려나... (먼산)
이것은
"How not to program in C++(2003)" 이라는 책에 소개되어 있는 문제입니다.
이 책의 Program 40. "Unbelievable Accuracy" 에 보면 정확히 같은 문제점과 해결책이 소개되어 있네요.
참고하세요.
아 참 그리고
키즈에 보니까-_- "원래 floating point는 오차가 생기는 거야"라고 훈계 아닌 훈계를 하는 사람들이 있는데...
오차가 생긴다는 의미는 "a에서 b를 뺐을 때 그 결과값이 수학적으로 a-b와 다를 수 있다."라는 뜻입니다. CPU floating point unit이 무슨 난수제너레이터도 아니고, 같은 값 a에서 같은 값 b를 뺐는데 그 값이 0이 됐다 안됐다 할 수가 없습니다. 그렇게 나오면 CPU가 망가졌거나 무슨 다른 사연이 있는 겁니다.
"a-b는 0인데 b-a가 0이 안나와요!" 하는 데다 대고 "원래 부동소숫점은 정확하지 않은 거 몰라?" 하던 사람들은 좀 반성하고 IEEE floating point standard 관련 문서 한번 정독하시기 바랍니다. 별 괴상한 NaN 같은 특수케이스 빼고 나면, 핵심 내용은 별로 복잡하지 않습니다.
80-bit에서 64-bit로 변환되면서 생기는 오차로 보기에는 좀 의심스러운 부분이 있습니다.
위의 뺄셈 계산 결과인 64비트 오차 값을 80비트로 변환한 뒤
floating point register stack내에서 덧셈하여 다시 빼 보았습니다만...
오차값이 그대로 계산되어 나옵니다.
제가 사용한 소스는 다음과 같습니다.
(GNU C로는 변수 다루기가 힘들어서 Visual C++용으로 만든 소스입니다)
#include <stdio.h> int main(void) { double a, b, c; double adj; union conv { double d; unsigned long l[2]; }; ((union conv *)&a)->l[0] = 0x26281980UL; ((union conv *)&a)->l[1] = 0x42489003UL; ((union conv *)&b)->l[0] = 0x26281980UL; ((union conv *)&b)->l[1] = 0x42489003UL; ((union conv *)&adj)->l[0] = 0x00000000UL; ((union conv *)&adj)->l[1] = 0xbee9b800UL; __asm finit; __asm fld a; // st(0)에 64-bit 변수 a의 값을 80-bit로 변환하여 push __asm fld adj; // st(0)에 64-bit 변수 adj의 값을 80-bit로 변환하여 push __asm faddp st(1), st(0); // 두 80-bit 값을 더하고 결과값을 st(0)에 저장 __asm fld b; // st(0)에 64-bit 변수 b의 값을 80-bit로 변환하여 push __asm fsubr st(0), st(1); // st(0)에서 st(1)을 뺀 뒤 결과를 st(0)에 저장 __asm fstp c; // 결과를 64-bit 변수 c에 저장 printf("RESULT: %f\n", c); printf("RESULT: 0x%08lx%08lx\n", ((union conv *)&c)->l[1], ((union conv *)&c)->l[0]); return 0; }결과는 0xbee9b80000000000이 나옵니다. 즉, 변환 오차로 보았던 값이 정상적으로 계산된 겁니다.
혹시 제가 잘못 알고 있는 부분이 있거나 실험 방법에 문제가 있다면 지적 부탁드립니다.
재현 코드를 보니 뭐가 잘못된 건지 알겠습니다.
이미 계산 결과로 나온 오차 값에는 아무 의미가 없었군요.
계산 중간에 사용된 값에서 mentissa의 64bit ~ 80bit 사이 구간이
오차로 나오는 것이라서 exponent의 값에 곱을 하면 의외로 큰 숫자도
나올 수 있다는 사실을 깨닫지 못했습니다.
좋은 정보 감사합니다.
댓글 달기