June 11, 2020 16:20

서론

지난 글에서는 부동소수점 자료형이 무엇이고, 이들의 대표적인 표준인 IEEE 754에서 규정한 이진법 16비트 자료형의 세부적인 구조 및 10진수와 서로 변환하는 법 등을 알아보았습니다. 이번 글에서는 지난 글에 이어, 부동소수점 자료형으로 나타낸 수들끼리의 사칙연산을 수행하는 과정을 알아보고, 실제 프로그래밍 언어에서의 수학 라이브러리 함수를 보며 비교적 간단한 연산들을 어떻게 구현하게 되는지에 대해 알아보도록 하겠습니다.

사칙연산

사칙연산과 같은 기본적인 연산은 CPU, 또는 FPU (floating-point unit)에 내장되어 있기 때문에 별도의 라이브러리 없이도 사용이 가능합니다. 이 문단에서는 부동소수점 자료형을 사용하여 어떤 방법으로 사칙연산을 수행할 수 있는지 예시와 함께 설명하도록 하겠습니다. 실제 하드웨어상의 계산 로직은 속도 향상을 위해 병렬적으로 수행하거나 더 개선된 알고리즘을 사용하는 것이 보통이지만 이 글에서는 자세히 다루지는 않겠습니다.

덧셈 및 뺄셈

덧셈이나 뺄셈을 수행하기 위해서는 먼저 지수부를 동일하게 맞춰주어야 합니다. 가수부가 나타내는 실제 값은 지수부에 따라 서로 다르기 때문입니다. 따라서 한 쪽의 지수가 다른 쪽과 같아지도록 지수를 다른 쪽으로 변경하고 가수부를 그 차이만큼 시프트한 뒤 가수부끼리를 더하면 됩니다. 뺄셈은 두 번째 수의 부호를 바꾸고 덧셈을 수행하면 되므로 원리는 같다고 할 수 있습니다.

예를 들어, 이진법의 수 $1.00101 \times 2^2$과 $1.11 \times 2^{-1}$을 서로 더하는 과정은 다음과 같습니다. 우선, 둘의 지수부를 동일하게 맞추기 위해 두 번째 수의 지수를 2의 제곱이 되게끔 바꾸어줍니다.

$1.11 \times 2^{-1} = 0.00111 \times 2^2$

이제 가수부끼리는 그냥 더할 수 있습니다.

$1.00101 + 0.00111 = 1.01100$

따라서 결과는 $1.01100 \times 2^2$로 표기할 수 있습니다. 덧셈의 결과로 유효숫자의 시작 자리가 바뀔 수도 있는데, 그런 경우 $1.xxx$의 형태가 되도록 지수를 정규화 (normalize) 해주면 됩니다.

부동소수점 자료형의 특성 때문에 덧셈이나 뺄셈을 수행한 결과에서는 낮은 자리수들이 크게 손실되는 경우가 자주 발생합니다. 두 수의 비가 큰 경우에는 덧셈 결과가 아예 변하지 않을 수도 있는데, 이러한 현상을 cancellation, 또는 loss of significance라고 합니다. 다시 말해, 어떤 수를 계속 더해나가면 수가 점점 커질 것을 기대하지만 실제로는 영원히 어떠한 변화도 생기지 않는 시점이 오게 됩니다.

예를 들어, 이진법의 수 $1.0 \times 2^{30}$과 $1.0 \times 2^0$를 더한 결과는 $1.000000000000000000000000000001 \times 2^{30}$이 되어야 하는데, 가수부에 할당된 비트 수에 따라서는 낮은 자리의 1을 표현하는 것이 불가능할 수도 있기 때문입니다. IEEE 754에서는 이러한 계산 결과가 이전 글에서 언급했던 반올림 규칙에 의해 표현 가능한 가장 가까운 수가 되도록 정하고 있습니다.

곱셈 및 나눗셈

곱셈의 경우도 비교적 간단합니다. 가수부끼리, 지수부끼리 각각 곱한 뒤 정규화를 해주면 됩니다. 다음의 예시를 통해서 보겠습니다.

