부동소수점과 오차

0.1 + 0.2 ≠ 0.3: 수치해석의 출발점

원인은 간단합니다. 컴퓨터는 수를 2진법으로 저장하는데, 10진법의 0.1은 2진법에서 끝나지 않는 소수입니다. 1/3이 10진법에서 $0.333\ldots$로 끝나지 않는 것과 같습니다. 0.1에 2를 거듭 곱하며 정수 부분을 떼어 내면 2진 자리가 하나씩 나옵니다.

이 현상은 프로그래밍 언어의 결함이 아니라 유한한 비트로 실수를 다루는 모든 계산의 성질입니다. 그리고 수치해석이 다루는 오차 가운데 하나일 뿐입니다. 실제 계산의 답이 참값과 다른 이유는 크게 네 가지입니다.

절단 오차는 단계 수를 늘리거나 간격 $h$를 줄이면 작아집니다. 3장부터 16장까지 이 책의 대부분은 절단 오차를 $O(h^p)$로 재고 줄이는 이야기입니다. 반올림 오차는 반대로, 단계를 늘릴수록 쌓이고 $h$를 줄일수록 커지기도 합니다(9장의 V자 그래프). 두 오차의 균형을 잡으려면 먼저 반올림 오차의 크기를 정확히 알아야 합니다. MathBook의 컴퓨터는 어떻게 틀리는가에서 맛본 내용을 여기서 엄밀하게 다시 세웁니다.

정수와 고정소수점의 한계

가장 단순한 표현은 정수입니다. 32비트 부호 있는 정수(2의 보수)는 $-2^{31}$부터 $2^{31} - 1 = 2\,147\,483\,647$까지를 간격 1로 빠짐없이 담습니다. 이 범위 안의 덧셈·곱셈은 정확합니다. 그러나 범위를 넘으면 경고 없이 값이 돌아갑니다. $2\,147\,483\,647 + 1$은 32비트 연산에서 $-2\,147\,483\,648$이 됩니다(자바스크립트로 (2147483647 + 1) | 0).

소수를 담으려면 고정소수점(fixed point) 표현을 씁니다. 32비트 중 아래 16비트를 소수 부분으로 정하면(Q16.16 형식) 저장된 정수 $N$이 값 $N/2^{16}$을 뜻합니다. 간격은 $2^{-16} \approx 1.5\times10^{-5}$로 일정하고, 범위는 $\pm 32768$입니다. 신호 처리 칩이나 부동소수점 장치가 없는 마이크로컨트롤러는 지금도 이 방식을 씁니다. 하지만 간격이 절대적으로 고정되어 있다는 것이 치명적입니다.

• $0.00001$을 저장하면 가장 가까운 값은 $1/65536 = 0.0000153$이라 상대 오차가 50%를 넘는다.
• $30000$을 저장하면 상대 오차는 $10^{-10}$ 수준으로 지나치게 정확하다.
• $10^{-31}$과 $10^{30}$을 함께 담으려면 소수점 앞뒤로 각각 100비트가 넘게 필요하다.

과학자들이 오래전부터 쓰는 해법은 과학적 표기법입니다. $9.1\times10^{-31}$처럼 유효 숫자 몇 자리와 지수를 따로 적으면, 크기와 상관없이 상대 정밀도가 일정합니다. 컴퓨터의 부동소수점은 이것을 2진법으로 옮긴 것입니다. 소수점이 지수에 따라 "떠다닌다"(floating)고 해서 붙은 이름입니다.

IEEE 754 배정밀도: 부호·지수·가수

$\beta = 2$이면 정규화된 수의 첫 자리 $d_0$은 언제나 1입니다. 그래서 저장할 필요가 없습니다. 이 "숨은 비트" 덕분에 52비트로 53자리의 정밀도를 얻습니다. IEEE 754의 배정밀도(binary64, double)가 바로 $\mathbb F(2, 53, -1022, 1023)$입니다.

배정밀도의 53비트는 10진법으로 $53\log_{10}2 \approx 15.95$자리입니다. "double은 십진 15–16자리"라는 말이 여기서 나옵니다. 단정밀도는 약 7자리입니다. 기계 학습 가속기가 즐겨 쓰는 16비트 형식은 2자리에서 3자리 남짓이라, 2장의 혼합 정밀도 기법과 함께 써야 합니다.

