54 실수 — 근사의 수학

54.1 타입 선택과 비교#

#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

실행 결과의 첫 네 줄이 7장 첫 사건의 실행 확인이다 — ==는 거짓이고, 20자리까지 찍어 보면 두 수가 마지막에서 갈린다. near_abs(절대 오차)를 쓰면 참이 된다.

뒷부분이 셋째 사건(흡수)의 확인이다 — ${10}^{16}$에 1을 더해도 값이 그대로다. 그 크기에서는 표현 가능한 이웃 사이의 간격이 이미 1보다 넓기 때문이다(7장의 눈금 계산). 여기서는 절대 오차가 아니라 near_rel(상대 오차)이 옳은 도구다.

실무 수칙으로 줄이면 셋이다. 기본은 double — 7장에서 본 대로 정밀도의 여유가 다르고, float는 메모리·대역폭이 아쉬운 곳에서만 고른다. ==는 쓰지 않는다 — 두 실수가 같은지 묻는 코드는 거의 언제나 의심 대상이다(정수와 비교할 때, 0과 비교할 때 등 예외는 있지만 의식적으로 판단해야 한다). 허용치는 문제에서 나온다 — 계산의 성격과 값의 크기를 보고 정하는 것이지 마법 상수가 있는 것이 아니다.

54.2 특수값 — 무한대와 NaN#

IEEE 754(7장)는 평범한 수 말고도 특별한 값을 정의해 두었다. 무한대 (양·음)는 넘침이나 0으로 나누기의 결과로 나오고, NaN(Not a Number)은 “수가 아님” — $\frac{0}{0}$이나 음수의 제곱근 같은 정의 불가 연산의 결과다. 여기서 0 으로 나누면 무엇이 나오는지를 정확히 해 두어야 한다. 정수의 0 나누기는 계약 밖(29장)이고, 부동소수점의 0 나누기는 흔히 “정의되어 무한대가 된다”고들 말한다. 그러나 그것은 C 언어 자체의 약속이 아니다. C 표준의 나눗셈 규칙은 제수가 0이면 동작을 정의하지 않는다 — 정수든 실수든 마찬가지다. 무한대가 나오는 것은 IEC 60559(즉 IEEE 754) 의미를 지원하는 구현에서다. 그런 구현에서는 유한한 0 아닌 값을 0으로 나누면 부호 있는 무한대가 나오고 “0 나누기” 예외가 일어나며, $\frac{0}{0}$은 NaN 쪽이다. 오늘의 주류 컴파일러와 x86-64·AArch64 에서는 그렇게 동작하지만, 표준이 모든 구현에 강제하는 바는 아니다.

54.2.1 비트를 직접 열어 본다#

7장에서 배치를 보았으니, 이제 실제 값의 비트를 찍어 확인한다. 표현을 옮길 때는 공용체가 아니라 memcpy를 쓴다 — 49장의 규칙대로 “값을 옮기는” 것이 가장 안전한 통로이고, 컴파일러는 이 복사를 대개 없애 버린다.

/* 실수의 속을 열어 본다 — 부호·지수·가수가 실제로 어떻게 들어 있는가.
   타입 퍼닝은 공용체가 아니라 memcpy 로 한다(49 · 56장의 규칙). */
#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만으로 $1.0\times 2^0$이 된다. 2.0은 지수만 하나 오르고, 0.5는 하나 내린다. 부호 비트만 뒤집으면 -1.0이 된다.

둘째, 0.1은 무한소수의 잘린 자국을 그대로 보여 준다. 가수가 999999999999a로 끝나는데, 마지막 자리의 a가 바로 반올림의 흔적이다. 7장 수학 박스의 “2진법으로 안 떨어진다”가 비트로 드러난 자리다.

셋째, 0.1 + 0.2와 0.3은 마지막 한 비트가 다르다. 두 값의 비트열은 ...3334와 ...3333 — 정확히 1만큼 어긋나 있다. 그래서 == 비교가 거짓이 되고, 그래서 이 장이 허용치를 이야기한다.

넷째, 한 눈금(ULP)의 정체가 보인다. 1.0의 비트에 정수 1을 더한 것이 바로 다음 실수이고, 그 차이가 DBL_EPSILON($2^{\left\{-52\right\}}$)이다. “1.0 근처에서 구분 가능한 가장 작은 차이”라는 말이 비트 하나라는 사실로 확인된다.

다섯째, 바닥에서 비정규수가 나타난다. 가장 작은 정규수를 반으로 나누면 지수가 더 내려갈 수 없어 대신 가수의 앞자리부터 0이 채워진다 — 지수 자리가 전부 0인 그 상태가 비정규수다. 정밀도를 조금씩 잃으면서 0까지 촘촘히 내려가고, 마지막 한 비트마저 사라지면 0이 된다. 갑자기 0으로 떨어지지 않고 서서히 잦아드는 이 설계를 점진적 언더플로(gradual underflow)라 부른다. 이 바닥이 어떻게 생겼는지, 그리고 그 위에서 무엇이 무너졌는지는 이 장 뒤쪽의 두 절에서 따로 연다.

