컴퓨터는 삼각함수를 어떻게 계산하는가

주의: 이 포스트는 약간의 미적분학 지식이 있어야 이해하기 쉽습니다. 그래도 설명을 위해 약간의 증명은 들어가 있으니, 아는 부분은 스킵하셔도 됩니다.

우리는 언제나 직각과 직사각형이 편합니다. 인류는 직교좌표계를 발명했습니다. 그리고 컴퓨터 모니터의 픽셀 배열은 대부분 가로세로 수백~수천 개 픽셀로 이루어진 격자로 이루어져 있고, 이는 직교좌표계로 접근할 수 있습니다. 직교좌표계로 화면을 그리면 픽셀 하나하나를 관리하기가 제일 쉽고 자유로워서가 아닐까 싶습니다.

만약에 누군가가 사칙연산만 제공되는 엔진으로 게임을 개발한다고 생각합시다. 직교좌표계 위에 그려진 어떤 도형이 있는데 이 도형을 회전시켜야 한다고 합니다. 어떻게 하면 될까요?

사실 수학은 이미 답을 알고 있습니다. 어떤 픽셀의 좌표 \(\left ( x, y\right )\)에 대해 회전된 좌표 \(\left ( x′, y′\right )\)는
\[\begin{bmatrix}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}x′\\ y′\end{bmatrix}\]
다른 말로 \[\begin{cases}x′=x \cos \theta-y\sin \theta\\ y′=x \sin \theta+y \cos \theta \end{cases}\] 임이 자명합니다. 이걸 회전하기 전의 이미지의 모든 픽셀의 \(\left ( x, y\right )\)에 대해 한 번씩 해 주면 회전한 이미지가 나옵니다. 참 쉽죠?

사칙연산만 갖고 삼각함수를 계산하라고요?

컴퓨터가 어떻게 삼각함수를 계산하는지를 논하기로 했으니까, 컴퓨터가 할 수 있는 연산만을 이용해야겠습니다. 아주 기본적인 연산만 사용 가능한 프로세서에서 말입니다.

그럼 대체 삼각함수는 어떻게 근사할 수 있을까요? \(\sin\)과 \(\cos\)를 구하면 \(\tan\) 등은 나눗셈으로도 구할 수 있으니까 \(\sin\)과 \(\cos\)를 구하는 데 집중해 봅시다.

삼각함수를 다항함수로 근사하면 되지 않을까요?

다항함수는 사칙연산만으로 계산할 수 있으니 삼각함수의 그래프와 비슷한 다항함수를 만들어서 거기다 집어넣으면 되겠네요!

근데 그 다항함수는 대체 어떻게 구할까요? 여기서 테일러 급수가 등장합니다.

테일러 급수

테일러 급수를 이용하면 어떤 미분 가능한 함수라도 다항함수로 근사해 버릴 수 있습니다. 세상에 그런 흑마법이 존재하냐고요? 네, 놀랍게도 존재합니다!

증명

