컴퓨터 속 삼각함수 - keompyuteo sog samgaghamsu

중3때 삼각비라는 부분에서 sin, cos, tan을 배웁니다. 특수각인 0°, 30°, 45°, 60°, 90°에서의 sin, cos, tan값을 배우죠. 고2에 오면 이과 과정 미적분 II에서 삼각함수 단원에서 좀 더 자세히 배웁니다. 60분법 대신에 호도법을 사용하며 정의역도 실수 전체로 확장하죠. 삼각함수 덧셈, 뺄셈정리를 통해 15°, 75° 등의 각을 삼각함수의 특수각 값을 이용하여 구하죠. 하지만 현실에서 각도는 특수각의 합, 차, 1/2만 존재하지 않습니다. 6°, 34.5° 등의 일반적인 각에서 우리는 삼각함수의 값을 구할 필요가 있을 때가 있죠. 그러나 우리가 배운 내용으로는 이러한 값을 구할 수 없습니다. 그러나 컴퓨터는 이러한 값들을 아주 손쉽게 계산해냅니다. 어떻게 계산해내는 걸까요?

그 과정을 알아보기 위하여 C언어의 math.h에서 사용하는 삼각함수 라이브러리(glibc)를 뜯어보았습니다. C언어에서는 math.h를 include하여 sin, cos, tan 등의 삼각함수를 사용할 수 있습니다. 이 링크에서 아래 코드의 원본을 보실 수 있습니다.

/* k_sinf.c -- float version of k_sin.c
 * Conversion to float by Ian Lance Taylor, Cygnus Support, .
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */

#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
#endif

#include <math.h>
#include <math_private.h>

static const float 
half =  5.0000000000e-01,/* 0x3f000000 */
S1  = -1.6666667163e-01, /* 0xbe2aaaab */
S2  =  8.3333337680e-03, /* 0x3c088889 */
S3  = -1.9841270114e-04, /* 0xb9500d01 */
S4  =  2.7557314297e-06, /* 0x3638ef1b */
S5  = -2.5050759689e-08, /* 0xb2d72f34 */
S6  =  1.5896910177e-10; /* 0x2f2ec9d3 */

