부동소수점은 어떻게 더할까

구현으로 알아보는 오차의 이유

10.0에 0.1을 곱해서 1.0이 되는 경우는 거의 없다.

10.0 times 0.1 is hardly ever 1.0.

— 브라이언 커니핸Brian Kernighan (1974)

프로그래밍을 하다 보면 0.1 + 0.20.3이 다른 현상을 마주치게 됩니다. 대표적으로 다음 파이썬Python 코드에서 볼 수 있는데요.

0.1 + 0.2 == 0.3 # False

이런 현상이 나타나는 이유로, 부동소수점floating point이 그런 숫자를 정확하게 표현할 수 없기 때문이라고 보통 설명합니다. 실제로 0.1은 이진법으로 쓰면 끝이 없습니다.

0.1=0.00011(2)=0.0001100110011(2)\begin{align*} 0.1 &= 0.0\overline{0011}_{(2)} \\ &= 0.0 \thinspace 0011 \thinspace 0011 \thinspace 0011 \dots \thinspace {}_{(2)} \end{align*}

그래서 컴퓨터가 가진 유한개의 비트로는 근사값을 표현할 수밖에 없습니다.

하지만 이런 표현의 부정확성이 곧바로 0.1 + 0.20.3을 다르게 만드는 것은 아닙니다. 왜냐면 이 두 숫자를 부정확하게 표현한 결과가 우연히 같을 수도 있기 때문입니다. 32비트 부동소수점으로 계산하면 이 둘은 보통은 같습니다. 다음 C++ 코드에서 볼 수 있는 것처럼요.

float f1 = 0.1;
float f2 = 0.2;
float f3 = 0.3;
std::cout << (f1 + f2 == f3) << "\n"; // == 1 (true)

한편, 앞서 언급한 파이썬 예시처럼, 실제 문제로 거론되는 것은 64비트 부동소수점의 경우입니다. 이는 C++의 경우 또한 두 숫자는 다릅니다.

double f1 = 0.1;
double f2 = 0.2;
double f3 = 0.3;
std::cout << (f1 + f2 == f3) << "\n"; // == 0 (false)

그러나 부동소수점이 항상 이 둘을 다르게 표현하는 것은 아닙니다. 덧셈을 비롯한 부동소수점 연산에는 결과를 반올림 처리하는 과정이 있는데요. 여기에는 라운딩 룰rounding rule, 또는 라운딩 모드rounding mode라는 이름으로 다양한 방법이 있습니다.

예를 들어, 라운딩 룰로 버림truncate을 하면 둘 다 비트 표현 상 똑같은 오차를 가집니다. 그래서 이 둘은 같게 됩니다.

fesetround(FE_DOWNWARD);
double f1 = 0.1;
double f2 = 0.2;
double f3 = 0.3;
std::cout << (f1 + f2 == f3) << "\n"; // == 1 (true)

따라서 0.1 + 0.20.3이 64비트 부동소수점으로서 다른 이유는, 유한한 비트로 인한 부정확한 표현과 더불어, 반올림 과정에서 만들어진 서로 다른 오차 때문입니다.

실제로 어떻게 덧셈이 연산되기에, 어떤 라운딩 룰에서는 두 숫자가 같으면서, 또 다른 경우에는 다르게 될까요? 여기서는 0.1과 같은 스트링string에서 부동소수점 비트 표현을 구하고, 이런 비트 표현끼리의 덧셈을 직접 구현하면서 알아보겠습니다.

이전 글에서 부동소수점 수를 어떻게 비교하는지 보았습니다. 여기에 덧셈 연산을 마저 만들텐데요. 그러면 0.1 + 0.2 == 0.3 과정에서 일어나는 일을 전부 다루게 됩니다.

리터럴 읽기

먼저 덧셈에 앞서, 부동소수점 리터럴literal 읽기를 구현할 것입니다. 예를 들면 0.1과 같은 스트링으로부터 부동소수점 비트 표현을 만들어내는 것입니다.