우선 어떤 함수 \(f^\prime\left(x\right )\)를 \(a\)부터 \(x\)까지 정적분해 봅시다. 미적분의 정의에 의해 아래 식과 같이 표현할 수 있습니다.
\[\int_{a}^{x}f^\prime\left(t\right )\mathrm{d}t=f\left(x\right )-f\left(a\right )\]
그리고 위 식을 변형해 봅시다.
\[\int_{a}^{x}f^\prime\left(t\right )\mathrm{d}t=\int_{a}^{x}\left(-1\right )\left(-f^\prime\left(t\right )\right )\mathrm{d}t\]
이제 \(-1\)을 적분하고 \(-f^\prime\left(t\right )\)를 미분해 부분적분법을 씁시다. 이 때 \(f\left(x\right )\)가 무한히 미분 가능하다면, 부분적분법도 무한히 써 버릴 수 있습니다. 그러니까
\[\cdots = \left.\begin{matrix}\left (-\left ( x-t \right )f^\prime\left ( t \right )- \dfrac{\left( x-t \right )^2}{2}f^{\prime\prime}\left ( t \right )- \dfrac{\left( x-t \right )^3}{6}f^{\prime\prime\prime}\left ( t \right )- \cdots\right )\end{matrix}\right|_{a}^{x}\]
가 됩니다. 이 때 적분변수 \(t\)와 관계없는 \(x\)는 상수취급됩니다. 이 식을 전개해 보면
\[\cdots =\left ( x-a \right )f^\prime\left ( a \right )+ \dfrac{\left( x-a \right )^2}{2!}f^{\prime\prime}\left ( a \right )+ \dfrac{\left( x-a \right)^3}{3!}f^{\prime\prime\prime}\left ( a \right )+ \cdots\\ \]
\[= f\left ( x \right )-f\left ( a \right )\]
마지막으로 \(f\left ( a \right )\)를 이항하면
\[f\left ( x \right ) =f\left ( a \right )+ \left ( x-a \right )f^\prime\left ( a \right )+ \dfrac{\left( x-a \right )^2}{2!}f^{\prime\prime}\left ( a \right )+ \dfrac{\left( x-a \right )^3}{3!}f^{\prime\prime\prime}\left ( a \right )+ \cdots\]
즉 \(f^{\left( n \right)}\)이 \(f\)를 \(n\)번 미분한 함수라고 할 때
\[f\left ( x \right ) = \sum_{n=0}^{\infty}\dfrac{f^{\left ( n \right )}\left ( a \right )}{n!}\left( x-a \right )^n\]
가 유도됩니다. 이 때 \(a\)에 어떤 수를 넣어도 근사됩니다.

그래서 이게 삼각함수랑 무슨 상관이에요!

삼각함수는 무한히 미분 가능합니다. 그러니까, \(f\) 대신 삼각함수를 넣으면 삼각함수를 다항함수로 근사할 수 있습니다. 어떤 삼각함수 \(f\left ( x \right )\)에 대해 \(a\)에 \(0\)을 대입해 보면
\[f\left ( x \right ) =f\left ( 0 \right )+ \left ( x-0 \right )f^\prime\left ( 0 \right )+ \dfrac{\left( x-0 \right )^2}{2!}f^{\prime\prime}\left ( 0 \right )+ \dfrac{\left( x-0 \right )^3}{3!}f^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
\[=f\left ( 0 \right )+ x f^\prime\left ( 0 \right )+ \dfrac{x^2}{2!}f^{\prime\prime}\left ( 0 \right )+ \dfrac{x^3}{3!}f^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
인데요, \(f\left ( x \right ) = \sin x\)라고 하고 대입해 보면
\[\sin x =\sin\left ( 0 \right )+ x \sin^\prime\left ( 0 \right )+ \dfrac{x^2}{2!}\sin^{\prime\prime}\left ( 0 \right )+ \dfrac{x^3}{3!}\sin^{\prime\prime\prime}\left ( 0 \right )+ \cdots\]
\[=0+ \left ( x \cdot 1 \right )+ \left ( \dfrac{x^2}{2!}\cdot 0 \right )+ \left ( \dfrac{x^3}{3!}\cdot -1 \right )+ \left ( \dfrac{x^4}{4!}\cdot 0 \right )+ \left ( \dfrac{x^5}{5!}\cdot 1 \right )+ \cdots\]
정리하면
\[\therefore \sin x =x- \dfrac{x^3}{3!}+ \dfrac{x^5}{5!}- \dfrac{x^7}{7!}+ \dfrac{x^9}{9!}+ \cdots\]
마찬가지로
\[\cos x =1-\dfrac{x^2}{2!}+\dfrac{x^4}{4!}-\dfrac{x^6}{6!}+\dfrac{x^8}{8!}+ \cdots\]
이 됩니다. 중간 과정에 비해서 정말 간단한 식이 되었습니다. 더구나, 우리의 궁극적인 목표인 사칙연산만으로 삼각함수 표현하기에 성공했습니다! (팩토리얼은 곱셈의 연속일 뿐이니까요!)

근데, 식이 참… 무한합니다. 이건 어떻게 하면 좋을까요? 우리는 어차피 floating-point 변수들은 굉장히 정확하지 않다는 걸 잘 알고 있습니다. 그러니까 저 식을 어디까지만 계산하고 그 이후의 오차는 무시해 버려도 됩니다.

위에서 구한 다항식을 최고차항이 \(n\)인 항까지만 계산하고, 그 이후는 무시해버립시다. 그리고 이걸 \(n\)차 근사식이라고 합시다. 이제부터 차수를 올려나가면서 \(n\)차 근사식의 그래프를 그려보겠습니다.

그래프

▶ 1차 근사식: \(f\left ( x \right ) = x\)

… 하나도 비슷해보이진 않습니다. 다음!

▶ 3차 근사식: \(f\left ( x \right ) = x-\dfrac{x^3}{3!}\)

앞에서는 감을 좀 잠은 거 같기도 한데, 겨우 \(\dfrac{\pi}{2}\)도 가기 전에부터 조금씩 이상하더니 아래로 곤두박질쳐버리는군요. 다음!

▶ 5차 근사식: \(f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}\)

이제 \(\dfrac{\pi}{2}\)에서도 비슷해졌네요! 다음!

▶ 7차 근사식: \(f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!}\)