NaN에는 유명한 성질이 하나 있다 — 자기 자신과도 같지 않다. x != x가 참이면 x는 NaN이라는 뜻이고, 이것이 NaN을 판별하는 고전적 관용구다(오늘은 isnan()을 쓴다). “같음”이라는 관계의 기본 성질이 깨지는 값이라, 정렬이나 검색 알고리즘에 NaN이 섞여 들면 기묘한 일이 벌어진다 — 실수 데이터를 다룰 때 NaN 검사를 경계에 두는 것이 관행인 이유다.

54.3 0 근처의 바닥 — 비정규수#

「비트를 직접 열어 본다」의 다섯째 항목에서 스치고 지나간 비정규수(subnormal number)를 이제 제대로 연다. 옛 문헌에서는 denormal(비정규화 수)이라고도 부르는데, 같은 것이다. 이 절은 세 가지를 묻는다. 왜 이런 것이 필요한가. 비트로는 어떻게 생겼는가. 그리고 그 대가는 무엇인가. 마지막 질문의 답이 다음 절의 보안 사고들로 곧장 이어진다.

54.3.1 정규수만 있다면 생기는 구멍#

정규수는 가수가 언제나 1.xxx 꼴이다. 그러니 double 의 가장 작은 양의 정규수는 지수가 바닥(−1022)이고 가수가 1.000…0 인 값, 곧 $2^{-1022}$ 다. <float.h> 는 이 값을 DBL_MIN(약 $2.2\times {10}^{-308}$)이라 부른다. 그 바로 위 이웃은 $2^{-1022}+2^{-1074}$ 이다 — 이 근처의 눈금 간격은 $2^{-1074}$ 다.

그런데 아래쪽을 보자. 정규수만으로는 0 과 DBL_MIN 사이에 아무것도 없다. 가수의 맨 앞이 1이어야 한다는 규칙 때문에 DBL_MIN 보다 작은 정규수는 만들 수 없기 때문이다. 그 빈 구간의 너비는 $2^{-1022}$ — 바로 위 눈금 간격의 $2^{52}$ 배, 약 4500조 배다. 눈금이 촘촘하게 이어지다가 0 앞에서 갑자기 낭떠러지가 나오는 꼴이다.

그림 54.1 — 정규수만 있으면 0 앞에 구멍이 생긴다. 비정규수는 그 구멍을 같은 간격의 눈금으로 메운다.

구멍이 왜 문제인가. 두 정규수 $𝑥=1.5\times 2^{-1022}$ 와 $𝑦=2^{-1022}$ 를 보자. 둘은 분명히 다르다. 그런데 차 $𝑥-𝑦=0.5\times 2^{-1022}$ 는 정확히 구멍 안에 떨어진다. 구멍이 비어 있다면 이 차는 가장 가까운 표현 가능한 값, 곧 0 으로 반올림될 수밖에 없다. 다음 코드를 보라.

double inverse_gap(double x, double y)
{
    if (x != y)
        return 1.0 / (x - y);   /* 다르다고 확인했으니 0 으로 나눌 일은 없다? */
    return 0.0;
}

「두 수가 다르면 그 차는 0이 아니다」는 너무 당연해서 아무도 의심하지 않는 전제다. 구멍이 있는 세계에서는 이 전제가 깨진다 — 위 함수는 x != y 를 확인하고도 0으로 나눈다. 비정규수가 있는 IEEE 754 에서는 이 전제가 언제나 성립한다. 두 유한한 수의 차가 0이면 두 수는 같고, 그 역도 참이다. 비정규수가 지키는 것은 몇몇 아주 작은 수가 아니라 이런 당연한 추론의 안전이다.

54.3.2 비정규수에서는 해석 규칙이 두 가지 바뀐다#

비정규수는 지수 칸이 전부 0 인 값이다. 이때 해석 규칙이 두 가지 동시에 바뀐다.

1. 숨은 비트가 1이 아니라 0이다. 가수는 1.xxx 가 아니라 0.xxx 로 읽는다. 앞자리가 0이 될 수 있으니 DBL_MIN 보다 작은 값을 만들 수 있다.
2. 실제 지수는 −1023 이 아니라 −1022 로 고정된다. 정규수의 규칙(저장된 지수 − 1023)을 그대로 적용하면 $0-1023=-1023$ 이어야 할 것 같지만, 비정규수의 지수는 $1-1023=-1022$ 다. 지수 칸이 1인 가장 작은 정규수와 같은 지수를 쓴다.

그래서 비정규수의 값은 이렇게 복원된다. 정규수의 식과 나란히 놓고 보면 차이가 정확히 두 곳이다.

그림 54.2 — 정규·비정규 경계 양쪽의 비트. 지수 칸은 1에서 0으로 내려가지만 실제 지수는 그대로다.