$1.101 \times 2^2$과 $1.01 \times 2^{-1}$의 곱을 수행하기 위해 먼저 가수부의 $1.101$과 $1.01$을 곱합니다. 그 결과로 $10.00001$을 얻습니다. 그 다음 지수부끼리를 곱해야 하는데, 밑이 같은 지수들이므로 지수끼리를 더한 것과 같습니다. 따라서 $2^2 \times 2^{-1} = 2^{2-1} = 2^1$이 됩니다.

이대로 쓰면 $10.00001 \times 2^1$인데, 이 수는 정규화된 수가 아니므로 소수점을 왼쪽으로 한 칸 옮기고 지수를 1 증가시켜 $1.000001 \times 2^2$으로 표시해야 합니다.

나눗셈도 곱셈과 비슷한 방법으로 수행할 수 있습니다. 가수부끼리, 지수부끼리를 각각 나눗셈을 수행한 뒤 정규화를 해주면 됩니다.

곱셈과 나눗셈에서는 cancellation 문제가 발생하지는 않지만, 역시 반올림에 의한 오차는 발생할 수 있습니다. 곱셈의 경우 두 수의 정확한 계산 결과는 가수부가 두 배가 될 것을 요구하지만, 낮은 자리수 반은 표현을 할 수 없기 때문에 반올림에 의해 손실되게 됩니다. 나눗셈의 경우에는 무한소수가 되어 일정 자릿수 아래를 반올림해야 하는 경우도 발생할 수 있습니다.

수학 라이브러리

상용 프로그래밍 언어들에서는 일반적으로 방대한 양의 수학 라이브러리를 제공합니다. 그 중 하나로 C, C++의 math.h가 있습니다. FPU는 이 라이브러리에 있는 함수들 중 일부를 직접 내장 명령어로 처리할 수 있기도 하지만, 이 글에서는 라이브러리 코드를 보면서 정수형들을 사용하여 비트 단위로 어떻게 연산이 가능한지를 중점으로 알아보도록 하겠습니다. 부호 비트, 지수부, 가수부들을 어떻게 분리하고 값을 조정하여 원하는 연산을 수행하게 되는지를 이해하는 것이 중요한 목표가 될 것입니다.

이 문단에서 참조하는 코드는 glibc입니다.

floor, ceil, round

• floor: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/flt-32/s_floorf.c
• ceil: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/flt-32/s_ceilf.c
• round: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/flt-32/s_roundf.c

매우 간단한 편에 속하는 floor, ceil, round부터 살펴보겠습니다. 이 함수들은 각각 주어진 float 값을 가까운 정수값으로 내림, 올림, 반올림하는 역할을 합니다. 세 함수 모두 전체적인 구조가 비슷한 것을 볼 수 있는데, 어느 방향으로 변환할지를 결정하는 부분을 제외하고는 거의 같은 과정을 거치기 때문입니다.

먼저 세 함수에 모두 공통적으로 있는 GET_FLOAT_WORD와 SET_FLOAT_WORD 매크로 함수가 눈에 띕니다. 두 매크로는 각각 float형의 값을 비트 구조 그대로 정수형의 레지스터에 넣는 것과 그 역연산을 수행합니다. 이 매크로들은 아키텍처별로 선언된 math_private.h에 그 아키텍처에서 가장 효율적으로 동작할 수 있는 코드 덩어리로 정의되어 있습니다. 정수형 레지스터에 담는 이유는 이전 글에서도 언급했듯이 비트 단위의 연산은 실수형에 대해서는 지원되지 않아 직접적으로 비트를 다룰 수 없기 때문입니다. 실제로 반환해야 하는 값은 float형이므로 연산이 끝나면 다시 이 값을 실수형 레지스터에 넣고 반환하게 됩니다.

다음 세 함수에 공통적으로 있는 j0 = ((i0>>23)&0xff)-0x7f;는 이 float 값의 지수부에서 bias를 제거한 값을 구하는 과정입니다. 가수부인 오른쪽의 23비트를 시프트 연산을 통해 제거하고, 0xff와 AND 연산을 해서 부호 비트를 제거하고, 32비트 float의 bias인 127 (0x7f)를 빼서 구한 것입니다.