\(pi\) 근처까지도 꽤 비슷해졌습니다. 몇 개만 더 해 봅시다.

▶ 9차 근사식: \(f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}-\dfrac{x^7}{7!} + \dfrac{x^9}{9!}\)

거의 원본 함수와 같아졌습니다. 사실 여기부터는 이제 \(0 < x < \pi\)일 때는 그냥 가져다 써도 오차는 크지 않을 것 같습니다. 어차피 \(sin\)은 주기함수이고 대칭함수이니까 \(0 < x < \dfrac{\pi}{2}\)에서만 정확해도 다른 범위에서는 함수의 특징을 이용해 계산해내면 됩니다.

▶ 11차 근사식: \(f\left ( x \right ) = x-\dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!}\)

차이가 얼마나 나는지 점점 이렇게 봐서는 확인하기가 어렵습니다. 마지막으로 13차 근사식까지만 그려보겠습니다.

▶ 13차 근사식: \(f\left ( x \right ) = x- \dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!} + \dfrac{x^{13}}{13!}\)

그래프 범위 내에서는 거의 똑같은 곡선이 그려집니다.

프로그래밍 언어에서의 구현 사례

OpenJDK의 StrictMath 네이티브 코드가 13차 근사식을 통해 삼각함수를 구현하고 있습니다. 정확성을 위해서 13차 근사식을 그대로 사용하진 않고, 이렇게 조금 변형해서 사용합니다.
\[\sin x \sim x- \dfrac{x^3}{3!} + \dfrac{x^5}{5!}- \dfrac{x^7}{7!} + \dfrac{x^9}{9!}- \dfrac{x^{11}}{11!} + \dfrac{x^{13}}{13!}\]
\[\frac{\sin x}{x} \sim 1- \dfrac{x^2}{3!} + \dfrac{x^4}{5!}- \dfrac{x^6}{7!} + \dfrac{x^8}{9!}- \dfrac{x^{10}}{11!} + \dfrac{x^{12}}{13!}\]
그러므로 \(r\)을 이렇게 정의할 때
\[r = x^3 \left(\dfrac{1}{5!} + x^2 \left(-\dfrac{1}{7!} + x^2 \left(\dfrac{1}{9!} + x^2 \left(-\dfrac{1}{11!} + {x^{2}} \cdot \dfrac{1}{13!}\right)\right)\right)\right)\]
\(\sin x\)는 이렇게 표현할 수 있습니다.
\[\sin x \sim x + x^3 \cdot \left(-\dfrac{1}{3!} + x^2 r \right)\]
예를 봅시다. sin 함수에서 __kernel_sin 함수를 호출할 때 코드에서 \(|x| \prec \pi/4\)라면 y, iy의 값은 모두 0입니다. 이외의 범위에서는 적절한 범위로 평행이동시켜 __kernel_sin 혹은 __kernel_cos 함수값을 구합니다.