둘째 규칙이 핵심이다. 처음 배우는 사람이 가장 자주 틀리는 자리이기도 하다 — 다음 절의 사고 가운데 자바의 변환 사고가 바로 이 −1022 를 −1023 으로 다룬 코드에서 났다.

float 과 double 두 형식의 경계값을 표로 모은다. 이름은 모두 <float.h> 에 있다. DBL_TRUE_MIN 처럼 「가장 작은 비정규수」를 가리키는 이름은 C11 에서 들어왔다.

54.3.3 비정규수의 대가 — 바닥에서는 정밀도가 샌다#

비정규수는 구멍을 메우지만 공짜가 아니다. 가수의 앞자리가 0으로 채워질수록 쓸 수 있는 유효 비트가 줄어든다. 정규수는 언제나 53비트의 정밀도를 갖지만, 비정규수는 가수 칸에서 맨 앞의 1부터 끝까지만 남는다. 가장 작은 비정규수 DBL_TRUE_MIN 에는 유효 비트가 단 하나뿐이다. 실행해 확인한다.

/* 0 근처의 바닥 — 비정규수의 비트와 성질.
   지수 칸이 전부 0이 되면 규칙이 두 가지 바뀐다: 숨은 비트가 0이 되고,
   실제 지수는 0 - 1023 이 아니라 1 - 1023 = -1022 로 고정된다.
   비트는 memcpy 로 옮겨 읽는다(bits.c 와 같은 방법). */
#include <float.h>
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static uint64_t bits_of(double d)
{
    uint64_t u;
    memcpy(&u, &d, sizeof u);
    return u;
}

static const char *kind(double d)
{
    switch (fpclassify(d)) {
    case FP_ZERO:      return "zero";
    case FP_SUBNORMAL: return "subnormal";
    case FP_NORMAL:    return "normal";
    case FP_INFINITE:  return "infinite";
    default:           return "NaN";
    }
}

/* 유효 비트 수: 정규수는 숨은 비트를 더해 언제나 53.
   비정규수는 가수 칸에서 맨 앞의 1 부터 끝까지만 남는다. */
static int significant_bits(double d)
{
    uint64_t u = bits_of(d);
    unsigned expo = (unsigned)((u >> 52) & 0x7FFu);
    uint64_t frac = u & 0xFFFFFFFFFFFFFu;
    if (expo != 0) return 53;
    int n = 0;
    while (frac) { n++; frac >>= 1; }
    return n;
}

/* 비트 전체·지수 칸·실제 지수·숨은 비트를 한 줄로 */
static void dump(const char *label, double d)
{
    uint64_t u = bits_of(d);
    unsigned expo = (unsigned)((u >> 52) & 0x7FFu);
    printf("%-15s %016" PRIx64 "  exp %4u (2^%-5d) hidden %d  %-9s %2d bits\n",
           label, u, expo, expo == 0 ? -1022 : (int)expo - 1023,
           expo == 0 ? 0 : 1, kind(d), significant_bits(d));
}

int main(void)
{
    puts("-- the floor of double: both sides of the border --");
    double min_normal = DBL_MIN;                     /* 2^-1022 */
    double max_sub    = nextafter(DBL_MIN, 0.0);     /* 바로 아래 이웃 */
    double min_sub    = DBL_TRUE_MIN;                /* 2^-1074 (C11) */
    dump("above DBL_MIN", nextafter(DBL_MIN, 1.0));
    dump("DBL_MIN", min_normal);
    dump("max subnormal", max_sub);
    dump("DBL_TRUE_MIN", min_sub);
    printf("DBL_MIN      = %.17g\n", DBL_MIN);
    printf("DBL_TRUE_MIN = %.17g\n", DBL_TRUE_MIN);
    printf("issubnormal(DBL_MIN) %d, issubnormal(largest subnormal) %d   (C23)\n",
           issubnormal(min_normal) != 0, issubnormal(max_sub) != 0);

    puts("\n-- the step is the same on both sides of the border --");
    printf("DBL_MIN - largest subnormal         == DBL_TRUE_MIN ? %s\n",
           min_normal - max_sub == DBL_TRUE_MIN ? "yes" : "no");
    printf("next above DBL_MIN - DBL_MIN        == DBL_TRUE_MIN ? %s\n",
           nextafter(DBL_MIN, 1.0) - min_normal == DBL_TRUE_MIN ? "yes" : "no");

    puts("\n-- x != y, so x - y is not zero (gradual underflow) --");
    double x = 1.5 * DBL_MIN, y = DBL_MIN;
    double diff = x - y;
    printf("x == y ? %s\n", x == y ? "yes" : "no");
    dump("x - y", diff);
    printf("x - y == 0 ? %s\n", diff == 0.0 ? "yes" : "no");

    puts("\n-- the price: 0.1 pushed down by 2^-k and lifted back --");
    for (int k = 1018; k <= 1074; k += 8) {
        double low = ldexp(0.1, -k);
        printf("k = %4d  %-9s %2d bits  back up: %.17g\n",
               k, kind(low), significant_bits(low), ldexp(low, k));
    }

    puts("\n-- the string that once hung PHP and Java --");
    const char *texts[] = { "2.2250738585072011e-308", "2.2250738585072012e-308" };
    for (int i = 0; i < 2; i++) {
        double d = strtod(texts[i], nullptr);
        printf("%s -> %016" PRIx64 "  %s\n", texts[i], bits_of(d), kind(d));
    }

    puts("\n-- one subnormal row poisons the running sum --");
    /* 열 개의 행 가운데 4번 행만 아주 작은 값을 낸다. 합계 d 가 한번
       비정규수가 되면, 그 뒤의 덧셈은 모두 비정규수를 피연산자로 받는다. */
    double d = 0.0;
    int slow = 0;
    for (int row = 0; row < 10; row++) {
        double v = (row == 3) ? 1e-310 : 0.0;
        if (issubnormal(d) || issubnormal(v)) slow++;
        d += v;
    }
    printf("sum %s, additions that touched a subnormal: %d of 10\n", kind(d), slow);
    return 0;
}

