Proven C BookEnglish GitHub

73 수 — <math.h>·<fenv.h>·<tgmath.h>

먼저 알아야 할 것

50장 실수 · 근사의 수학
8장 수의 표현 · IEEE 754

돌아보기

50장에서 실수를 ==로 비교하지 말라고 했고, 8장에서는 0.1이 정확히 표현되지 않는다고 했다. 그러면 수학 함수가 “계산할 수 없는 입력”을 받으면 어떻게 알려 주는가?

답. 두 경로가 있다. 반환값으로 NaN이나 무한대를 주고, 동시에 errno에 이유를 남긴다 — 정의역 밖이면 EDOM, 결과가 표현 범위를 넘으면 ERANGE. 다만 구현이 errno 대신 부동소수점 예외(floating-point exception) 플래그(<fenv.h>)로 알리도록 선택할 수도 있어서, 이식성 있게 확인하려면 둘 다 볼 준비를 해야 한다. 실무에서는 대개 반환값을 isnan·isinf로 검사하는 쪽이 간단하다.

이 장의 필요성과 맥락

8장(IEEE 754)과 50장(근사의 수학)에서 두 번 다진 바닥 위에 이제 함수들이 올라선다. 세 번째이자 마지막 층인 셈이다. 여기서야 다루는 이유는 하나 — 반올림 모드 같은 숨은 전역 상태는 「계약」이라는 낱말(51장) 없이는 설명할 수 없기 때문이다.

이 장이 끝나면

실수 계산이 실패를 알리는 방식을 본다. 8장과 50장에서 배운 근사의 수학이 여기서 함수의 계약이 된다 — 정의역 밖 호출, 범위를 넘는 결과, NaN과 무한대의 성질, 그리고 반올림 모드라는 숨은 전역 상태까지.

이 장에서 답할 질문

  1. math.h의 함수는 실패를 어떻게 알리는가 — 반환값만 봐서는 알 수 없지 않은가?

73.1 NaN과 무한대의 성질

examples/ch73/math.c

#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <string.h>

int main(void)
{
    /* ① NaN 은 자기 자신과도 같지 않다 */
    double nan_v = nan("");
    printf("NaN == NaN : %s\n", nan_v == nan_v ? "true" : "false");
    printf("isnan(NaN) : %s\n", isnan(nan_v) ? "true" : "false");

    /* ② 무한대와 0 나누기(실수는 UB 가 아니다) */
    double inf_v = 1.0 / 0.0;
    printf("1.0/0.0    : %f (isinf: %d)\n", inf_v, isinf(inf_v));
    printf("0.0/0.0    : %s\n", isnan(0.0 / 0.0) ? "NaN" : "number");

    /* ③ 정의역 밖 호출은 errno 와 NaN 으로 알려 온다 */
    errno = 0;
    double r = sqrt(-1.0);
    printf("sqrt(-1)   : %s, errno=%s\n",
           isnan(r) ? "NaN" : "number", errno == EDOM ? "EDOM" : "0");

    /* ④ 범위를 넘으면 ERANGE 와 무한대 */
    errno = 0;
    double big = exp(1000.0);
    printf("exp(1000)  : %s, errno=%s\n",
           isinf(big) ? "inf" : "number", errno == ERANGE ? "ERANGE" : "0");

    /* ⑤ 분류 함수로 값의 종류를 가른다 */
    double vals[] = { 1.0, 0.0, -0.0, inf_v, nan_v, 1e-320 };
    const char *names[] = { "1.0", "0.0", "-0.0", "inf", "NaN", "1e-320" };
    for (size_t i = 0; i < sizeof vals / sizeof vals[0]; i++) {
        const char *kind = "?";
        switch (fpclassify(vals[i])) {
            case FP_NORMAL:    kind = "normal";      break;
            case FP_SUBNORMAL: kind = "subnormal";    break;
            case FP_ZERO:      kind = "0";         break;
            case FP_INFINITE:  kind = "infinite";      break;
            case FP_NAN:       kind = "NaN";       break;
        }
        printf("  %-7s -> %-6s (sign bit %d)\n", names[i], kind, signbit(vals[i]) != 0);
    }

    /* ⑥ 0.0 과 -0.0 은 같다고 비교되지만 부호는 다르다 */
    printf("0.0 == -0.0 : %s\n", 0.0 == -0.0 ? "true" : "false");
    return 0;
}