이 작업은 컴파일러compiler가 보통 맡아서 해주는데요. 직접 구현하려면 크게 다음 두 작업을 해야합니다.

  1. 리터럴 파싱: 먼저 리터럴을 읽는 작업, 즉 간단한 파싱parsing이 필요합니다. 이 파싱 결과에서 부호sign와 지수exponent, 소수점fraction을 얻을 것입니다. 이전 글에서 다룬 것처럼, 부동소수점은 이 세 가지 정보가 필요합니다.
  2. 비트 표현 만들기: 이로부터 부동소수점 비트 표현을 만들어내야 합니다. 이 때 근사값이 필요할 수 있기 때문에 반올림도 구현해야 합니다.

그러면 차례로 이것들을 구현해봅시다.

리터럴 파싱

가장 먼저 0.1이나 -2 같은 부동소수점 리터럴을 파싱할 것입니다. 즉 부호나 소수 부분은 생략될 수 있습니다. 간단한 구현을 위해 1e2와 같은 경우는 생략하겠습니다.

리터럴을 스트링으로 받아 부호와 숫자를 결과로 주도록 함수를 만들어봅시다. 다음 파이썬 코드처럼, 읽으려는 문자의 위치를 위해 변수 i를 만듭니다.

def readNumberLiteral(literal: str) -> tuple[bool, float]:
    i = 0

첫 문자가 +-인 부호라면 이를 읽습니다.

    # read sign
    positive = literal[i] != "-"
    if literal[i] in "+-":
        i += 1

그 다음으로 정수 부분을 읽습니다. 리터럴이 끝나거나 소수점을 만날때까지 반복해 정수 부분을 계산합니다.

    # read whole numbers
    num = 0
    while i < len(literal) and literal[i].isdigit():
        num = num*10 + (ord(literal[i]) - ord('0'))
        i += 1

소수점이 없다면 여기서 결과를 만들어줍니다.

    # return number if end
    if i == len(literal) or literal[i] != ".":
        return (positive, num)

소수점이 있다면 더 읽고 결과를 냅니다.

    # read decimal part
    i += 1
    factor = 1
    while i < len(literal) and literal[i].isdigit():
        num = num*10 + (ord(literal[i]) - ord('0'))
        i += 1
        factor *= 10
    divided = num / factor

    return (positive, divided if positive else -divided)

이제 이 함수는 다음처럼 리터럴을 읽을 수 있게 됩니다.

positive, num = readNumberLiteral("12.25") # True, 12.25

positive0이 양수인지 음수인지 구분하기 위해 사용할 것입니다. 부동소수점은 0에도 부호가 있기 때문에 이 정보가 필요합니다.

비트 표현 만들기

본격적으로 이진법 표현을 구해야 하는데요. 여기서 사용할 방법은 이렇습니다. 숫자 nn이 주어지면, 그 숫자보다 크지 않은 2의 제곱수 2m2^m을 얻습니다. 이는 nn의 이진법 표현에서 2m2^m가 가장 큰 자릿수라는 뜻입니다. 그러면 다음 자릿수를 구해가며 이진법 표현을 구할 수 있습니다.

예를 들면 숫자 5.55.5가 주어졌다고 해봅시다. 이보다 크지 않은 제곱수로 22=42^2=4를 얻습니다. 그렇다면 가장 큰 자릿수는 222^2의 자리가 됩니다.

222120215.5=1???(2)\begin{align*} &2^2 &&2^1 &&2^0 &&2^{-1} &&\dots \\ 5.5 ={} &1 &&\text{?} &&\text{?} &&\text{?} &&\dots \medspace {}_{(2)} \end{align*}

다음 자릿수를 11로 했을 때의 이진수를 구해봅니다. 그러면 110(2)110_{(2)}, 즉 66인데요. 주어진 숫자 5.55.5보다 크기 때문에, 자릿수 212^1는 이진법 표현에서 00이 됩니다.

222120215.5=10??(2)\begin{align*} &2^2 &&2^1 &&2^0 &&2^{-1} &&\dots \\ 5.5 ={} &1 &&0 &&\text{?} &&\text{?} &&\dots \medspace {}_{(2)} \end{align*}