그 다음은 j0가 23보다 작은지를 기준으로 분기를 나누고 있는데, 이는 쉽게 말해서 이 값이 이미 온전한 정수인지를 구분하기 위함입니다. 가수부가 23자리이므로, 지수부가 23 이상이라면 가수부의 가장 낮은 자리도 정수부에 속하기 때문입니다. 정수인 경우에는 floor, ceil, round에 관계 없이 항상 자기 자신을 반환하면 되기 때문에 별도로 처리를 해줍니다.

j0가 23보다 작다면 소수부에 속하는 가수부 비트가 있다는 뜻이므로 추가 작업이 필요해집니다. 여기서도 다시 분기를 나누는데, j0가 0보다 작다면 모든 가수부가 실수부에 속한다는 의미이므로, 즉, 수의 절댓값이 1보다 작기 때문에 각 함수에 맞게 -1, 0, 1 중 하나의 값을 만들어 반환하게 됩니다.

그렇지 않은 경우에는 가수부가 정수부와 소수부에 걸쳐있기 때문에 소수부가 시작되는 위치를 구해서 추가적인 작업이 수행되어야 합니다. 각 함수는 다음과 같이 작업할 수 있습니다.

• floor: 소수부의 모든 비트를 0으로 변경해주면 됩니다.
• ceil: 소수부의 모든 비트가 0이면 그대로 반환하고, 그렇지 않으면 소수부를 0으로 만든 뒤 정수부를 1 증가시키면 됩니다.
• round: 소수부의 모든 비트가 0이면 그대로 반환하고, 그렇지 않으면 소수 첫째 자리에 1을 더한 값에서 소수부를 0으로 만들면 됩니다.

modf

• 코드: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/flt-32/s_modff.c

modf 함수는 실수를 정수부와 소수부로 분리해주는 역할을 합니다. 이 라이브러리 함수의 원형은 다음과 같습니다.

float       modff( float arg, float* iptr );
double      modf( double arg, double* iptr );
long double modfl( long double arg, long double* iptr );

arg는 분리할 값이며, 정수부는 iptr가 가리키는 위치에 저장되고 반환값으로 실수부가 전달됩니다. 예를 들어 arg에 123.45를 전달했다면 *iptr는 123, 반환값은 0.45가 되는 방식입니다.

분기를 나누는 과정은 위의 floor, ceil, round에서와 동일합니다. 소수부가 있는 경우와 없는 경우로 나누고, 없는 경우 모든 값을 *iptr에 넣고 반환값은 ±0, 있는 경우 지수부가 23 미만이면 모든 값을 반환값으로 하고 *iptr은 ±0이 됩니다.

정수부와 소수부가 모두 존재하는 경우에도 역시 마찬가지로 둘이 나뉘는 정확한 지점을 찾고, 이 값이 순수한 정수인지 아니면 소수점 아래에도 값이 있는지를 판별하여 분기를 나눕니다. 정수값인 경우에는 모든 값이 *iptr에 들어가고 반환값은 ±0이 됩니다. 그렇지 않은 경우, i0&(~i)를 통해 정수부만을 떼어 *iptr에 저장하고 원래 값에서 *iptr을 뺀 값을 반환하게 됩니다.

fmod

• 코드: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/flt-32/e_fmodf.c

fmod 함수는 modf와 이름이 비슷한 만큼 역할도 비슷합니다. fmod는 하나의 float값을 다른 float값으로 나눈 나머지를 반환하는 함수인데, 나누는 값을 1로 설정한 것이 곧 modf와 같다고 볼 수 있습니다. 즉, fmod 함수는 modf 함수의 일반화 버전입니다.

fmod 함수의 원형은 다음과 같습니다.

float       fmodf( float x, float y );
double      fmod( double x, double y );
long double fmodl( long double x, long double y );

쉽게 말해 x를 y로 나눈 나머지가 반환됩니다. 몫을 구하고 싶다면 x에서 반환값을 빼고 다시 나누기를 수행하면 됩니다.

이 함수의 내부 구조는 위에서 살펴본 함수들과는 사뭇 다를 뿐 아니라 훨씬 복잡한 과정을 거치는 것을 볼 수 있습니다. 천천히 부분별로 따라가 봅시다.