-- the floor of double: both sides of the border --
above DBL_MIN   0010000000000001  exp    1 (2^-1022) hidden 1  normal    53 bits
DBL_MIN         0010000000000000  exp    1 (2^-1022) hidden 1  normal    53 bits
max subnormal   000fffffffffffff  exp    0 (2^-1022) hidden 0  subnormal 52 bits
DBL_TRUE_MIN    0000000000000001  exp    0 (2^-1022) hidden 0  subnormal  1 bits
DBL_MIN      = 2.2250738585072014e-308
DBL_TRUE_MIN = 4.9406564584124654e-324
issubnormal(DBL_MIN) 0, issubnormal(largest subnormal) 1   (C23)

-- the step is the same on both sides of the border --
DBL_MIN - largest subnormal         == DBL_TRUE_MIN ? yes
next above DBL_MIN - DBL_MIN        == DBL_TRUE_MIN ? yes

-- x != y, so x - y is not zero (gradual underflow) --
x == y ? no
x - y           0008000000000000  exp    0 (2^-1022) hidden 0  subnormal 52 bits
x - y == 0 ? no

-- the price: 0.1 pushed down by 2^-k and lifted back --
k = 1018  normal    53 bits  back up: 0.10000000000000001
k = 1026  subnormal 45 bits  back up: 0.10000000000000142
k = 1034  subnormal 37 bits  back up: 0.1000000000003638
k = 1042  subnormal 29 bits  back up: 0.10000000009313226
k = 1050  subnormal 21 bits  back up: 0.10000002384185791
k = 1058  subnormal 13 bits  back up: 0.100006103515625
k = 1066  subnormal  5 bits  back up: 0.1015625
k = 1074  zero       0 bits  back up: 0

-- the string that once hung PHP and Java --
2.2250738585072011e-308 -> 000fffffffffffff  subnormal
2.2250738585072012e-308 -> 0010000000000000  normal

-- one subnormal row poisons the running sum --
sum subnormal, additions that touched a subnormal: 7 of 10

출력을 위에서부터 읽는다.

첫 묶음은 그림의 비트 그대로다. DBL_MIN 은 지수 칸이 1, 가장 큰 비정규수와 DBL_TRUE_MIN 은 0이다. 그런데 셋 다 실제 지수는 2^-1022 로 같고, 바뀌는 것은 숨은 비트뿐이다. 유효 비트는 53 → 52 → 1 로 줄어든다. C23 에 새로 들어온 issubnormal 이 경계 양쪽을 정확히 가른다.

둘째 묶음은 수학 박스의 확인이다. 경계 아래 한 칸과 경계 위 한 칸이 모두 DBL_TRUE_MIN, 곧 $2^{-1074}$ 이다.

셋째 묶음은 구멍이 메워졌다는 증거다. $𝑥\ne 𝑦$ 이고, $𝑥-𝑦$ 는 0이 아닌 비정규수로 남는다. 앞의 inverse_gap 은 이 기계에서 안전하다.

넷째 묶음이 대가다. 0.1 을 $2^{-𝑘}$ 로 눌러 바닥에 넣었다가 $2^{𝑘}$ 로 다시 들어 올렸다. 정규수로 남은 k = 1018 에서는 0.1 이 그대로 돌아온다. 비정규수 구간으로 들어가면 8단계마다 유효 비트가 8개씩 빠지고, 돌아온 값이 0.1 에서 점점 멀어진다. 5비트만 남으면 0.1015625 가 되고, 끝내 0이 된다. 「상대 오차는 언제나 DBL_EPSILON 의 절반 이하」라는 정규수의 약속은 바닥에서 지켜지지 않는다 — 앞 절의 상대 오차 비교(near_rel)도 이 구간에서는 기대만큼 일하지 못한다.