아래 편집기에서 비트를 직접 눌러 보세요. 64칸 각각이 값에 어떻게 기여하는지, 지수 비트를 바꾸면 값이 2배씩 뛰고 소수 마지막 비트를 바꾸면 값이 ulp 한 칸만 움직이는 것을 볼 수 있습니다.

수직선 위의 부동소수점 수: 간격은 크기에 비례한다

부동소수점 수는 수직선 위에 고르게 놓여 있지 않습니다. 지수가 같은 수들, 곧 구간 $[2^e, 2^{e+1})$ 안의 수들은 가수의 마지막 자리 한 칸씩 떨어져 있습니다. 지수가 하나 커지면 같은 개수의 수가 두 배 넓은 구간에 퍼집니다.

배정밀도에서 1 다음 수는 $1 + 2^{-52}$, 1000 근처의 간격은 $1.14\times10^{-13}$, $10^{16}$ 근처의 간격은 2입니다. 그래서 $10^{16} + 1$은 $10^{16}$으로 반올림되어 사라집니다. 정수를 빠짐없이 담는 한계는 $2^{53} = 9\,007\,199\,254\,740\,992$이고, $2^{53} + 1$은 double로 표현할 수 없습니다. 자바스크립트의 Number.MAX_SAFE_INTEGER가 $2^{53} - 1$인 이유입니다.

아래 시뮬레이터는 정밀도 $p$와 지수 범위를 작게 줄인 "장난감 부동소수점 체계"를 수직선에 그립니다. 64비트 체계도 정확히 같은 모양이고, 단지 눈금이 $2^{52}$배 촘촘할 뿐입니다.

기계 엡실론, 반올림 단위, ulp

IEEE 754의 기본 반올림은 가장 가까운 수로, 동률이면 마지막 비트가 0인(짝수) 쪽으로(round to nearest, ties to even)입니다. 동률을 늘 위로 올리면 긴 계산에서 오차가 한쪽으로 쏠리기 때문에, 짝수 쪽을 골라 평균적으로 상쇄되게 한 것입니다. 이 반올림에서 오차 상한은 간격의 절반입니다.

이 정리는 수 하나를 저장할 때의 이야기입니다. IEEE 754는 여기에 결정적인 약속을 더합니다. $+, -, \times, \div, \sqrt{\ }$의 결과는 정확한 결과를 반올림한 값이어야 합니다(correctly rounded). 그러므로 부동소수점 수 $x, y$에 대해

기계 엡실론은 실험으로 찾을 수 있습니다. 1에 $e$를 더해 1보다 커지는지 확인하며 $e$를 반씩 줄이면 됩니다.

// 기계 엡실론 찾기: 1 + e/2 가 1과 구별되지 않을 때까지 e를 반씩 줄인다
e ← 1
while fl(1 + e/2) > 1 do
    e ← e / 2
end
return e          // double: 52번 반복 후 e = 2^(-52)

오버플로·언더플로·비정규 수, Inf와 NaN

지수의 범위가 유한하므로 표현할 수 있는 크기에도 위아래 끝이 있습니다.

비정규 수는 0과 $x_{\min}$ 사이를 간격 $2^{-1074}$로 고르게 채웁니다(그림의 주황 눈금). 이것을 점진적 언더플로(gradual underflow)라 합니다. 덕분에 "$x \ne y$이면 $x - y \ne 0$"이 언제나 성립합니다. 대가는 정밀도입니다. 비정규 수에서는 앞자리 0들이 가수를 차지하므로 상대 오차 상한 $u$가 성립하지 않고, $2^{-1074}$ 근처의 수는 유효 숫자가 한 자리뿐입니다.

연산이 수학적으로 정의되지 않으면 결과는 NaN(Not a Number)입니다. 특수 값의 규칙은 대부분 극한의 상식을 따릅니다.