다음 자릿수 202^011로 한다면 101(2)101_{(2)}, 즉 55입니다. 5.55.5보다 작으므로, 이진법 표현에서 11이 됩니다. 이런 식으로 반복하면 101.1(2)101.1_{(2)}라는 이진수를 얻게 됩니다.

222120215.5=1011(2)\begin{align*} &2^2 &&2^1 &&2^0 &&2^{-1} &&\dots \\ 5.5 ={} &1 &&0 &&\text{1} &&\text{1} &&\dots \medspace {}_{(2)} \end{align*}

부호와 지수는 이로부터 간단히 얻을 수 있습니다. 결과를 노말라이즈하면 1.011×221.011 \times 2^2가 되는데요. 즉 부동소수점 비트 표현이 가질 지수는, 아까 구한 가장 큰 자릿수 222^2의 지수 22가 됩니다.

이렇게 구한 값으로부터, 부호 SS, 지수 EE, 소수 FF의 비트 표현을 만들어낼 수 있습니다.

S=0E=10000001F=01100\begin{align*} & S = \texttt{0} \\ & E = \texttt{10000001} \\ & F = \texttt{011\thinspace{}0\dots0} \end{align*}

이전 글에서 다룬 것처럼, 지수는 127 만큼 바이어스biased 되었고, 소수는 맨 앞의 비트 1이 생략되었습니다.

위 결과는 C++ 코드로도 볼 수 있습니다.

float f = 5.5;
string b = bitset<32>(*reinterpret_cast<int *>(&f)).to_string();
cout << "S: " << b[0] << ", "
     << "E: " << b.substr(1, 8) << ", "
     << "F: " << b.substr(9) << "\n";
// S: 1, E: 10000001, F: 01100000000000000000000

구현: 리터럴을 읽고 부동소수점 만들기

위에서 만든 함수를 이용해, 리터럴로부터 부동소수점을 만들 것입니다. 이전에 만든 Float 클래스의 정적 메소드로 구현해보겠습니다.

class Float:
    @staticmethod
    def fromLiteral(literal: str) -> "Float":
        positive, num = readNumberLiteral(literal)

        if num == 0:
            return Float(0 if positive else 1, 0, 0)

        return Float.fromNumber(num)

만약 읽은 숫자가 0이면 바로 부동소수점을 만들어줍니다. 아니면 나머지 과정을 곧 구현할 fromNumber() 메소드에 넘깁니다. 이 메소드는 숫자로부터 부동소수점을 만들 것입니다.

먼저 부호를 읽습니다. 음수라면 양수로 만들어둡니다.

    @staticmethod
    def fromNumber(number: float) -> "Float":
        sign = 1 if number < 0 else 0
        if sign == 1:
            number *= -1

앞서 소개했던 대로 이진법 표현을 구합니다. 즉 주어진 숫자를 넘지 않는 가장 큰 2의 제곱수를 구한 뒤, 이진법 표현을 구합니다. 이로부터 지수와 소수 부분을 구합니다.

        exp = 0
        digit = 1.0
        if number > 1:
            while digit < number/2:
                exp += 1
                digit *= 2
        if number < 1:
            while digit > number:
                exp -= 1
                digit /= 2

        acc = 0.0
        bits = []
        for _ in range(1+23+3):
            if acc+digit <= number:
                acc += digit
                bits.append(1)
            else:
                bits.append(0)
            digit /= 2

        binaryStr = "".join(map(lambda n: str(n), bits[1:]))
        binary = int(binaryStr, base=2)

부동소수점이 실제로 가질 소수는 23비트입니다. 그런데 실제로는 27비트를 구했는데요. 이전 글에서 본 것처럼, 부동소수점의 소수 비트는 소수점 앞 비트 1을 생략해서 사실상 24비트의 정보를 담습니다. 여기에 이후 근사값을 구하기 위해 3비트를 더 구해놓았습니다.

나머지 부분에서 반올림을 한 뒤 부동소수점을 만듭니다.

        rounded = roundToNearestEven(binary)

        if rounded >= (1 << 23):
            rounded -= (1 << 23)
            exp += 1

        return Float(sign, exp+127, rounded)