나머지 두 묶음은 다음 절의 사고 두 개를 미리 보여 준다.

C 표준이 비정규수를 어떻게 다루는지도 적어 둔다. 표준은 비정규수를 0으로 다룰지 말지를 구현 정의로 남긴다 (§5.2.5.3.3). 그러나 IEC 60559(곧 IEEE 754) 산술을 따르는 구현은 비정규수를 언제나 0이 아닌 값으로 다룬다. 비정규수가 있는지를 알리던 DBL_HAS_SUBNORM 같은 매크로는 C23 에서 낡은 기능(obsolescent)이 되었다. 오늘은 값을 분류하는 fpclassify 의 FP_SUBNORMAL(C99)과 issubnormal(C23, §7.12.3.9)을 쓴다.

54.4 바닥에서 터진 보안 사고들#

비정규수는 아주 작은 수의 이야기라 실무와 멀어 보인다. 그런데 이 바닥은 세 번, 서로 다른 방식으로 보안 사고의 무대가 되었다. 세 사고는 각각 비정규수의 다른 성질을 찌른다.

54.4.1 사고 1 — 숫자 하나로 서버를 멈추다#

2011년 1월 초, 부동소수점 변환을 연구하던 릭 리건(Rick Regan)이 이상한 현상을 블로그에 적었다(PHP 의 버그 추적기에는 그 며칠 전에 보고되어 있었다). PHP 에 2.2250738585072011e-308 이라는 문자열을 숫자로 바꾸라고 시키면 영원히 돌아오지 않았다. 웹 프로그램은 사용자가 보낸 값을 숫자로 바꾸는 일을 끊임없이 한다. 그러니 공격자는 이 문자열 하나를 요청에 실어 보내기만 하면 서버의 일꾼 프로세스를 하나씩 무한 반복에 빠뜨릴 수 있었다. 이런 공격을 서비스 거부(denial of service, DoS)라 부른다 — 무엇을 훔치는 것이 아니라, 서비스가 다른 사람에게 응답하지 못하게 만드는 것이다. 이 결함은 CVE-2010-4645가 되었고, PHP 5.3.5 와 5.2.17 에서 고쳐졌다.

변환이 끝나지 않는 원리를 차근차근 따라가 보자. 출발점은 이것이다 — 십진 문자열을 double 로 바꾸는 일은 생각보다 어렵다. 정답은 「그 십진수에 가장 가까운 double」 하나로 정해져 있지만, 그것을 찾으려면 정밀하게 계산해야 한다. PHP 의 변환 함수(zend_strtod)는 데이비드 게이 (David Gay)가 1990년에 공개한 dtoa.c 를 바탕으로 했고, 그 알고리즘은 이렇게 일한다.

1. 먼저 어림값을 구한다. 대개 정답이거나 정답의 바로 옆 이웃이다.
2. 아주 큰 정수 계산으로 입력한 십진수와 어림값의 차이를 정확히 잰다.
3. 차이를 보고 어림값을 한 눈금쯤 고친다.
4. 정답이 확인될 때까지 2–3을 되풀이한다.

넷째 단계가 반복이라는 데 주목하자. 이 반복은 「고치면 언젠가 정답에 닿는다」는 믿음 위에서만 끝난다.

이제 문제의 수가 어디에 있는지 보자. 「0 근처의 바닥」 절의 비정규수 시연 다섯째 묶음이 보여 주듯, 이 문자열의 정답은 가장 큰 비정규수 000fffffffffffff 다. 이 십진수는 DBL_MIN(…2014e-308)과 가장 큰 비정규수(…2009e-308) 사이, 거의 한가운데에 있다. 그리고 가장 큰 비정규수 쪽에 아주 조금 더 가깝다. 어림값은 DBL_MIN, 곧 $2^{-1022}$ 로 나왔다. 정답보다 한 눈금 위다. 그러니 셋째 단계는 「한 눈금($2^{-1074}$) 내려라」가 되어야 한다.

고치는 코드의 핵심은 두 줄이었다. 뜻만 살려 옮기면 이렇다.

double correct_once(double rv, double aadj1, double step)
{
    double adj = aadj1 * step;   /* 고칠 양: 한 눈금(step)의 몇 배인가 */
    rv += adj;                   /* 어림값을 고친다 */
    return rv;
}

여기서 adj 는 정확히 한 눈금이 아니라 「반 눈금보다 아주 조금 많이 내려라」에 해당하는 음수였다($-1.001\dots \times 2^{-1075}$). 정상이라면 이 덧셈의 결과는 반올림을 거쳐 한 눈금 아래, 곧 가장 큰 비정규수에 떨어지고 반복이 끝난다. 그런데 32비트 x86 에서는 일이 이렇게 흘렀다.