실행 결과

NaN == NaN : false
isnan(NaN) : true
1.0/0.0    : inf (isinf: 1)
0.0/0.0    : NaN
sqrt(-1)   : NaN, errno=EDOM
exp(1000)  : inf, errno=ERANGE
  1.0     -> normal (sign bit 0)
  0.0     -> 0      (sign bit 0)
  -0.0    -> 0      (sign bit 1)
  inf     -> infinite (sign bit 0)
  NaN     -> NaN    (sign bit 0)
  1e-320  -> subnormal (sign bit 0)
0.0 == -0.0 : true

출력에서 네 가지를 짚는다.

① NaN은 자기 자신과 같지 않다. IEEE 754가 그렇게 정했다. 그래서 x != x가 참이면 x는 NaN이라는 오래된 관용구가 성립하고, 표준 함수로는 isnan(x)를 쓴다. 이 성질 때문에 NaN이 든 배열을 qsort로 정렬하면 비교자의 전순서가 깨져 결과가 무너진다(66장).

② 실수의 0 나누기는 계약 밖이 아니다. 정수 나눗셈과 달리(28장) IEEE 754 환경에서는 무한대나 NaN을 낳는다. 다만 0으로 나눈다는 사실 자체가 대개 버그라는 점은 같다.

sqrt(-1)EDOM, exp(1000)ERANGE. 전자는 정의역 밖, 후자는 결과가 표현 범위를 넘은 경우다. errno를 볼 생각이라면 호출 직전에 0으로 놓아야 한다(75장).

④ 0.0과 −0.0은 ==로 같다. 그러나 부호 비트는 다르고, 1/0.01/-0.0은 각각 +∞와 −∞다. 부호를 구별해야 하면 signbit을 쓴다.

흔한 오해. “실수 비교는 엡실론을 쓰면 안전하다”

50장에서 배운 엡실론 비교는 만능이 아니다. 절대 오차(fabs(a-b) < eps)는 값이 클 때 무의미해지고 — 1e9 근처에서 1e-9는 표현조차 안 된다 — 상대 오차는 0 근처에서 무너진다. 실무의 처방은 상황에 맞는 허용 오차를 정하는 것이지 보편 상수를 쓰는 것이 아니다. 그리고 애초에 비교가 필요 없도록 정수나 고정소수점으로 다룰 수 있는지 먼저 묻는 것이 낫다(8장의 금액 계산 이야기).

문. math.h의 함수는 실패를 어떻게 알리는가 — 반환값만 봐서는 알 수 없지 않은가?

답. 세 갈래로 알린다. 정의역 밖이면(예: sqrt(-1)) NaN을 돌려주고 errnoEDOM을 놓는다. 범위를 넘으면(예: exp(1000)) 무한대를 돌려주고 ERANGE를 놓는다. 그리고 부동소수점 예외 플래그(<fenv.h>)가 함께 선다.

문제는 이 셋이 구현에 따라 어디까지 지켜지는지 다르다는 점이다. 그래서 실무의 관용구는 errno = 0으로 지운 뒤 호출하고, 곧바로 확인하는 것이다 (75장). 값 자체를 검사하려면 isnan·isinf를 쓴다 — x != x 같은 재주보다 뜻이 분명하다.

73.2 자주 틀리는 함수들

함수하는 일함정
pow(x, y)거듭제곱정수 거듭제곱에 쓰면 느리고 부정확할 수 있다
round·nearbyint반올림round는 0에서 먼 쪽, nearbyint는 현재 모드
floor·ceil·trunc정수로 깎기음수에서 방향이 다르다
fmod·remainder나머지부호 규칙이 서로 다르다
abs·fabs절댓값abs는 정수용. 실수에 쓰면 잘린다
atan2(y, x)각도인자 순서가 y, x
isnan·isinf분류매크로다 — 함수 포인터로 쓸 수 없다