변환

먼저 x와 y를 비트 단위로 다루기 위해 변환 과정을 거쳐야 합니다. 이 두 값을 각각 hx와 hy라는 32비트 정수형 변수에 담는데, 추가로 sx라는 변수에 x의 부호값을 담고, hx와 hy에는 절댓값을 취하는 작업을 하고 있습니다. 이렇게 부호를 제거하고 연산 후 마지막에 다시 부호를 추가하는 과정은 라이브러리에서 자주 볼 수 있는데, 보다 일반화된 코드로 효율적으로 연산하기에 용이하기 때문입니다.

예외 처리

나누기 연산에는 중간에 발생할 수 있는 예외 상황이 몇 가지 있기 때문에 그런 경우들을 먼저 처리해주고 시작해야 합니다. 크게 다섯 가지 경우입니다.

• y가 0일 때
• x가 유한하지 않을 때¹
• y가 NaN일 때

이 세 경우에는 마땅한 계산을 수행할 수 없으므로 NaN이 반환되도록 만들어 줍니다. 그 외에도 아래 두 경우를 특수하게 예외로 처리합니다.

• x의 절댓값이 y의 절댓값보다 작을 때: x를 그대로 반환
• x의 절댓값과 y의 절댓값이 같을 때: x부호에 따라 ±0을 반환

이러한 예외 처리 역시 이후 연산들을 보다 일반성 있게 할 수 있도록 만들어주는 역할을 한다고 볼 수 있습니다.

지수부 구하기

다음은 x와 y의 지수부를 구하는 과정인데, 일반적인 정규화된 표기법이라면 위에서 본 식과 똑같이 ix = (hx>>23)-127;을 해서 구할 수 있습니다. 그런데 여기서는 예외로 처리해야 할 사항이 있습니다. 바로 지난 글에서 살펴본 비정규화된 수인 경우입니다. 이런 경우에는 지수부가 더 작아질 수 있으므로, 가수부에서 가장 먼저 1이 나오는 비트 위치를 찾아야 합니다.

for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;

이 반복문의 역할은 쉽게 말해 비정규화된 수가 가지는 실제 최대 지수 (-126)부터 시작해서, 가수부를 가장 높은 자리부터 보면서 처음으로 1이 나타나는 비트의 위치를 찾는 과정입니다. 여기서 신기한 테크닉을 사용하는 것을 볼 수 있는데, 조건문에서 i>0이라고 해놓고 증감문에서는 반대로 왼쪽 시프트를 해서 값을 점점 키워나가고 있다는 점입니다. 이 루프가 끝나는 시점은 가장 높은 자리에 있던 1이 부호 비트까지 올라가는 순간이 된다는 것을 이용한 테크닉입니다.

가수부 구하기

지수부를 구했으니 이번에는 가수부를 분리해내야 합니다. 단순히 hx, hy의 하위 23비트를 사용하면 되지 않겠느냐고 생각할 수 있으나, 여기서도 수가 정규화된 경우와 비정규화된 경우를 구분해야 하는 문제가 있습니다. 두 경우 모두 동일한 방법으로 사용할 수 있도록 통일성 있는 형태를 만들어 주는 것이 이 부분의 목적입니다. 또한, 기존의 가수부는 유효숫자의 맨 앞의 1을 제외하고 쓰였지만 연산 시에는 고려를 해줘야 하므로 그 수가 보이도록 가수부 앞의 비트가 1이 되게 만들어줍니다.²

x와 y 각각에 대해, 위에서 계산한 ix나 iy를 이용하여 정규화된 수인지 비정규화된 수인지를 판별합니다. 정규화된 수인 경우 우선 지수부의 비트들을 모두 지우고, 가수부 직전의 비트, 즉, 본래는 가수부의 가장 낮은 비트였던 자리를 1로 만들어주게 됩니다. 비정규화된 수인 경우에는 비트값이 1인 가장 높은 위치가 이 자리로 오도록 시프트를 해주는데, 이 간격이 바로 -126-ix, -126-iy입니다.

본격적인 연산