1. x87 은 계산을 80비트 레지스터에서 한다(88장). 이 레지스터 형식은 지수의 범위가 훨씬 넓어서 $2^{-1022}$ 근처는 바닥이 아니다. 그곳의 수도 그냥 정규수다.
2. 그래서 adj 는 double 로 반올림되지 않고 레지스터에 그대로 남는다. 그 상태로 $2^{-1022}+\text{adj}$ 를 계산해 53비트로 반올림하면, 결과는 $2^{-1022}-2^{-1075}$ 이다. 가장 큰 비정규수와 DBL_MIN 의 정확히 한가운데다.
3. 이 값을 메모리의 double 변수 rv 에 저장하는 순간 두 번째 반올림이 일어난다. 메모리의 double 에서 이 값은 비정규수다. 비정규수는 유효 비트가 52개뿐이니 마지막 한 비트를 더 버려야 한다. 그런데 값이 정확히 한가운데라 동률이고, IEEE 754 의 기본 규칙(동률이면 짝수 쪽)은 위쪽, 곧 DBL_MIN 을 고른다.
4. rv 는 출발점 $2^{-1022}$ 로 돌아왔다. 다음 반복은 똑같은 차이를 재고, 똑같은 adj 를 만들고, 똑같이 제자리로 돌아온다. 끝나지 않는다.

이것을 두 번 반올림(double rounding)이라 부른다. 첫 반올림이 값을 정확히 가운데에 올려놓았고, 둘째 반올림이 동률 규칙으로 그것을 위로 되돌렸다. 정규수였다면 둘째 반올림에서 버릴 비트가 없으니 이런 일이 생기지 않는다. 비정규수가 정밀도를 한 비트 잃는 바로 그 경계에서만 생기는 사고다.

수정은 한 낱말이었다. adj 를 volatile double 로 선언한 것이다. volatile 은 「이 변수는 매번 메모리에서 읽고 메모리에 써라」는 약속이다(14장). 그러면 adj 가 덧셈 전에 메모리를 한 번 거치면서 double 로 반올림되어 정확히 $-2^{-1074}$ 가 된다. 이제 덧셈 결과가 처음부터 가장 큰 비정규수라 두 번째 반올림이 아무것도 바꾸지 않고, 반복이 끝난다. 다만 이것은 이 자리 하나를 막은 처방이지, volatile 이 두 번 반올림 일반을 없애 주는 것은 아니다.

이 사고에서 C 프로그래머가 가져갈 것은 넷이다.

• 숫자를 읽는 코드도 입력의 경계다. 문자열을 받는 모든 곳이 공격면이듯, 그 문자열을 숫자로 바꾸는 함수도 공격면이다. 변환 함수를 직접 짜지 말고 검증된 strtod 를 쓴다.
• 「수렴하면 끝난다」에 기대는 반복에는 상한을 둔다. 반복 횟수의 상한은 수학이 틀렸을 때 피해를 무한에서 유한으로 바꿔 준다.
• 경계값을 시험에 넣는다. DBL_MIN, 그 바로 아래 이웃, DBL_TRUE_MIN, 그리고 그 사이에 떨어지는 십진 문자열. 앞의 비정규수 시연이 바로 그런 시험의 모양이다.
• 계산의 폭이 하나라고 가정하지 않는다. 중간 계산을 더 넓은 형식으로 하는 구현이 있다 (FLT_EVAL_METHOD, 88장). 넓은 중간값이 메모리에 저장되며 한 번 더 반올림되는 자리는 이런 사고의 단골 무대다.

54.4.2 사고 2 — 느린 계산이 비밀을 흘린다#

두 번째 사고를 이해하려면 「부채널」이라는 낱말부터 알아야 한다. 부채널(side channel)은 프로그램의 정식 출력이 아닌 곁가지 — 걸린 시간, 전력 소모, 소리 — 로 정보가 새어 나가는 통로다. 금고 다이얼을 돌릴 때 나는 미세한 딸깍 소리로 번호를 알아내는 영화 속 장면이 바로 부채널 공격이다. 그 가운데 걸린 시간으로 비밀을 알아내는 것을 타이밍 부채널(timing side channel)이라 한다.

암호 분야는 1996년 폴 코처(Paul Kocher)의 논문 이후 이 위험을 잘 안다. 그래서 암호 코드에는 「비밀 값에 따라 걸리는 시간이 달라지면 안 된다」(상수 시간, constant time)는 규칙이 있다. 보통은 두 가지를 조심한다. 비밀에 따라 갈라지는 if 와, 비밀에 따라 달라지는 메모리 주소다. 2015년 캘리포니아 대학교 샌디에이고 연구진(앤드리스코 외, IEEE 보안·프라이버시 심포지엄)은 세 번째를 보였다. 같은 명령이 같은 메모리를 다뤄도, 값에 따라 걸리는 시간이 다르다. FTZ·DAZ 상자에서 본 대로 비정규수가 끼면 곱셈 하나가 수십 배 느려지기 때문이다. 연구진은 이것으로 실제 시스템 두 개를 뚫었다.

