오차 없는 변환

1. 개요

오차 없는 변환 Error-free transformation (EFT)
정의	$\mathrm{fl}(a\circ b)$ 와 그 반올림 오차를 둘 다 부동소수점 수로 정확히 얻는 것
TwoSum	$a+b=s+e$ — 6 flops, 분기 없음 (묄러 1965 · 커누스)
FastTwoSum	3 flops, 단 $\lvert a\rvert\ge\lvert b\rvert$ 필요 (데커 1971)
TwoProduct	$ab=p+e$ — FMA 로 2 연산, 없으면 분할로 17 flops
전제	최근접 반올림, 오버플로·(곱셈의) 언더플로 없음
쓰임	보정 합·내적·호너, 더블-더블 산술, 정확 기하 술어

반올림 오차를 버리지 말고 주워 담아라. 놀랍게도 그게 정확히 담긴다.

오차 없는 변환(error-free transformation, EFT)은 부동소수점 연산 한 번의 결과 $\mathrm{fl}(a \circ b)$ 와 그 연산이 만든 반올림 오차 $e$ 를 둘 다 부동소수점 수로, 그것도 정확히 얻어 내는 짧은 연산열이다. 덧셈과 곱셈에 대해

$$a + b = s + e, \qquad a\times b = p + e, \qquad s = \mathrm{fl}(a+b),\ p = \mathrm{fl}(ab)$$

가 실수 등식으로 성립한다. 오차 $e$ 가 근사치가 아니라 정확한 값이라는 점이 핵심이다. IEEE 754 의 최근접 반올림 아래에서 두 부동소수점 수의 합의 반올림 오차는 언제나 부동소수점 수로 표현 가능하고(오버플로가 없는 한), 곱의 오차도 마찬가지다(언더플로가 없는 한).

이 사실 하나로 할 수 있는 일이 의외로 많다. 오차들을 따로 모아 두 번째 누산기로 굴리면 배정밀도 하드웨어로 두 배 정밀도 계산의 정확도를 얻는다 — 호너 방법의 보정 호너, 클렌쇼 점화식의 보정판, 오기타·럼프·오이시(2005)의 정확한 합과 내적이 모두 그렇다. 오차를 끝까지 전개하면 정확 술어의 셔척 확장 산술이 되고, 두 수의 쌍을 한 수처럼 다루면 더블-더블 산술이 된다. 혼합 정밀도 문서의 카한 보상 합도 사실 이 계열의 원조 격이다.

2. TwoSum 과 FastTwoSum

FastTwoSum(데커 1971)은 연산 세 번이다:

s = a + b
e = b - (s - a)        # |a| >= |b| 일 때 정확

$s - a$ 가 정확히 계산된다는 것(스터베노 보조정리 계열의 논리)이 핵심인데, 이것은 $\lvert a\rvert \ge \lvert b\rvert$(엄밀하게는 $a$ 의 지수가 $b$ 의 지수 이상)일 때만 보장된다. 순서를 어기면 틀린다 — $a = 1$, $b = 10^{100}$ 으로 넣으면 $s = 10^{100}$, $e = 0$ 이 나와 1 이 통째로 사라진다.

TwoSum(묄러 1965, 커누스)¹은 순서를 몰라도 되는 대신 연산이 여섯 번이다:

s  = a + b
bb = s - a
e  = (a - (s - bb)) + (b - bb)

같은 $a = 1$, $b = 10^{100}$ 에서 $e = 1$ 을 정확히 돌려준다. 실무에서는 비교·분기 한 번을 아끼는 대신 연산 세 번을 더 쓰는 TwoSum 이 오히려 빠른 경우가 많다. 현대 CPU 에서 예측 실패한 분기가 덧셈 몇 번보다 비싸고, 분기가 없으면 SIMD 로 벡터화하기도 쉽기 때문이다.

3. TwoProduct — FMA 가 있으면 한 줄

곱의 오차는 FMA(fused multiply-add, $\mathrm{fl}(ab + c)$ 를 반올림 한 번으로 계산하는 명령)가 있으면 그냥

p = a * b
e = fma(a, b, -p)      # ab - p 를 정확히

이다. $ab - p$ 는 정확히 부동소수점 수이고 FMA 는 그것을 반올림 한 번으로 계산하므로 반올림할 것이 없다. 2008년판 IEEE 754 가 FMA 를 필수 연산으로 넣었고 x86(Haswell 이후)·ARMv8·GPU 가 모두 하드웨어로 지원하니 지금은 이것이 표준이다.