NaN과 Inf는 "계산을 멈추지 않고 문제를 끝까지 운반"하도록 설계되었습니다. 계산 도중 한 번 NaN이 생기면 이후 거의 모든 연산 결과가 NaN이 되므로, 마지막에 한 번만 검사하면 됩니다. 반대로 말하면 NaN이 나왔을 때는 처음 생긴 곳을 찾아야 합니다.

빗변 문제의 처방은 척도 조정(scaling)입니다. 큰 쪽 $m = \max(|x|, |y|)$으로 나누어 계산하면 중간값이 1 근처에 머뭅니다.

// 오버플로·언더플로 없는 빗변 sqrt(x^2 + y^2)
m ← max(|x|, |y|)
if m = 0 then return 0
r ← min(|x|, |y|) / m          // 0 ≤ r ≤ 1
return m · sqrt(1 + r·r)

절대 오차, 상대 오차, 유효 숫자

오차의 크기를 말하는 방법은 두 가지입니다. 0.001의 오차는 지구–달 거리에서는 놀라운 정확도이고, 반도체 게이트 길이(수십 나노미터)를 미터로 잴 때는 쓸모없는 값입니다. 크기에 견준 오차가 필요합니다.

참값이 0에 가까우면 상대 오차는 쓸 수 없으므로 절대 오차를 씁니다. 실무의 멈춤 조건이 흔히 $|\hat x - x| \le \mathrm{atol} + \mathrm{rtol}\cdot|x|$처럼 둘을 섞는 이유입니다(3장, 11장). 상대 오차의 중요한 성질은 곱셈과 나눗셈에서 거의 더해진다는 것입니다. $\hat x = x(1 + \delta_1)$, $\hat y = y(1 + \delta_2)$이면 $\hat x\hat y = xy(1 + \delta_1 + \delta_2 + \delta_1\delta_2)$. 곱셈·나눗셈은 상대 오차를 키우지 않습니다. 문제는 덧셈과 뺄셈입니다.

자릿수 상쇄: 비슷한 두 수의 뺄셈

뺄셈 자체는 정확합니다. 아주 가까운 두 부동소수점 수의 차는 반올림 없이 정확히 계산됩니다(스테르베즈 보조정리(Sterbenz lemma): $y/2 \le x \le 2y$이면 $x - y$는 부동소수점 수다). 문제는 뺄셈에 들어가는 값이 이미 오차를 품고 있다는 것입니다. 앞자리가 같은 두 수를 빼면 같은 앞자리가 사라지고, 오차가 섞인 뒷자리만 남아 앞으로 올라옵니다. 이것을 자릿수 상쇄(catastrophic cancellation)라 합니다.

근의 공식에서는 $-b = 10^8$과 $\sqrt{b^2 - 4ac} = 99999999.99999998$이 앞 16자리가 거의 같습니다. $\sqrt{\ }$의 결과에 있던 $10^{-16}$ 수준의 상대 오차가 $10^8/10^{-8}$배 증폭되어 결과를 망칩니다. 처방은 빼지 않는 식으로 바꾸는 것입니다. 두 근의 곱이 $c/a$라는 비에트 공식을 쓰면, 상쇄가 없는 큰 근을 먼저 구하고 작은 근을 나눗셈으로 얻을 수 있습니다.

// 상쇄 없는 이차방정식 ax^2 + bx + c = 0 (실근, a ≠ 0)
d ← sqrt(b·b − 4·a·c)
q ← −(b + sign(b)·d) / 2     // b와 sign(b)·d는 부호가 같다 → 더해도 상쇄 없음
x1 ← q / a                   // 크기가 큰 근
x2 ← c / q                   // 작은 근: 비에트 공식 x1·x2 = c/a

표의 둘째 열은 $b$가 100배 커질 때마다 오차가 약 $10^4$배 커집니다. 증폭 인자 $\approx b^2$과 맞습니다. 같은 일이 $1 - \cos x$에서도 일어납니다. $x = 10^{-8}$이면 $\cos x = 1 - 5\times10^{-17}$인데 이 값은 1과 $u$보다 가까워 $\mathrm{fl}(\cos x) = 1$이고, 결과는 0입니다. 참값은 $5\times10^{-17}$입니다. 삼각 항등식 $1 - \cos x = 2\sin^2(x/2)$는 뺄셈이 없어 모든 $x$에서 정확합니다.