긴 사전 작업이 끝나고, 드디어 본격적인 연산에 들어가게 됩니다…만, 이 부분은 오히려 길지 않고 간단하게 수행되는 것을 볼 수 있습니다. 별도로 복잡한 과정 없이, 학교에서 배우는 그 나눗셈 방식을 그대로 사용합니다. 위에서의 사전 작업들을 통해 둘의 비트 구성을 동일하게 맞추어주었기 때문에 이런 단순화가 가능해집니다.

두 수의 자릿수 차이가 ix - iy이므로 이 횟수만큼 루프를 돌면서, hx가 hy 이상이라면 그 수를 빼고 자릿수를 조정해나가는 방식입니다. 이진수이기 때문에 몫의 각 자리 역시 0 또는 1이기 때문에 이렇게 두 가지 경우만으로 나누어 진행을 할 수 있습니다.

예시로 $1.11 \times 2^2$을 $1.10 \times 2^0$으로 나눈 몫을 구하면 다음과 같이 진행됩니다. 첫 번째 루프에서 $hz = 0.01$이 되고, 이는 $0$ 이상이므로 두 배를 하여 $hx = 0.10$이 됩니다. 두 번째 루프에서 $hz = -1.00$이 되므로 $hx$를 두 배를 하여 $hx = 1.00$이 됩니다. 이 값에 $2^{iy(=0)}$을 곱한 값이 실제 나머지의 값이 됩니다.

루프 중간에 hx가 0이 되면 hx가 hy의 배수라는 뜻이므로 나머지가 0이 됩니다. 이 경우 더 연산을 수행하지 않고 곧바로 ±0을 반환해서 불필요한 연산의 수를 줄이는 최적화를 수행합니다.

float 형태로 복원하기

연산 전 사전 작업이 많았던 것처럼, 연산 후 원래대로 되돌리는 것 역시 할 일들이 조금 있습니다. 다행히도(?) 이전만큼 복잡하지는 않고, 수를 하나만 변환하면 되기 때문에 훨씬 간단합니다. 여기서도 몇 가지 분기를 나누어가며 예외로 처리해야 할 사항들이 있습니다.

우선 hx가 0인 경우에 대한 예외 처리를 합니다. 위에서 루프 도중이 아닌 루프가 끝나는 순간이 된 것으로 역시 똑같이 ±0을 반환해 줍니다.

그 다음에는 hx를 정규화하는데, 가장 높은 1인 비트가 가수부의 가장 낮은 자리까지 올라올 때까지 한 칸씩 왼쪽으로 옮기고 iy를 그에 맞게 조정해 줍니다.

그 결과 만들어진 지수가 float의 정규화 범위에 해당되면 정규화에 맞게 지수부를 설정해주고, 비정규화 범위에 해당되면 비정규화에 맞게 설정해주게 됩니다. 마지막으로, 여기에 처음에 따로 빼놓았던 부호 비트 sx를 다시 추가해주고 float 레지스터에 옮기는 것으로 모든 작업이 끝나게 됩니다.

마치며

학교에서, 실생활에서 사용해오던 표기법과는 조금 다른 표기법인 부동소수점 표기법을 이용하여 연산하는 방법을 공부하는 것은 상당히 새로운 경험이었을 것입니다. 특히, 이러한 구조를 10진수가 아닌 2진수 구조로, 비트 단위로 연산하는 과정을 보는 것이 쉽지는 않았을 것입니다.

이 글에서는 먼저 가장 기본적인 연산들은 사칙 연산들을 어떻게 수행할 수 있는지를 따라가 보았습니다. 그리고 수학 라이브러리의 함수들 중 가장 기본적인 것들을 몇 가지 분석해보며, 이 함수들은 단순하게 한 두 문장으로 이루어지는 것이 아닌, 많은 양의 변환 과정과 단계적인 연산들을 거쳐 수행되는 복잡적인 연산이라는 것을 볼 수 있었습니다. 또한 비트 단위로 미리 정해진 필드들을 세밀하게 분리하고 각각의 의미를 정확하게 해석해야 올바른 식을 세울 수 있다는 것, 많은 예외 처리가 필요하다는 것 등도 볼 수 있었습니다.