#include "fdlibm.h"
#ifdef __STDC__
    static const double
#else
    static double
#endif
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */

#ifdef __STDC__
    double __kernel_sin(double x, double y, int iy)
#else
    double __kernel_sin(x, y, iy) double x,y; int iy; /* iy=0 if y is zero */
#endif
{
    double z,r,v; int ix;
    ix = __HI(x)&0x7fffffff; /* high word of x */
    if(ix < 0x3e400000) /* |x| < 2**-27 */ {
        if((int)x==0) return x;
    } /* generate inexact */
    z = x*x;
    v = z*x;
    r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
    if(iy==0) return x+v*(S1+z*r);
    else return x-((z*(half*y-v*r)-y)-v*S1);
}

소스에서
\(S1 = -\dfrac{1}{3!} \approx -1.66666667 \times {10}^{-1}\)
\(S2 = \dfrac{1}{5!} \approx 8.33333333 \times {10}^{-2}\)
\(S3 = -\dfrac{1}{7!} \approx -1.98412698 \times {10}^{-4}\)
\(S4 = \dfrac{1}{9!} \approx 2.75573137 \times {10}^{-6}\)
\(S5 = -\dfrac{1}{11!} \approx -2.50521084 \times {10}^{-8}\)
\(S6 = \dfrac{1}{13!} \approx 1.58969099 \times {10}^{-10}\)
으로 근사식의 계수들이 근사되어 하드코딩되어 있는 것을 알 수 있습니다. 또한 \(z=x^2\), \(v=zx=x^3\)으로 계산하고 있습니다. 13차 근사식을 그대로 쓴다면 항마다 1번, 3번, 5번, …, 13번의 곱셈을 해야 하지만 \(x^2\)를 새 상수로 두면 곱셈을 줄일 수 있어서 이런 방법을 사용하지 않았을까 하는 생각입니다.

다른 방법으로도 계산할 수 있나요?

물론 다른 방법들도 있습니다. 예를 들어 볼더가 1956년에 고안한 CORDIC(COordinate Rotation DIgital Computer) 알고리즘을 이용해 삼각함수의 값을 근사하는 것도 가능합니다.

간단하게 설명하자면, 맨 위에서 회전을 하기 위해 곱하는 행렬을
\[\begin{bmatrix}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta\end{bmatrix}\]
으로 소개했습니다만, CORDIC 알고리즘은 이를 \(\tan\) 함수만으로 표현해
\[R_i = \dfrac{1}{\sqrt{1 + \tan^2 \left(\gamma_i\right)}}\begin{bmatrix}1 & -\tan \left(\gamma_i\right) \\\tan \left(\gamma_i\right) & 1\end{bmatrix}\]
로 변형하고, 컴퓨터가 빠르게 계산할 수 있도록 \(\tan \left(\gamma_i\right) = \pm2^{-i}\)가 되는 \(\gamma_i\), 즉 \(\arctan\left(2^{-i}\right)\)의 값들과 이 때의 \(\dfrac{1}{\sqrt{1 + \tan^2 \left(\gamma_i\right)}}\)의 값들을 하드코딩해 두고 원하는 각도에 가까워질 때까지 \(i\)를 하나씩 증가시키면서
\(\arctan\left(2^{-i}\right)\)를 더하거나 빼면서 행렬 계산을 해나가는 알고리즘입니다.

이 알고리즘은 계산기에 주로 쓰이고 있고, 인텔의 8087 ~ 80486 CPU에도 채용되어 왔습니다.

컴퓨터는 이렇게 삼각함수를 계산합니다.