첫째 표적은 웹 브라우저였다(Firefox 23–27). 브라우저는 동일 출처 정책(same-origin policy)으로 사이트를 서로 격리한다. 공격자의 페이지가 은행 사이트를 페이지 안의 작은 창 (iframe)으로 띄울 수는 있어도, 그 창에 그려진 내용을 스크립트로 읽을 수는 없다. 이 벽이 무너지면 남의 화면에 뜬 이름, 계좌 정보, 로그인 상태가 공격자에게 넘어간다. 공격은 이렇게 진행되었다.

1. 공격자 페이지가 피해자 사이트를 iframe 으로 띄우고, 가로세로 1픽셀짜리 창으로 픽셀 하나만 보이게 한다. 그 픽셀을 색 필터로 검정 아니면 흰색으로 바꾼 뒤, 200×200 크기로 크게 늘린 상자를 만든다. 공격자는 이 상자의 색을 볼 수 없다. 다만 그 위에 필터를 걸 수는 있다.
2. 그 상자에 SVG 필터 하나를 건다. 이웃 픽셀 값들에 가중치를 곱해 더하는 합성곱 필터인데, 가중치를 전부 1e-42 로 준다. 이 값은 float 의 비정규수다(FLT_MIN 이 약 $1.18\times {10}^{-38}$ 이니까).
3. 픽셀이 흰색이면 0이 아닌 값 × 1e-42 = 비정규수가 되어, 상자 전체의 곱셈과 덧셈이 비정규수 연산이 된다 — 느리다. 픽셀이 검은색이면 0 × 1e-42 = 0 이라 정상 속도로 끝난다 — 빠르다.
4. 스크립트는 화면을 다시 그리는 데 걸린 시간을 잰다(requestAnimationFrame 과 고해상도 타이머). 측정에서 검정과 흰색의 차이는 약 40밀리초였다. 미리 검정과 흰색을 번갈아 그려 기준선을 잡아 두면, 걸린 시간만 보고 픽셀의 색을 맞힐 수 있다.
5. 이것을 픽셀마다 되풀이하면 화면을 한 점씩 읽어 낸다. 링크가 방문한 곳이면 검정, 아니면 흰색이 되도록 꾸며 사용자의 방문 기록을 초당 10개 넘게 알아낼 수 있었고, 남의 사이트에 뜬 글자도 읽어 냈다.

이 공격에는 if 문이 하나도 없다는 점이 무섭다. 필터 코드는 픽셀 색에 따라 갈라지지 않는다. 같은 명령을 같은 순서로 실행한다. 다른 것은 오직 값, 그리고 값에 따라 달라지는 곱셈 한 번의 시간이다. Firefox 는 28판에서 SVG 필터 계산을 그래픽 처리 장치로 옮기면서 이 길을 막았다. 그 사이 Firefox 24 를 바탕으로 한 토르 브라우저(Tor Browser)도 영향을 받았다.

둘째 표적은 C 로 짠 한 줄의 덧셈이었다. 차등 프라이버시(differential privacy) 데이터베이스는 「통계는 알려 주되, 어느 한 사람이 그 안에 있는지는 알 수 없게」 하는 시스템이다. 결과에 계산된 잡음을 섞어 개인을 숨긴다. 연구 대상인 Fuzz 는 바로 타이밍 공격을 막으려고 설계된 시스템이었다. 질의는 행마다 정해진 시간 안에 끝나도록 강제되고, 전용 기계에서 돌았다. 그런데 행별 결과를 모두 더하는 마지막 함수는 C 로 이렇게 짜여 있었다(뜻만 줄였다).

double bag_sum(const char *rows[], int n)
{
    double d = 0;
    for (int i = 0; i < n; i++)
        d += atof(rows[i]);      /* 각 행의 질의 결과를 더한다 */
    return d;
}

공격자는 질의를 이렇게 짠다. 「행마다 0을 내되, 찾는 사람의 행이면 1e-310 을 내라.」 1e-310 은 double 의 비정규수다. 그 사람이 데이터베이스에 없으면 합계 d 는 끝까지 0이다. 0 + 0 은 빠르다. 그 사람이 있으면 그 행에서 d 가 비정규수가 되고, 그 뒤의 모든 덧셈이 비정규수 + 0 이 된다. 「0 근처의 바닥」 절 비정규수 시연의 마지막 묶음이 정확히 이 모양이다. 4번 행 하나가 비정규수를 냈을 뿐인데, 열 번의 덧셈 가운데 일곱 번이 비정규수를 만졌다. 행이 100만 개면 찾는 사람의 행 하나가 그 뒤의 덧셈 수십만 번을 느리게 만든다. 누적 변수 하나가 증폭기가 된 것이다.

