부정확한 부동소수점수 계산에 대하여

 

C 언어로 배정밀도 부동소수점수(floating point number of double type)로 높은 차수의 다항식 값을 계산을 하는 경우의 오차에 대하여 생각해보기로 하자. 우선 0.7 을 (배정밀도이든 단정밀도이든 상관없이) 부동소수점수로 정확하게 기억시키는 방법은 없다. 0.7 을 배정밀도 부동소수점수(double 타입의 수)로 처리하여 십진법수로 표현하면 유효숫자가 겨우 14개 내지 15개 정도 밖에 되지 않는다.  수 1 에 비하여 오차 1.0e-15 는 상대적으로 매우 작은 오차이다. 그러나 차수가 매우 높은 다항식 함수는 아주 짧은 구간(간격이 약 1.0e-15  인 구간) 안에서도 매우 많은 진동을 하여 함수값의 차이는 매우 커질 수가 있으니 주의해야 한다.

식 f(x) = 2*x*x - 1 로 정의되는 함수 f 를 생각해 보자. x 가 폐구간 [-1,1] 안에 있으면 f(x) 도 이 구간 안에 있다. 즉 f 는 f : [-1, 1] -> [-1, 1] 인 함수이다. 더군다나 이 함수는 두 개의 부동점(fixed point, 즉 f(x) = x 인 x)을 갖는 함수이다. 방정식 2*x*x - 1 = x 플 풀면 근이 -0.5 와 1 이다, 실제로 f(1) = 1 이고 f(-0.5) = -0.5 이다. 이제 53 개의 f 를 반복해서 합성한 함수를 g 로 놓고  g(0.7) 을 계산해 보자. g(0.7) 의 값의 오차가 어떠할런지... 

아래는 C 언어로 g(0.7) = f^(53) (0.7) 을 계산해 본 소스이다. 결과는 참값에 비하여 오차가 너무 큰 근사값(즉 유효하지 않은 근사값))이 된다는 것이다. 참고로 g(x) 는 2^53 차(즉, 9007199254740992 차)의 다항식 함수이다.

 

//  Filename: AccumulativeError07.c
//
//          This illustrates that a tiny error gives a big error
//          in the evaluation of a polynomial of very high degree,
//
//   Compile: gcc AccumulativeError07 AccumulativeError07.c
//   Execute: ./AccumulativeError07
//
//     Or
//
//   Compile: cl AccumulativeError07.c
//   Execute: AccumulativeError07
//
//  Date: 2013. 1. 22.
//  Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)

#include <stdio.h>

#define REPEAT 53

double f(double x)
{
    return 2.0*x*x - 1.0;
}

int main()
{
    double x = 0.7;
    double y, t;
    int i;

    printf("x = %g\n", x);
    printf("x = %f\n", x);
    printf("x = %.15f\n", x);
    printf("x = %.16f\n", x);
    printf("x = %.20f\n", x);
    printf("\n");

    printf("f is a function defined by f(x) = 2*x*x - 1\n");
    printf("f^(n) (x) means the n times composite of the function f(x),\n");
    printf("say, f^(2) (x) = f(f(x)), f^(2) (x) = f(f(f(x))), and etc\n");
    printf("\n");

    t = x;
    for (i = 1; i <= REPEAT; i++) {
        y = f(t);
        t = y;
    }

    printf("y = f^(53) (x) = f^(53) (%g) = %g\n", x, y);
    printf("y = f^(53) (x) = f^(53) (%g) = %f\n", x, y);
    printf("y = f^(53) (x) = f^(53) (%g) = %.15f\n", x, y);
    printf("y = f^(53) (x) = f^(53) (%g) = %.16f\n", x, y);
    printf("y = f^(53) (x) = f^(53) (%g) = %.20f\n", x, y);
    printf("Is it true?\n");
    printf("\n");

    printf("Its error is too big.\nThe iterated value shoul be -0.5671977142553364749...\n");
    printf("\n");

    return 0;
}

/*
Output:
x = 0.7
x = 0.700000
x = 0.700000000000000
x = 0.7000000000000000
x = 0.69999999999999996000

f is a function defined by f(x) = 2*x*x - 1
f^(n) (x) means the n times composite of the function f(x),
say, f^(2) (x) = f(f(x)), f^(2) (x) = f(f(f(x))), and etc

y = f^(53) (x) = f^(53) (0.7) = 0.538614
y = f^(53) (x) = f^(53) (0.7) = 0.538614
y = f^(53) (x) = f^(53) (0.7) = 0.538614158130321
y = f^(53) (x) = f^(53) (0.7) = 0.5386141581303205
y = f^(53) (x) = f^(53) (0.7) = 0.53861415813032054000
Is it true?

Its error is too big.
The iterated value shoul be -0.5671977142553364749...
*/