표 73.1

pow(x, 2)는 널리 쓰이지만, 정수 제곱이라면 x * x가 빠르고 정확하다. 컴파일러가 최적화해 주는 경우도 많지만 언제나 그런 것은 아니다.

abs를 실수에 쓰는 실수는 특히 조용하다. <stdlib.h>absint를 받으므로 abs(-1.5)는 −1.5를 1로 만들어 버린다. 요즘 컴파일러는 경고하지만 <math.h>를 포함하지 않은 파일에서는 놓치기 쉽다.

73.3 반올림 모드와 부동소수점 예외 — <fenv.h>

부동소수점 연산에는 숨은 전역 상태가 둘 있다.

반올림 모드 — 기본은 “가장 가까운 값, 동점이면 짝수”다. fesetround로 바꿀 수 있고, 바꾸면 그 뒤의 모든 실수 연산이 영향을 받는다.

예외 플래그 — 나눗셈에서 0, 넘침, 부정확 등이 일어나면 플래그가 켜진다. fetestexcept로 읽고 feclearexcept로 지운다. errno보다 정밀하지만, 이 기능을 쓰려면 #pragma STDC FENV_ACCESS ON을 켜야 하고 — 그러면 컴파일러가 실수 연산 재배열을 자제한다 — 최적화가 줄어든다.

반례. -ffast-math를 켜고 NaN을 검사하기

cc -O2 -ffast-math app.c        # 컴파일러에게 "NaN·무한대는 없다고 쳐라"
if (isnan(x)) { /* 이 분기가 통째로 사라질 수 있다 */ }

-ffast-math 계열 선택지는 결합 법칙을 가정하고 NaN·−0.0의 존재를 무시해도 좋다고 컴파일러에 알린다. 속도는 얻지만 검사 코드가 최적화로 사라질 수 있다 — 17장에서 본 “릴리스에서만 나는 버그”의 실수 판이다. 수치의 정확성이 중요한 프로그램에서는 켜지 않는 것이 기본이다.

73.4 타입 제네릭 — <tgmath.h>

sqrtdouble용, sqrtffloat용, sqrtllong double용이다. <tgmath.h>를 포함하면 sqrt(x) 하나로 인자의 타입에 맞는 판이 골라진다 — 58장에서 본 _Generic이 표준 라이브러리에 쓰인 대표 사례다.

편하지만 대가가 있다. 매크로이므로 함수 포인터로 넘길 수 없고, 인자를 두 번 평가하는 구현이 있을 수 있어 부수효과가 있는 표현식을 넣으면 위험하다.

실제 사례. 같은 계산, 다른 답 — 초과 정밀도의 역사

x86의 옛 부동소수점 장치(x87)는 내부적으로 80비트로 계산했다. 그래서 같은 double 연산이 레지스터에 남아 있을 때와 메모리에 저장된 뒤가 달라지는 일이 생겼다 — 최적화 수준을 바꾸면 결과가 미세하게 달라지고, x == y가 참이었다가 거짓이 되기도 했다.

C99가 FLT_EVAL_METHOD로 이 사정을 명시하게 만들었고, 오늘날 64비트 x86은 SSE를 쓰므로 문제가 크게 줄었다. 그러나 “같은 코드, 다른 답”의 가능성은 여전히 컴파일 선택지와 대상 기계에 남아 있다 — 50장에서 “실수 계산은 재현 가능성을 따로 챙겨야 한다”고 한 이유다.

복습 정리

수 요약.

상황쓸 것조심할 것
NaN 검사isnanx == NaN은 언제나 거짓
무한대 검사isinf실수의 0 나누기는 UB가 아니다
값의 종류fpclassify비정규 수의 존재
정의역·범위 오류반환값 + errno호출 직전 errno = 0
정수 제곱x * xpow(x, 2)
실수 절댓값fabsabs(정수용)
타입별 함수<tgmath.h>매크로 — 부수효과 인자 금지
빠른 수학 옵션끄는 것이 기본검사 코드가 사라진다

표 73.2

수를 지났다. 다음 장은 시간 — 표준이 정해 주지 않은 것이 유난히 많은 자리다.