실측에서 100만 행 질의는 그 사람이 없을 때 50.309초, 있을 때 50.489초가 걸렸다. Fuzz 가 온갖 노력으로 줄여 놓은 시간 흔들림은 0.009초였으니, 0.18초의 차이는 네트워크 너머에서도 또렷하게 보인다. 「어느 한 사람이 있는지 알 수 없게」라는 이 시스템의 존재 이유가 덧셈 한 줄에서 무너졌다.

연구진이 내놓은 대책은 7장의 고정소수점으로 돌아가는 것이었다. 정수 명령만 써서 언제나 같은 시간에 끝나는 소수 연산 라이브러리(libfixedtimefixedpoint)를 만들었고, 덧셈 15클록, 곱셈 43클록이 걸렸다. FTZ·DAZ 를 켜면 시험한 프로세서에서는 시간 차이가 사라졌지만, 연구진은 그것이 앞으로도 보장된다는 근거는 없다고 적었다. 2017년의 후속 연구(콜브레너·샤참, USENIX 보안 심포지엄)는 주요 브라우저들이 그 사이 넣은 방어가 여전히 모자라고, 비정규수가 아닌 경우에도 부동소수점 연산 시간이 값에 따라 달라진다는 것을 보였다.

C 프로그래머가 가져갈 규칙은 분명하다. 비밀 값 — 암호 키, 개인 기록, 인증 정보 — 은 부동소수점 연산에 넣지 않는다. 암호 라이브러리가 정수만 쓰는 데는 이유가 있다. 그리고 「상수 시간」을 따질 때는 분기와 메모리 주소뿐 아니라 값에 따라 속도가 달라지는 명령까지 보아야 한다. 정수 나눗셈도 프로세서에 따라 그런 명령이다.

54.4.3 사고 3 — 남의 라이브러리가 내 계산을 바꾼다#

세 번째 사고에는 공격자가 없다. 그래도 보안의 한 기둥인 무결성 — 계산 결과를 믿을 수 있는가 — 이 흔들린 사건이다.

GCC 와 Clang 의 -ffast-math(그리고 그것을 포함하는 -Ofast)는 「IEEE 754 의 엄밀함을 조금 포기하고 빠르게」라는 옵션이다. 이 옵션으로 링크하면 컴파일러는 작은 시작 코드 (crtfastmath.o)를 끼워 넣는다. 그 코드는 프로그램이 시작할 때 — 또는 그 코드가 든 공유 라이브러리가 적재되는 순간 — FTZ 와 DAZ 를 켠다. FTZ·DAZ 상자에서 본 대로 이 스위치는 스레드 전체의 상태다.

2022년, 보안 연구자 브렌던 돌런개빗(Brendan Dolan-Gavitt)은 파이썬으로 기계학습 모델을 불러오다가 numpy 가 내는 이상한 경고를 보았다. 「가장 작은 비정규수의 값이 0이다.」 그가 그 프로세스에 적재된 공유 라이브러리 158개를 하나씩 불러 보며 범인을 찾자, 네트워크 라이브러리 gevent 가 나왔다. 수치 계산과 아무 상관 없는 라이브러리가 -Ofast 로 빌드되어, 적재되는 순간 프로세스 전체의 비정규수를 꺼 버린 것이다. 그가 파이썬 패키지 저장소(PyPI)를 훑어보니 -ffast-math 로 빌드된 공유 라이브러리가 5,830개, 49개 패키지에서 나왔고, 그것에 기대는 패키지가 2,500개를 넘었다.

이 상태에서는 비정규수가 지켜 주던 것이 모두 거꾸로 돈다. x != y 인데 x - y == 0 이 되고, inverse_gap 은 0으로 나누고, DBL_TRUE_MIN 은 0처럼 군다. 그것도 내가 짠 코드는 한 줄도 바뀌지 않았는데 말이다. 이 문제는 GCC 버그 55522 로 2012년부터 보고되어 있었고, GCC 13 부터 공유 라이브러리(-shared)를 만들 때는 이 시작 코드를 넣지 않게 되었다.

가져갈 것은 둘이다. 라이브러리를 -ffast-math·-Ofast 로 빌드해 배포하지 않는다. 그리고 부동소수점 환경은 내 것이 아니라 프로세스 전체의 것이라고 여긴다. 정밀도가 중요한 프로그램이라면 시작할 때 DBL_TRUE_MIN 같은 값으로 바닥이 살아 있는지 스스로 점검해 볼 수 있다. numpy 가 경고를 낸 것도 바로 그런 자기 점검 덕분이었다.

54.5 IEEE 754 만 있는 것은 아니다#

이 장은 IEEE 754 를 중심으로 이야기했다. 끝으로 IEEE 754 말고 실수를 담는 다른 방법들을 상자 하나로 훑어 둔다.

근사의 세계를 C에서 다루는 법을 익혔다. 다음 장은 이 부의 중심 주제 — 계산이 실패할 수 있다는 사실을 프로그램이 어떻게 다루는가, 오류와 계약의 이야기다.