50 실수 — 근사의 수학
먼저 알아야 할 것
돌아보기
8장에서 0.1 + 0.2와 0.3은 마지막 비트 하나(1 ulp) 차이의 이웃이라 했고, 비교는 “충분히 가까운가”로 바꿔야 한다고 했다. 그러면 그 “충분히”의 기준 — 엡실론 — 은 어떻게 정하는가?
답. 값의 크기에 따라 달라져야 한다는 것이 8장 셋째 사건의 교훈이다. 0 근처에서는 눈금이 촘촘하니 작은 고정값이면 되지만, 쯤에서는 눈금 간격 자체가 1을 넘으므로 같은 기준이 무의미해진다. 그래서 실무는 두 가지를 갖춘다 — 절대 오차(0 근처용)와 상대 오차(크기에 비례). 이 장의 시연이 둘을 나란히 보인다.
이 장의 필요성과 맥락
이 장이 끝나면
float 와 double의 선택, 비교의 올바른 방법(엡실론 — 절대와 상대), 그리고 특수값(무한대·NaN)까지. 8장의 세 사건이 실행 결과로 확인된다.이 장에서 답할 질문
- 그러면 돈 계산 같은 곳에는 실수를 쓰면 안 되는가?
50.1 타입 선택과 비교
examples/ch50/eps.c
#include <math.h>
#include <stdio.h>
/* 절대 오차: 0 근처의 비교에 쓴다 */
bool near_abs(double a, double b, double eps)
{
return fabs(a - b) < eps;
}
/* 상대 오차: 값이 커지면 눈금도 커지므로 크기에 비례한 허용치를 쓴다 */
bool near_rel(double a, double b, double rel)
{
double scale = fabs(a) > fabs(b) ? fabs(a) : fabs(b);
return fabs(a - b) <= rel * scale;
}
int main(void)
{
double sum = 0.1 + 0.2;
printf("0.1 + 0.2 == 0.3 ? %d\n", sum == 0.3);
printf("compare with absolute eps? %d\n", near_abs(sum, 0.3, 1e-9));
printf("%.20f\n", sum);
printf("%.20f\n", 0.3);
double big = 1e16;
printf("1e16 + 1 == 1e16 ? %d (gap wider than 1)\n", big + 1 == big);
printf("compare with relative eps? %d\n", near_rel(big + 1, big, 1e-12));
return 0;
}
실행 결과
0.1 + 0.2 == 0.3 ? 0
compare with absolute eps? 1
0.30000000000000004441
0.29999999999999998890
1e16 + 1 == 1e16 ? 1 (gap wider than 1)
compare with relative eps? 1
첫 세 줄이 8장 첫 사건의 실행 확인이다 — ==는 거짓이고, 20자리까지 찍어 보면 두 수가 마지막에서 갈린다. near_abs(절대 오차)를 쓰면 참이 된다.
뒷부분이 셋째 사건(흡수)의 확인이다 — 에 1을 더해도 값이 그대로 다. 그 크기에서는 표현 가능한 이웃 사이의 간격이 이미 1보다 넓기 때문이다(8장의 눈금 계산). 여기서는 절대 오차가 아니라 near_rel(상대 오차)이 옳은 도구다.
실무 수칙으로 줄이면 셋이다. 기본은 double — 8장에서 본 대로 정밀도의 여유가 다르고, float는 메모리·대역폭이 아쉬운 곳에서만 고른다. ==는 쓰지 않는다 — 두 실수가 같은지 묻는 코드는 거의 언제나 의심 대상이다(정수와 비교할 때, 0과 비교할 때 등 예외는 있지만 의식 적으로 판단해야 한다). 허용치는 문제에서 나온다 — 계산의 성격과 값의 크기를 보고 정하는 것이지 마법 상수가 있는 것이 아니다.
50.2 특수값 — 무한대와 NaN
IEEE 754(8장)는 평범한 수 말고도 특별한 값을 정의해 두었다. 무한대 (양·음)는 넘침이나 0으로 나누기의 결과로 나오고, NaN(Not a Number)은 “수가 아님” — 이나 음수의 제곱근 같은 정의 불가 연산의 결과다. 여기서 한 가지를 정확히 짚어야 한다. 정수의 0 나누기는 계약 밖(28장)이고, 부동소수점의 0 나누기는 흔히 “정의되어 무한대가 된다”고들 말한다. 그러나 그것은 C 언어 자체의 약속이 아니다. C 표준의 나눗셈 규칙은 제수가 0이면 동작을 정의하지 않는다 — 정수든 실수든 마찬가지다. 무한대가 나오는 것은 IEC 60559(즉 IEEE 754) 의미를 지원하는 구현에서다. 그런 구현에서는 유한한 0 아닌 값을 0으로 나누면 부호 있는 무한대가 나오고 “0 나누기” 예외가 서며, 은 NaN 쪽이다. 오늘의 주류 컴파일러와 x86-64·AArch64 에서는 그렇게 동작하지만, 표준이 모든 구현에 강제하는 바는 아니다.
플랫폼 노트. 이 구별이 실무에서 갈리는 자리
구현이 이진 부동소수점에 대해 IEC 60559 를 따른다고 선언하는 표시는 __STDC_IEC_60559_BFP__ 다(BFP = binary floating point). 이 매크로가 정의돼 있으면 부속서 F 의 의미가 계약이 되고, 위의 동작을 기대할 수 있다. 십진 부동소수점 쪽은 __STDC_IEC_60559_DFP__ 로 따로 표시한다.
이름이 비슷한 __STDC_IEC_60559_BF16_TYPES__ 와 헷갈리지 않아야 한다 — 그것은 bfloat16 타입을 제공하는지에 관한 별개의 기능 표시이지, 일반 이진 부동소수점의 연산 의미와는 상관이 없다. 덧붙여, 타입의 형식이 IEC 60559 와 같다는 것과 연산 의미까지 부속서 F 를 따른다는 것은 서로 다른 질문이다 — 앞은 __STDC_IEC_60559_TYPES__ 쪽 이야기이고, 여기서 필요한 것은 뒤쪽이다. 그렇지 않은 환경 — 일부 임베디드 툴체인, 그리고 -ffast-math 같이 부속서 F 의미를 일부러 끄는 빌드 — 에서는 무한대·NaN 이 나온다는 보장이 사라진다. 이식 가능한 코드는 0 제수를 먼저 걸러 낸다.
50.2.1 비트를 직접 열어 본다
8장에서 배치를 보았으니, 이제 실제 값의 비트를 찍어 확인한다. 표현을 옮길 때는 공용체가 아니라 memcpy를 쓴다 — 48장의 규칙대로 “값을 옮기는” 것이 가장 안전한 통로이고, 컴파일러는 이 복사를 대개 없애 버린다.
examples/ch50/bits.c
/* 실수의 속을 열어 본다 — 부호·지수·가수가 실제로 어떻게 들어 있는가.
타입 퍼닝은 공용체가 아니라 memcpy 로 한다(37·48장의 규칙). */
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
static uint64_t bits_of(double d)
{
uint64_t u;
memcpy(&u, &d, sizeof u); /* 표현을 그대로 옮긴다 */
return u;
}
static uint32_t bits_of_f(float f)
{
uint32_t u;
memcpy(&u, &f, sizeof u);
return u;
}
/* double: 부호 1 + 지수 11 + 가수 52 */
static void dump(const char *label, double d)
{
uint64_t u = bits_of(d);
unsigned sign = (unsigned)(u >> 63);
unsigned expo = (unsigned)((u >> 52) & 0x7FFu);
uint64_t frac = u & 0xFFFFFFFFFFFFFu;
printf("%-12s %016" PRIx64 " sign %u exponent %4u(=%+5d) significand %013" PRIx64,
label, u, sign, expo,
expo == 0 ? -1022 : (int)expo - 1023, frac);
if (expo == 0x7FF) printf(" <- %s", frac ? "NaN" : "infinity");
else if (expo == 0) printf(" <- %s", frac ? "subnormal" : "zero");
printf("\n");
}
int main(void)
{
puts("-- the bits of a double (sign 1 + exponent 11 + significand 52) --");
dump("1.0", 1.0);
dump("-1.0", -1.0);
dump("0.5", 0.5);
dump("2.0", 2.0);
dump("0.1", 0.1);
dump("0.3", 0.3);
dump("0.1+0.2", 0.1 + 0.2);
dump("0.0", 0.0);
dump("-0.0", -0.0);
dump("inf", INFINITY);
dump("NaN", NAN);
puts("\n-- 0.1 + 0.2 and 0.3 have different bits --");
printf("0.1+0.2 = %.20f\n", 0.1 + 0.2);
printf("0.3 = %.20f\n", 0.3);
printf("bits of the difference: %016" PRIx64 " vs %016" PRIx64 " (the last bit)\n",
bits_of(0.1 + 0.2), bits_of(0.3));
puts("\n-- one step (ULP) up from 1.0 raises the last bit of the significand --");
double one = 1.0;
double next = nextafter(1.0, 2.0);
printf("1.0 %016" PRIx64 "\n", bits_of(one));
printf("next value %016" PRIx64 " difference %.17g\n", bits_of(next), next - one);
printf("same as DBL_EPSILON? %s\n",
(next - one) == 0x1p-52 ? "yes" : "no");
puts("\n-- a float uses the same structure, narrower (sign 1 + exponent 8 + significand 23) --");
float f = 0.1f;
uint32_t fu = bits_of_f(f);
printf("0.1f %08" PRIx32 " sign %u exponent %3u(=%+4d) significand %06" PRIx32 "\n",
fu, fu >> 31, (fu >> 23) & 0xFFu, (int)((fu >> 23) & 0xFFu) - 127,
fu & 0x7FFFFFu);
printf("printing (double)0.1f again gives %.17g - the narrowing to float leaves a trace\n",
(double)f);
puts("\n-- below the smallest normal number you get subnormals --");
double small = 0x1p-1022; /* 가장 작은 정규수 */
dump("2^-1022", small);
dump("half", small / 2); /* 비정규수 */
dump("2^-1074", 0x1p-1074); /* 가장 작은 비정규수 */
dump("half again", 0x1p-1074 / 2); /* 0 으로 가라앉는다 */
return 0;
}
실행 결과
-- the bits of a double (sign 1 + exponent 11 + significand 52) --
1.0 3ff0000000000000 sign 0 exponent 1023(= +0) significand 0000000000000
-1.0 bff0000000000000 sign 1 exponent 1023(= +0) significand 0000000000000
0.5 3fe0000000000000 sign 0 exponent 1022(= -1) significand 0000000000000
2.0 4000000000000000 sign 0 exponent 1024(= +1) significand 0000000000000
0.1 3fb999999999999a sign 0 exponent 1019(= -4) significand 999999999999a
0.3 3fd3333333333333 sign 0 exponent 1021(= -2) significand 3333333333333
0.1+0.2 3fd3333333333334 sign 0 exponent 1021(= -2) significand 3333333333334
0.0 0000000000000000 sign 0 exponent 0(=-1022) significand 0000000000000 <- zero
-0.0 8000000000000000 sign 1 exponent 0(=-1022) significand 0000000000000 <- zero
inf 7ff0000000000000 sign 0 exponent 2047(=+1024) significand 0000000000000 <- infinity
NaN 7ff8000000000000 sign 0 exponent 2047(=+1024) significand 8000000000000 <- NaN
-- 0.1 + 0.2 and 0.3 have different bits --
0.1+0.2 = 0.30000000000000004441
0.3 = 0.29999999999999998890
bits of the difference: 3fd3333333333334 vs 3fd3333333333333 (the last bit)
-- one step (ULP) up from 1.0 raises the last bit of the significand --
1.0 3ff0000000000000
next value 3ff0000000000001 difference 2.2204460492503131e-16
same as DBL_EPSILON? yes
-- a float uses the same structure, narrower (sign 1 + exponent 8 + significand 23) --
0.1f 3dcccccd sign 0 exponent 123(= -4) significand 4ccccd
printing (double)0.1f again gives 0.10000000149011612 - the narrowing to float leaves a trace
-- below the smallest normal number you get subnormals --
2^-1022 0010000000000000 sign 0 exponent 1(=-1022) significand 0000000000000
half 0008000000000000 sign 0 exponent 0(=-1022) significand 8000000000000 <- subnormal
2^-1074 0000000000000001 sign 0 exponent 0(=-1022) significand 0000000000001 <- subnormal
half again 0000000000000000 sign 0 exponent 0(=-1022) significand 0000000000000 <- zero
출력에서 다섯 가지를 짚어 둔다.
첫째, 1.0은 놀랍도록 단정하다. 지수 자리가 1023(바이어스 그대로라 실제 지수 0)이고 가수는 전부 0이다 — 숨은 비트 1만으로 이 된다. 2.0은 지수만 하나 오르고, 0.5는 하나 내린다. 부호 비트만 뒤집으면 -1.0이 된다.
둘째, 0.1은 무한소수의 잘린 자국을 그대로 보여 준다. 가수가 999999999999a로 끝나는데, 마지막 자리의 a가 바로 반올림의 흔적이다. 8장 수학 박스의 “2진법으로 안 떨어진다”가 비트로 드러난 자리다.
셋째, 0.1 + 0.2와 0.3은 마지막 한 비트가 다르다. 두 값의 비트열은 ...3334와 ...3333 — 정확히 1만큼 어긋나 있다. 그래서 == 비교가 거짓이 되고, 그래서 이 장이 허용치를 이야기한다.
넷째, 한 눈금(ULP)의 정체가 보인다. 1.0의 비트에 정수 1을 더한 것이 바로 다음 실수이고, 그 차이가 DBL_EPSILON()이다. “1.0 근처에서 구분 가능한 가장 작은 차이”라는 말이 비트 하나라는 사실로 확인된다.
다섯째, 바닥에서 비정규수가 나타난다. 가장 작은 정규수를 반으로 나누면 지수가 더 내려갈 수 없어 대신 가수의 앞자리부터 0이 채워진다 — 지수 자리가 전부 0인 그 상태가 비정규수다. 정밀도를 조금씩 잃으면서 0까지 촘촘히 내려가고, 마지막 한 비트마저 사라지면 0이 된다. 갑자기 0으로 떨어지지 않고 서서히 잦아드는 이 설계를 점진적 언더플로(gradual underflow)라 부른다.
플랫폼 노트. 비정규수는 느릴 수 있다
비정규수 연산은 하드웨어에 따라 정규수보다 훨씬 느리다(수십 배가 되는 기계도 있다). 그래서 신호 처리나 게임 엔진에서는 비정규수를 0으로 밀어 버리는 모드(flush-to-zero)를 켜기도 한다 — 정확성을 조금 포기하고 최악의 지연을 없애는 거래다. 표준 C에는 이 모드를 켜는 이식 가능한 방법이 없다(컴파일러 옵션이나 플랫폼 API를 쓴다).NaN에는 유명한 성질이 하나 있다 — 자기 자신과도 같지 않다. x != x가 참이면 x는 NaN이라는 뜻이고, 이것이 NaN을 판별하는 고전적 관용구다(오늘은 isnan()을 쓴다). “같음”이라는 관계의 기본 성질이 깨지는 값이라, 정렬이나 검색 알고리즘에 NaN이 섞여 들면 기묘한 일이 벌어진다 — 실수 데이터를 다룰 때 NaN 검사를 경계에 두는 것이 관행인 이유다.
실제 사례. 0.1초의 누적 — 패트리어트 미사일 사고
8장에서 배운 “작은 어긋남이 쌓인다”가 실제로 사람의 목숨과 이어진 사건이 있다. 1991년 걸프전에서 패트리어트 방공 시스템이 날아오는 미사일을 요격하지 못해 28명이 사망한 사고인데, 원인 분석의 핵심이 부동소수점 오차였다. 시스템은 시간을 0.1초 단위로 세었는데 — 8장에서 본 대로 0.1은 2진법으로 무한소수라 담을 때마다 미세한 오차가 생긴다. 그 시스템이 24비트 그릇을 쓴 탓에 오차가 상대적으로 컸고, 재부팅 없이 100시간을 연속 가동하자 누적 오차가 약 0.34초에 이르렀다. 그 0.34초 동안 목표는 500미터 넘게 이동하고, 추적 창은 엉뚱한 하늘을 보게 된 것이다. “근사는 성실하지만 무해하지는 않다” — 8장의 교훈이 가장 무겁게 확인된 사례다.문. 그러면 돈 계산 같은 곳에는 실수를 쓰면 안 되는가?
답. 쓰지 않는 것이 정석이다 — 8장에서 배운 고정소수점의 자리가 바로 여기다. 금액을 원 단위 실수로 다루면 0.1원 단위의 어긋남이 쌓여 장부가 맞지 않게 되므로, 전 단위 정수로 계산하고 표시할 때만 소수점을 찍는 것이 금융 소프트웨어의 관행이다. 규칙으로 줄이면 — 정확한 십진 값이 중요한 곳에는 정수(고정소수점), 물리량과 과학 계산에는 부동 소수점. 도구를 문제에 맞추는 것이고, 그 판단의 근거가 8장과 이 장에서 배운 표현의 성격이다.
근사의 세계를 C에서 다루는 법을 익혔다. 다음 장은 이 부의 중심 주제 — 계산이 실패할 수 있다는 사실을 프로그램이 어떻게 다루는가, 오류와 계약의 이야기다.