이제 마지막으로 반올림만 남았습니다.

라운딩 룰 구현하기

리터럴을 읽거나 덧셈 등의 사칙 연산을 할 때, 그 결과를 반올림해 근사값이 계산되는데요. 앞에서 언급한 것처럼, 반올림을 처리하는 방법은 여러가지가 있습니다.

그 중에 기본적으로는 가까운 짝수로 반올림round to nearest even하는 방법이 사용됩니다. 한편 앞에서 예시로 본 것처럼 버림을 사용할 수도 있습니다. 여기서는 이 두 가지를 만들어보겠습니다.

가까운 짝수로 반올림

먼저, 보통 때는 일반적인 반올림을 하면서, 정확히 반인 경우는 짝수로 반올림하는 방법입니다.

0.01(2)0(2)0.1(2)0(2)(to even)0.11(2)1(2)1.01(2)1(2)1.1(2)10(2)(to even)1.11(2)10(2)\begin{align*} &0.01_{(2)} &\rightarrow{} &&0_{(2)} &{} \\ &0.1_{(2)} &\rightarrow{} &&0_{(2)} &{} \quad \text{(to even)} \\ &0.11_{(2)} &\rightarrow{} &&1_{(2)} &{} \\ &1.01_{(2)} &\rightarrow{} &&1_{(2)} &{} \\ &1.1_{(2)} &\rightarrow{} &&10_{(2)} &{} \quad \text{(to even)} \\ &1.11_{(2)} &\rightarrow{} &&10_{(2)} &{} \end{align*}

여기서는 마지막 3개의 비트를 사용해 반올림을 할 것인데요. 즉 0101 (0.101(2)0.101_{(2)}) 같은 경우는 1 (1(2)1_{(2)})로 만드는 식입니다. 이 3개의 비트는 큰 자릿수부터 각각 부동소수점 용어로 가드guard, 라운드round, 스티키sticky 비트라고도 불립니다.

그러면 받은 비트 중에 마지막 3개의 비트를 이용해 반올림하도록 만듭시다.

def roundToNearestEven(binary: int) -> int:
    """
    round using the last three bits
    """

    least = (binary & 0b1000) >> 3
    last = binary & 0b111 # guard, round, sticky bits

    if last < 0b100:
        return binary >> 3
    if last > 0b100:
        return (binary >> 3) + 1

    # round to nearest even
    if least == 0:
        return binary >> 3
    else:
        return (binary >> 3) + 1

이 함수는 다음과 같이 쓸 수 있습니다.

roundToNearestEven(0b0101) # == 1 (0.101 rounded to 1)

이제 리터럴로부터 부동소수점을 얻을 수 있습니다.

f = Float.fromLiteral("0.1")
assert f.sign == 0
assert f.exp  == 0b01111011
assert f.frac == 0b10011001100110011001101

이 결과는 C++ 코드로 얻은 결과와도 일치합니다.

float f = 0.1;
string b = bitset<32>(*reinterpret_cast<int *>(&f)).to_string();
cout << "S: " << b[0] << ", "
     << "E: " << b.substr(1, 8) << ", "
     << "F: " << b.substr(9) << "\n";
// S: 1, E: 01111011, F: 10011001100110011001101

버림

버림은 마지막 세 개의 비트를 단순히 버림으로써 간단히 구현할 수 있습니다.

def truncate(binary: int) -> int:
    """
    truncate the last three bits
    """

    return (binary >> 3)

여기서 만든 함수들은 다음 덧셈 구현에서 사용할 것입니다.

부동소수점 덧셈 구현

이제 실제로 덧셈을 할 일만 남았는데요. 부동소수점은 그 덧셈 방법이 알려져 있습니다. 이를 간단히 알아보고 코드로 옮겨봅시다.

방법

덧셈은 두 수 중에 지수 부분이 큰 쪽을 기준으로 자리를 맞춰 이루어집니다. 그리고 결과를 반올림합니다.

0.1 + 0.2를 예로 들어봅시다. 두 리터럴을 읽었을 때, 각각 다음과 같은 부동소수점으로 파싱하게 됩니다.