상쇄가 어디서 일어나는지 한 단계씩 보면 더 분명합니다. 아래 표는 근의 공식의 중간값을 17자리로 보여 주고, 뺄셈에서 사라지는 같은 앞자리를 밑줄로 표시합니다.

오차의 누적과 합산 순서, 카한 합산

$n$개의 수를 차례대로 더하는 재귀 합산 $s_1 = x_1$, $s_k = \mathrm{fl}(s_{k-1} + x_k)$의 오차를 표준 모형으로 추적해 봅시다. 각 덧셈에서 $s_k = (s_{k-1} + x_k)(1 + \delta_k)$, $|\delta_k| \le u$입니다.

이 정리에서 두 가지 교훈이 나옵니다. 첫째, 먼저 더한 항일수록 더 많은 반올림을 겪으므로, 크기가 작은 항부터 더하면 큰 $m$이 작은 $|x_i|$에 곱해져 실제 오차가 줄어듭니다. 둘째, 상한이 $\sum|x_i|$에 비례하므로, 부호가 섞여 $|S| \ll \sum|x_i|$인 합(교대 급수)은 상대 오차가 클 수 있습니다. 뺄셈 상쇄의 합산판입니다.

순서를 바꾸는 것보다 강력한 방법은 잃어버린 부분을 따로 모으는 것입니다. 카한(Kahan)이 1965년에 발표한 보정 합산은 덧셈마다 반올림으로 떨어져 나간 아랫자리를 보정 변수 $c$에 정확히 회수해 다음 항에 되돌려 줍니다.

// 카한 보정 합산 (compensated summation)
s ← 0;  c ← 0                // c: 지금까지 잃어버린 아랫자리(부호 반대)
for i = 1 to n do
    y ← x[i] − c             // 지난번에 잃은 만큼 보정한 항
    t ← s + y                // 큰 s에 작은 y를 더하면 y의 아랫자리가 떨어진다
    c ← (t − s) − y          // (t − s)는 실제로 더해진 부분 → 빼면 떨어진 부분
    s ← t
end
return s

보정 합산의 둘째 줄 $c \leftarrow (t - s) - y$는 수학적으로는 0입니다. 컴파일러가 "대수적으로 단순화"하면 알고리즘이 무너지므로, 실무에서는 공격적인 부동소수점 최적화 옵션을 끄고 씁니다. 비용은 덧셈 네 번으로 재귀 합산의 4배입니다. NumPy의 sum은 짝짓기 합산을, 파이썬의 math.fsum은 이보다 더 정확한 방법을 씁니다.

실제 사고: 시간 누적 오차와 변환 오버플로

반올림 오차는 교과서의 장난감이 아닙니다. 공식 조사 보고서로 원인이 확인된 두 사건을 봅시다.

계산해 봅시다. 0.1을 소수점 아래 23비트에서 자르면 $\lfloor 0.1\times2^{23}\rfloor/2^{23} = 0.09999990463\ldots$이고, 한 번의 곱마다 오차는 $9.54\times10^{-8}$초입니다. 100시간은 0.1초 단위로 $3.6\times10^{6}$틱이므로 누적 오차는 $3.6\times10^6\times9.54\times10^{-8} = 0.343$초입니다. 스커드는 초속 약 1.7 km로 날아오므로, 0.34초면 레이더가 미사일을 찾는 위치가 500 m 넘게 어긋납니다. 오차 하나하나는 $10^{-7}$이었지만 오래 쌓이고, 빠른 물체의 위치로 곱해져 치명적이 되었습니다.

두 사건의 교훈은 다릅니다. 패트리어트는 반올림 오차의 누적, 아리안 5는 표현 범위(오버플로)의 문제였습니다. 둘 다 계산 하나하나는 명세대로 정확했고, 문제는 "이 오차·이 범위가 이 사용 조건에서 괜찮은가"를 묻지 않은 데 있었습니다. 수치해석의 질문이 바로 그것입니다.

이 도구가 쓰이는 곳

반올림 단위 $u$와 상쇄·누적의 원리는 이 책의 거의 모든 장에서 다시 나옵니다. 몇 가지 예를 미리 짚어 둡니다.