FMA 가 없던 시절에는 데커의 TwoProduct 를 썼다. 먼저 벨트캄프 분할로 53비트 가수를 위·아래 26비트 정도의 두 조각 $a = a_h + a_l$ 로 쪼갠다(배정밀도면 상수 $2^{27} + 1 = 134217729$ 를 곱해 쓰는 4 연산)². 조각끼리의 곱 $a_hb_h$, $a_hb_l$, $a_lb_h$, $a_lb_l$ 은 비트 수가 53 을 안 넘어 정확하므로, 그것들로 $ab - p$ 를 정확히 조립할 수 있다. 분할 두 번과 조립을 합쳐 17 flops. 단 $\lvert a\rvert$ 가 $10^{300}$ 근처면 분할 상수를 곱하는 순간 넘치므로 큰 수는 미리 스케일해야 한다.

무작위 $(a, b)$ 10만 쌍씩(지수 $2^{\pm30}$ 범위)에 대해 TwoSum 과 분할판 TwoProduct 의 결과를 유리수로 검산하면 $s + e = a + b$, $p + e = ab$ 가 모든 쌍에서 비트 단위로 성립한다.

4. 보정 내적 — 정밀도 두 배의 정확도

오기타·럼프·오이시(2005)의 Dot2 는 내적 $\sum x_iy_i$ 를 이렇게 계산한다: 각 곱을 TwoProduct 로 $(h, r)$ 로 쪼개고, 주 누산기 $p$ 에 $h$ 를 TwoSum 으로 더해 거기서 나온 오차 $q$ 와 곱의 오차 $r$ 을 보조 누산기 $s$ 에 그냥(보통 부동소수점으로) 모은다. 끝에 $p + s$ 를 반환한다. 오차 상한은

$$\lvert \mathrm{Dot2}(x,y) - x^\mathsf{T}y\rvert \;\le\; u\,\lvert x^\mathsf{T}y\rvert + \gamma_n^2\,\lvert x\rvert^\mathsf{T}\lvert y\rvert, \qquad \mathrm{cond} = \frac{2\,\lvert x\rvert^\mathsf{T}\lvert y\rvert}{\lvert x^\mathsf{T}y\rvert}$$

이고($u = 2^{-53}$, $\gamma_n \approx nu$), 상대오차로 쓰면 $u + \tfrac12\gamma_n^2\,\mathrm{cond}$ 다. 이것은 두 배 정밀도($u^2$)로 계산한 뒤 배정밀도로 반올림한 것과 같은 형태다. 일부러 조건수를 키운 길이 100 의 내적(오기타 등의 생성법을 따라 만듦)으로 재면:

조건수	순진한 내적 상대오차	Dot2 상대오차
$6.1\times10^{8}$	$1.7\times10^{-8}$	$3.4\times10^{-17}$
$1.0\times10^{17}$	6.6	$2.9\times10^{-15}$
$7.2\times10^{24}$	$5.4\times10^{8}$	$1.9\times10^{-7}$
$4.8\times10^{32}$	$2.3\times10^{16}$	2.0

순진한 내적은 오차가 대략 $u\cdot\mathrm{cond}$ 라 조건수 $10^{16}$ 에서 이미 부호조차 못 믿는다. Dot2 는 $u^2\cdot\mathrm{cond}$ 라 $10^{16}$ 에서도 열다섯 자리, $10^{24}$ 에서도 일곱 자리가 살아 있고, $1/u^2 \approx 10^{32}$ 을 넘으면 비로소 무너진다. 정확히 후진 오차 해석이 예측하는 그림이다. flops 는 원소당 순진한 내적의 2 에 비해 FMA 가 있으면 10 안팎, 없으면 25 안팎이지만, 메모리 대역폭에 묶인 큰 내적에서는 실측 시간 차이가 그보다 훨씬 작다.

이것을 $K$ 번 겹친 SumK·DotK 는 $K$ 배 정밀도의 정확도를 낸다. 또 결과를 「충실 반올림」(정답에 인접한 두 부동소수점 수 중 하나)까지 보장하는 정확 합산 알고리즘들도 같은 저자들이 이어서 냈다. 파이썬의 math.fsum 은 셔척식 부분합 목록으로 정확히 반올림된 합을 돌려준다.