0.1:
  sign: 0
  exp:  01111011 (-4)
  frac: 10011001100110011001101

0.2:
  sign: 0
  exp:  01111100 (-3)
  frac: 10011001100110011001101

0.2 쪽의 지수가 더 크기 때문에, 이를 기준으로 자리를 맞춰 덧셈합니다. 한편, 덧셈 시에는 소수점 앞에 생략된 비트 1도 고려해야 합니다.

  0.2:  1.10011001100110011001101
  0.1:  0.110011001100110011001101
       ---------------------------
   + : 10.011001100110011001100111
                                ^^^

이 결과에서 25번째부터 세 개의 비트, 여기서는 110이 반올림에 사용됩니다. 만약 가까운 짝수로 반올림하는 라운딩 룰을 적용한다면 1을 더하게 될 것입니다.

이제 덧셈 결과를 노말라이즈하고, 소수점 앞 비트 1을 생략하면 23개의 소수 비트는 다음과 같습니다.

   + : 1.0011001100110011001100111
                                ^^^
       ----------------------------
         00110011001100110011010

이는 실제로 C++에서 계산한 결과와 같습니다.

float f1 = 0.1;
float f2 = 0.2;
float f = f1 + f2;
string b = bitset<32>(*reinterpret_cast<int *>(&f)).to_string();
cout << "F: " << b.substr(9) << "\n";
// F: 00110011001100110011010

구현하기

이제 Float 클래스에 덧셈을 구현할 것입니다. 구현을 간단히 하기 위해 부호가 다른 경우는 제외하겠습니다.

    def __add__(self, other) -> "Float":
        if not isinstance(other, Float):
            raise Exception("not comparable")
        if self.sign != other.sign:
            raise Exception("not implemented for different signs")

먼저 지수가 작은 쪽을 고른 다음, 0과 같으면 큰 수를 새로 리턴합니다. 이 때 이전에 만든 isZero() 메소드를 사용합니다.

        greater, smaller = (self, other) if self.exp >= other.exp else (other, self)

        if smaller.isZero():
            return Float(greater.sign, greater.exp, greater.frac)

한편 덧셈에서는 맨 앞에 생략된 비트 1을 고려해야 합니다. 그리고 나중의 반올림을 위해, 뒤에 세 개의 비트를 붙입니다.

        # prepend omitted bit
        greaterFrac = greater.frac | (1 << 23)
        smallerFrac = smaller.frac | (1 << 23)

        # append three bits
        greaterFrac <<= 3
        smallerFrac <<= 3

작은 쪽의 수를 자리에 맞추고 더합니다. 덧셈 결과를 노말라이즈 하기 위해, 24 비트가 넘어간 만큼 지수에 반영합니다. 마지막으로 맨 앞의 비트 1을 생략하면, 부동소수점 비트 표현을 만들 수 있습니다.

        # align by exp
        smallerFrac >>= greater.exp - smaller.exp

        added = greaterFrac + smallerFrac
        rounded = roundToNearestEven(added)

        exp = greater.exp

        large = 1 << 24
        while rounded >= large:
            exp += 1
            rounded >>= 1

        # get last 23 bits to drop the first bit
        rounded &= ((1 << 23) - 1)

        return Float(sign, exp, rounded)

이제 다음과 같이 덧셈이 가능합니다.

f = Float.fromLiteral("0.1") + Float.fromLiteral("0.2")
assert f.sign == 0
assert f.exp  == 0b01111101
assert f.frac == 0b00110011001100110011010

이전에 만든 비교 연산으로, 0.3과 같음을 알 수 있습니다. 이 결과는 앞에서 소개한 예시와 일치합니다.

assert f == Float.fromLiteral("0.3")

반올림 모드 추가하기

반올림은 한 가지 방법만 있는 것은 아닙니다. 위에서 만든 버림도 사용해볼 수 있는데요. 이를 위해 Float 클래스에 반올림 모드를 다음처럼 추가합시다. 이는 기본적으로 가까운 짝수로 반올림하도록 만듭니다.

class Float:
    roundMode = "NEAREST_EVEN"