float __kernel_sinf(float x, float y, int iy)
{
    float z,r,v;
    int32_t ix;
    GET_FLOAT_WORD(ix,x);
    ix &= 0x7fffffff;           /* high word of x */
    if(ix<0x32000000)           /* |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);
}

이 kernel_sin을 사용하여 실제의 sin함수를 구현하는 부분인 s_sinf.c 코드를 보면 $-{\pi \over 4} < x < {\pi \over 4}$ 에서는 간단히 iy = 0으로 정의되어 있을 것을 볼 수 있으니 iy = 0일 때 알고리즘을 살펴봅시다. 크기가 크지 않으니 z, v, r을 계산하여 x + v * (S1 + z * r)을 계산하여 넘기는 것으로 볼 수 있습니다. 이를 식으로 정리해보면 아래와 같습니다.

$$ \sin x = x + S1 x^3 + S2 x^5 + S3 x^7 + S4 x^9 + S5 x^{11} + S6 x^{13} $$

이제 위에 하드코딩된 값을 대입하여 이 그래프를 그려보겠습니다. Legend에 나와있는 것처럼 f(x)가 위의 함수고, sin(x)는 Mathmatica 기본 내장 함수입니다.

사실상 일치해서 sin(x) 값이 다 덮어씌운 것만 보입니다. 이제 정의역을 f(x)에 정의된 것보다 좀 더 늘려볼까요?

정의역을 $-2\pi < x < 2\pi$ 까지 늘려봤습니다. 이제는 어느정도 차이가 두드러지게 나타나네요! 실제 sin 함수를 구현하는 s_sinf.c 코드를 보시면 $-{\pi \over 4} < x < {\pi \over 4}$ 이외에는 오류처리거나 reduction을 통해 표현하는 것을 볼 수 있습니다. 왜냐면 sin은 $0 \le x \le {\pi \over 4}$ 만 알면 정의역이 모든 실수일 때 정의할 수 있기 때문입니다.

위와 같이 다항식으로 해석함수(고등학교 과정에서는 무한히 미분 가능한 함수 정도? 자세한 건 맨 아래에 있습니다)를 표현하는 것을 테일러 급수(Taylor series)라고 합니다. 이제 테일러 급수가 어떻게 미분 가능한 함수를 표현하는지 알아봅시다. 일단 미분 가능한 함수에 대한 n차 근사식을 $p_n(x)$ 라고 합시다. 그리고 우리는 한 점 $\left(x_0,\ f(x_0)\right)$ 과 그 점에 대한 미분계수( $f'(x)$ )가 주어졌을 때 그 점에서의 접선을 나타낼 수 있습니다 공식은 아래와 같죠.

$$ p_1(x) = f(x_0) + f'(x_0)(x - x_0)$$

이 공식에서 우리는 아래와 같은 성질을 찾을 수 있습니다.

$$ p_1(x_0) = f(x_0),\ p_1'(x_0) = f'(x_0) $$

따라서 이와 같은 성질을 이용하여 이번엔 2차 다항식에서 근사식을 만들어봅시다. 위의 성질을 사용하면 아래와 같은 조건을 만족하는 $p_2(x)$ 를 구해야겠죠.

$$ p_2(x_0) = f(x_0),\ p_2'(x_0) = f'(x_0),\ p_2''(x_0) = f''(x_0) $$

이에 대해 만족하는 식을 구하면 아래와 같습니다.

$$ p_2(x) = f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 $$

이와 같은 방식을 계속 반복하면 n차 미분한 식을 $f^n(x)$ 라고 했을 때 아래와 같은 n차 근사식을 구할 수 있습니다.

$$ p_n(x) = \sum^n_{k=0}{f^k(x - x_0)^k \over k!} $$

이때 아래와 같은 표현으로 $f(x)$ 를 나타내는데, 이 표현을 테일러 정리라고 하며 $R_n$ 을 나머지 항이라고 부릅니다.

$$ \begin{align*} f(x) &= p_n(x) + R_n \\&= f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 + \cdots + {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1} \end{align*} \\ R_n(x) = {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1}$$

나머지항은 일종의 테일러 정리의 식과 실제 $f(x)$ 의 식 간의 차이(오류)를 보여주므로, $R_n(x)$ 가 $n = \infty$ 일 때 0으로 수렴하면 테일러 정리의 식과 실제 식과 같아지므로 이때 우리는 위 무한급수를 테일러 급수라고 하며, 특별히 $x_0 = 0$ 인 경우에는 매클로린 급수라고 합니다.

이러한 방식으로 테일러 급수를 유도해낼 수 있는데, 수학적으로 엄밀하게 증명해보려면 미분법과 적분법을 사용하는 방법이 있습니다. 굳이 쓸 필요는 없어 보여서 이 링크를 참고하시기 바랍니다. 검색해도 잘 나와요.

자, 이제 대략적으로 테일러 급수에 대한 내용을 정리해보았습니다. 이제 그러면 적용을 해봅시다. 위에서 예시를 든 것처럼 $\sin x$ 에 대해서 테일러 급수를 써봅시다. 아래 그래프는 $\sin x$ 의 그래프(검정)와, $x_0 = 0$ 에서 테일러 급수의 1차(빨강), 3차(주황), 5차(노랑), 7차(초록), 9차(파랑), 11차(보라) 항까지 합한 것의 그래프입니다.

보이는 것처럼 합하는 항수가 커질수록 실제 $\sin x$ 에 다가가는 것을 볼 수 있습니다.

이제 다음으로 파데 근사(Padé approximant)에 대해 써보려고 합니다. 어떤 함수를 다항식으로 표현한 것이 테일러 급수라면, 어떤 함수를 유리식으로 표현한 것이 파데 근사입니다. 테일러 급수의 일반화이며 상위호환이라고 하는 거 같네요. 음이 아닌 정수 $m,\ n$에 대하여 파데 근사의 꼴은 아래와 같습니다. 이때 $n = 0$ 이 되면 매클로린 급수와 같게 됩니다.

$$ R(x)_{m/n} = {\sum_{j=0}^m{a_jx^j} \over 1 + \sum_{i=1}^n{b_kx^k}} $$

...파데 근사 계산을 한 번 해서 정리해보고 싶었는데 내용 이해가 잘 안 되는 관계로(쿨럭) 링크는 남겨둘테니 궁금하신 분들은 여기서 살펴보세요. 계산법을 제대로 이해하긴 힘들거 같고, 대신에 Mathematica에서 PadeApproximant를 함수로 제공하기 때문에 이것으로 $x_0=0$ 일 때 계산한 $\sin x$ 를 아래 표로 정리해보았습니다.

m\n	1(빨강)	3(주황)	5(초록)	7(파랑)
0(0)	$x$	$x-\frac{x^3}{6}$	$x-\frac{x^3}{6}+\frac{x^5}{120}$	$x-\frac{x^3}{6}+\frac{x^5}{120}-\frac{x^7}{5040}$
2(0.3)	$\dfrac{x}{1 + \frac{x^2}{6}}$	$\dfrac{x - \frac{7x^3}{60}}{1 + \frac{x^2}{20}}$	$\dfrac{x - \frac{x^3}{7} + \frac{11x^5}{2520}}{1 + \frac{x^2}{42}}$	$\dfrac{x - \frac{11x^3}{72} + \frac{13x^5}{2160} - \frac{x^7}{12096}}{1 + \frac{x^2}{72}}$
4(0.7)	$\dfrac{x}{1 + \frac{x^2}{6} + \frac{7x^4}{360}}$	$\dfrac{x - \frac{31x^3}{294}}{1 + \frac{3x^2}{49}+\frac{11x^4}{5880}}$	$\dfrac{x - \frac{53x^3}{396} + \frac{551x^5}{166320}}{1 + \frac{13x^2}{396} + \frac{5x^4}{11088}}$	$\dfrac{x - \frac{241x^3}{1650} + \frac{601x^5}{118800} - \frac{121x^7}{2268000}}{1 + \frac{17x^2}{825} + \frac{19x^4}{118800}}$