5. 더블-더블 — 두 수를 한 수처럼

수 하나를 겹치지 않는 부동소수점 수 두 개의 합 $x = x_h + x_l$($\lvert x_l\rvert \le \tfrac12\mathrm{ulp}(x_h)$)로 들고 다니면 가수가 약 106비트인 「더블-더블」 수가 된다. 덧셈은 TwoSum 과 FastTwoSum 몇 번, 곱셈은 TwoProduct 와 교차항 정리로 구현하고, 한 연산에 배정밀도 연산 십수~수십 번이 든다. 히다·리·베일리(2001)의 QD 라이브러리(더블-더블·쿼드-더블)가 대표적이고, 소프트웨어 4배 정밀도(__float128)보다 대개 빠르다. 다만 지수 범위는 배정밀도 그대로이고, $x_h$ 가 비정규화 근처로 가면 $x_l$ 이 언더플로해 정밀도가 소리 없이 떨어진다. 혼합 정밀도 반복 정제에서 잔차만 높은 정밀도로 계산하고 싶을 때 이 방식이 자주 쓰인다.

6. 깨지는 조건 — 컴파일러가 제일 위험하다

오차 없는 변환은 IEEE 754 를 글자 그대로 지킨다는 전제 위에 서 있다. 그 전제를 깨는 것들:

• -ffast-math / /fp:fast. 컴파일러가 결합법칙을 가정하면 (a - (s - bb)) 같은 식을 대수적으로 0 으로 접어 버린다. 보정항이 통째로 사라지고, 코드는 경고 없이 순진한 계산으로 퇴화한다. EFT 가 든 파일은 따로 엄격한 부동소수점 옵션으로 컴파일해야 한다.³
• x87 확장 정밀도. 중간값이 80비트 레지스터에 머물다가 나중에 64비트로 반올림되면(이중 반올림) 「$s$ 는 $a + b$ 를 한 번 반올림한 값」이라는 전제가 깨진다. SSE2 이후의 x86-64 에서는 기본이 아니지만, 32비트 레거시 빌드에서 재현되는 고전 버그다.
• 자동 FMA 축약. 컴파일러가 a*b - p 를 멋대로 FMA 로 바꾸는 것은 TwoProduct 에는 오히려 이득이지만, 분할판 TwoProduct 의 중간식을 FMA 로 바꾸면 정확성 증명이 깨질 수 있다(-ffp-contract=off).
• 오버플로·언더플로. TwoSum 은 $s$ 가 넘치면 $e$ 가 NaN 이 되고, TwoProduct 의 $e$ 는 곱의 오차가 비정규화 범위로 내려가면 정확하지 않다.
• 방향 반올림. 구간 산술처럼 반올림 모드를 위·아래로 바꿔 둔 상태에서는 합의 오차가 부동소수점 수라는 보장이 없다.

7. 관련 문서

• 부동소수점 연산 · 혼합 정밀도 · 후진 오차 해석 · 조건수 · IEEE 754
• 호너 방법 · 클렌쇼 점화식 · 정확 술어 · 구간 해석
• 더블-더블 연산 · 재현 가능한 합산 · 임의 정밀도 연산

8. Footnotes

1. 묄러(1965)의 원 논문 제목은 「부동소수점 덧셈의 준(準)배정밀도」다. 하드웨어에 배정밀도가 없던 시절, 단정밀도 두 개로 배정밀도를 흉내 내려던 것이 출발점이었다. 60년이 지나 지금은 GPU 의 단정밀도·반정밀도로 같은 짓을 다시 하고 있다. ↩

2. 53비트를 26+26 으로 나누면 1비트가 남는데, 그 1비트는 아래 조각의 부호가 가져간다. 아래 조각이 음수일 수 있어서 26비트로 충분하다는 게 이 트릭의 백미다. 처음 보면 다들 비트 수를 세다가 한 번씩 멘붕한다. ↩

3. 「빨라진다길래 켰더니 정확도 보정 코드가 0 이 됐어요」는 수치 라이브러리 이슈 트래커의 국룰 제목이다. 컴파일러 입장에서는 $(a+b)-a-b = 0$ 이라는 중학교 수학을 적용했을 뿐이라 억울할 것이다. ↩