그리고 __add__() 메소드에서 반올림 모드에 따라 다르게 반올림하도록 바꿉시다.

    def __add__(self, other) -> "Float":
        # ...

        # add and round
        added = greaterFrac + smallerFrac
        rounded = 0
        if Float.roundMode == "NEAREST_EVEN":
            rounded = roundToNearestEven(added)
        else:
            rounded = truncate(added)

이렇게 버림을 이용해 아까와 같이 0.1 + 0.2를 계산하면, 이번에는 0.3과는 다른 결과가 나타납니다.

Float.roundMode = "TRUNCATE"
f = Float.fromLiteral("0.1") + Float.fromLiteral("0.2")
assert f != Float.fromLiteral("0.3")

이는 C++로 같은 결과를 구했을 때와 마찬가지입니다.

fesetround(FE_DOWNWARD);
float f1 = 0.1;
float f2 = 0.2;
float f3 = 0.3;
float f = f1 + f2;
std::cout << (f1 + f2 != f3) << "\n"; // 1

이렇게 반올림 모드와 함께 덧셈 연산을 완성했습니다. 여기까지 작성한 코드는 지스트Gist에서 확인할 수 있습니다.

그 외의 것들

본문에서 생략했지만 언급할 만한 것들로 이런 것이 있습니다.

덧셈의 결합성

모든 숫자 a,b,ca, b, c에 대해, 다음을 만족하면 덧셈이 결합법칙을 만족한다associative고 합니다.

(a+b)+c=a+(b+c) (a + b) + c = a + (b + c)

즉 어디서부터 더하든지 같은 결과를 내는 성질입니다.

부동소수점은 이런 간단한 성질도 만족시키지 않습니다. 아래 파이썬 코드에서 확인할 수 있듯이요.

(0.1+0.2)+0.3 == 0.1+(0.2+0.3) # False

이는 본문에서 구현한 반올림 때문이라는 것을 알 수 있습니다.

실제 구현

여기서는 파이썬으로 만들어봤지만, 실제로는 대부분 CPU의 명령어 세트instruction set에서 부동소수점 연산을 지원합니다. 즉 하드웨어 수준에서 이미 구현되어 있습니다. 따라서 이렇게 소프트웨어 상에서 구현한 것보다 훨씬 빠르게 동작하게 됩니다.

특별한 숫자를 위한 비트 표현

부동소수점에는 무한과 ‘숫자가 아님’, 즉 NaNnot-a-number을 표현하는 비트 표현도 있는데요. 본문에서는 간단한 구현을 위해 생략했지만, 실제 부동소수점 사용에서 종종 마주칠 수 있는 값입니다. 이것들은 특별한 비트 표현으로 예약되어 있습니다.

비트의 개수와 실수

32비트로 표현할 수 있는 경우의 수는 2322^{32}개 뿐입니다. 하지만 0011 사이의 수만 보더라도 이보다 훨씬 많은데요. 이는 사실 집합론set theory에서 알 수 있듯이 실수real number만큼이나 많습니다. 그래서 얼마나 많은 비트를 사용하든지 무관하게, 비트는 실수를 표현하기에는 항상 압도적으로 부족하게 됩니다.

마치며

여기까지 32비트 부동소수점의 비교 연산과 덧셈 연산을 간단히 구현해보았습니다. 컴파일러와 같은 프로그램이 부동소수점 리터럴로부터 비트 표현을 얻는 과정을 따라해보고, 이렇게 얻은 부동소수점 비트 표현끼리 덧셈하는 함수를 만들어보았습니다. 이 과정에서 원래의 숫자를 제한된 비트의 개수로 표현하기 위해 반올림을 했는데요. 라운딩 룰 중에 가까운 짝수로 반올림하는 방법은 오차를 최소화하려는 방법이지만, 경우에 따라 어쩔 수 없이 나타나는 오차를 보았습니다.

이 내용은 64비트의 경우 또한 이와 비슷하므로 쉽게 확장할 수 있습니다.

레퍼런스

  • Computer Organization and Design (David Patterson, John Hennessy, 2013), 또는 컴퓨터 구조 및 설계 (2015)