다항식

이 단원에서는 다항식들의 값을 계산하거나 푸는 함수들을 기술합니다. 해석적 방법을 통해, 2차 3차 식에 대해 실수, 허수 근을 찾는 기능들을 제공합니다. 일반적인 실수 계수의 다항식들에 대해서도 반복적인 방법을 통해 근을 찾는 방법 또한 제공합니다.

이 함수들은 헤더 파일 gsl_poly.h 정의되어 있습니다.

다항식의 계산

다항식 값의 계산은 안전성을 위해 호너 방법(Horner’s method)을 사용합니다.

\[P(x) = c[0] + c[1] x + c[2] x^2 + \dotsb + c[len-1] x^{len -1}\]

HAVE_INLINE 가 정의되어 있으면 inline 함수가 사용됩니다.

double gsl_poly_eval(const double c[], const int len, const double x): 실수 값 x 대해, 실수 계수를 가지는 다항식의 값을 계산합니다.

gsl_complex gsl_poly_complex_eval(const double c[], const int len, const gsl_complex z): 허수 값 z 대해, 실수 계수를 가지는 다항식의 값을 계산합니다.

gsl_complex gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex z): 허수 값 z 대해, 허수 계수를 가지는 다항식의 값을 계산합니다.

int gsl_poly_eval_derivs(const double c[], const size_t lenc, const double x, double res[], const size_t lenres): 다항식과 그 도함수를 계산해 길이 lenres 의 배열 res 에 저장합니다. 이 배열은 주어진 x 값에 대한, \(d^k P(x)/ d x^k\) 값을 각각의 원소로 가지고 있습니다. \(k\) 는 \(0\) 부터 시작합니다.

다항식의 분할 차분 표현

이 단락에서 기술하는 함수들은 뉴턴의 분할 차분법으로 표현된 다항식들을 다룹니다. 분할 차분법은 Abramowitz & Stegun의 25.1.4와 25.2.26 단원, 그리고 Burden and Faires, 3단원에 기술되어 있습니다. 이 단락에서는 분할 차분법의 개요를 서술합니다.

주어진 함수 \(f(x)\) 와 \(n\) 차수의 다항식 \(P_n(x)\) 는 함수 \(f\) 와 \(n+1\) 개의 서로 다른 지점 \(x_0, x_1, \dots, x_n\) 들로 표현할 수 있습니다. 다음과 같은 다항식의 표현법을 뉴턴의 분할 차분 표현이라 합니다.

\[P_n(x) = f(x_0) + \sum_{k=1}^n [x_0, x_1, \dots, x_k] (x-x_0)(x-x_1) \cdots (x-x_{k-1})\]

분할 차분 \([x_0, x_1, \dots, x_k]\) 은 Abramowitz & Stegun 25.1.4에 정의되어 있습니다. 추가적으로 \(x_0, x_1, \dots, x_k\) 에서 \(f\) 의 \(1\) 계 도함수 값이 주어졌다면, 차수 \(2n+1\) 의 보간 다항식을 만드는 데 사용할 수도 있습니다. 이 보간 다항식은 에르미트 보간 다항식으로 불리며 다음과 같이 정의됩니다.

\[H_{2n+1} (x) = f(z_0) + \sum_{k=1}^{2n+1} [z_0, z_1, \dots, z_k] (x-z_0)(x-z_1) \cdots (x-z_{k-1})\]

\(z = \{ x_0, x_0, x_1, x_1, \dots, x_n, x_n \}\) 는 \(z_{2k} = z_{2k+1} = x_k\) 로 정의됩니다. 분할 차분 \([z_0, z_1, \dots, z_k]\) 는 Burden and Faires의 3.4 단원을 참고할 수 있습니다.

int gsl_poly_dd_init(double dd[], const double xa[], const double ya[], size_t size): 길이 size 의 배열 xa 와 ya 에 저장된 지점 \((x,y)\) 에 대해, 보간 다항식의 분할 차분 표현을 계산합니다. ( xa ya )의 분할 차분 계산 결과는 길이 size 배열 dd 저장됩니다. 개요에서 설명한 서술법을 따라 \(dd[k] = [x_0, x_1,\dots, x_k]\) 입니다.

double gsl_poly_dd_eval(const double dd[], const double xa[], const size_t size, const double x): 길이 size 의 배열 xa 와 ya 에 저장된 분할 차분 표현된 다항식의 값을 주어진 변수 x 대해 계산합니다. 매크로 HAVE_INLINE 가 정의되어 있으면 인라인 버전이 사용됩니다.

int gsl_poly_dd_taylor(double c[], double xp, const double dd[], const double xa[], size_t size, double w[]): 분할 차분 표현된 다항식을 테일러 전개로 바꾸어줍니다. 분할 차분 표현은 길이 size 배열 dd 와 xa 로 표현됩니다. xp 지점의 테일러 계수들은 배열 c 에 저장됩니다. c 배열도 size 크기의 길이를 가집니다. 배열 w 는 길이가 size 로 같습니다.

int gsl_poly_dd_hermite_init(double dd[], double za[], const double xa[], const double ya[], const double dya[], const size_t size): 길이 size 배열 xa 와 ya 에 저장된 지점 \((x,y)\) 들에 대해, 에르미트 보간 다항식의 분할 차분 표현을 계산합니다. 에르미트 보간법으로 만들어지는 다항식은 \(1\) 계 도함수 \(dy/dx\) 의 값을 필요로 합니다. 이 값은 길이 size 의 배열 dya 로 주어져 있습니다. \(1\) 계 도함수 값들은 새로운 자료 집합 \(z= \{ x_0, x_0, x_1, x_1 \dots\}\) 를 정의해서, 일반적인 분할 차분법에 통합시킬 수 있습니다. 이 값들은 길이 \(2 \cdot\) size 의 배열 za 에 저장되어 있습니다. 계산 결과들은 \(2 \cdot\) size 길이를 가지는 배열 dd 저장됩니다. 개요에서 설명한 서술법을 따라 \(dd[k] = [z_0, z_1, \dots, z_k]\) 로 표현됩니다. 계산된 에르미트 다항식은 gsl_poly_dd_eval() 함수를 호출해 xa 에 대한 za 값을 넘겨 계산될 수 있습니다.

2 차 다항식 (Quadratic Equations)

int gsl_poly_solve_quadratic(double a, double b, double c, double *x0, double *x1)

2차 다항식

\[a x^2 + b x + c = 0\]

의 실수 근을 찾습니다. 근의 갯수(0,1,2)를 반환하며, 각 근의 위치는 x0 와 x 에 저장됩니다. 만약, 실수 근이 존재하지 않는다면 x0 와 x1 의 값을 수정하지 않습니다. 한 개의 근만이 있는 경우(예를 들어 \(a=0\) )는 x0 에 저장됩니다. 두 개의 근이 존재하면 x0 와 x1 는 각각 오름차순으로 저장됩니다. 중근의 경우는 특별히 취급되지 않습니다. 예를 들어, \((x-1)^2 = 0\) 은 값이 같은 두 개의 근을 가지는 방정식으로 취급됩니다.

근의 갯수는 \(b^2 -4ac\) 의 부호로 판별됩니다. 배 정밀도의 계산에서 이 방법은 반올림과 소거 오차의 영향을 받으며, 다항식의 계수가 정확하지 않을 때 비슷한 오류를 가질 수 있습니다. 하지만, 작은 정수 계수를 가지는 다항식에서는 정확하게 계산할 수 있습니다.

int gsl_poly_complex_solve_quadratic(double a, double b, double c, gsl_complex *z0, gsl_complex *z1)

2차 다항식

\[az^2 + bz +c =0\]

의 복소수 근을 계산합니다.

함수의 반환 값은 복소수 근의 숫자를 의미합니다. ( \(1\) 이거나 \(2\) 입니다.) 각 근의 위치는 z0 z1 에 저장됩니다. 저장되는 순서는 오름차순으로 저장되고, 실수부를 우선으로 판정하고 그 다음 허수부의 크기를 기준으로 배열합니다. 만약 한 개의 실수 근만 존재하면 ( 예를 들어 \(a=0\) ) z0 에 저장됩니다.

3 차 다항식 (Cubic Equations)

int gsl_poly_solve_cubic(double a, double b, double c, double *x0, double *x1, double *x2)

최고 차항의 계수가 1인 3차 다항식

\[x^3 + a x^2 + b x + c = 0\]

의 실수 근을 계산합니다. 실수 근의 숫자 (1-3)을 반환합니다. 이 근들의 위치는 x0 , x1 그리고 x2 에 저장됩니다. 만약 한 개의 실수 근만이 존재한다면, x0 에 저장됩니다. 세 개의 근이 존재한다면, 오름차순으로 x0 , x1 그리고 x2 에 저장됩니다. 중근은 특별하게 취급하지 않습니다. 예로 \((x-1)^3 =0\) 인 경우, 같은 값을 가지는 세 개의 근을 가지는 것으로 취급됩니다.

2차 다항식의 경우와 같이, 유한한 정밀도로 인해 밀접한 실수 근들이 실수 축에서 복소수 평면으로 이동해 근의 숫자가 달라질 수 있습니다.

int gsl_poly_complex_solve_cubic(double a, double b, double c, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2)

3차 다항식

\[z^3 + a z^2 + b z + c = 0\]

의 복소수 근을 찾습니다. 복소수 근의 숫자를 의미합니다(항상 3 입니다). 각 근의 위치는 z0 , z1 그리고 z2 에 저장됩니다. 각 근은 오름차순으로 실수부를 우선 판정하고, 허수부를 판정해 결정합니다.

일반 다항식

일반적으로 2차, 3차 그리고 4차 다항식같은 특수한 경우를 제외하면, 다항식 근은 해석적으로 찾을 수 없습니다. 이 단원에서 서술하는 알고리즘은 이러한 고차 다항식들의 근들을 반복적인 방법을 이용해 근사적인 위치를 구해줍니다.

type gsl_poly_complex_workspace: 일반적인 다항식의 근들을 찾기 위한 인자들을 가진 작업 공간입니다.

gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n)

gsl_poly_complex_workspace 구조체를 할당합니다. 이 작업 공간은 \(n\) 개의 계수를 가지는 다항식을 푸는 함수 gsl_poly_complex_solve() 를 위한 공간입니다.

오류가 생기지 않는다면, 새로 할당된 gsl_poly_complex_workspace 를 가르키는 포인터를 반환하고, 오류가 생기면 NULL 포인터를 반환합니다.

void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w): 작업 공간 w rk 할당된 모든 메모리를 해제합니다.

int gsl_poly_complex_solve(const double *a, size_t n, gsl_poly_complex_workspace *w, gsl_complex_packed_ptr z)

일반 다항 함수

\[P(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_{n-1} x^{n-1}\]

의 근들을 계산합니다. 동반 행렬(companion matrix)의 균형-QR 차원 감소를 이용합니다 1 . 인자 n 는 계수 배열의 길이를 나타냅니다. 가장 높은 차수의 계수는 반드시 \(0\) 이 아니여야 합니다. 적절한 크기의 작업 공간 w 를 필요로 합니다. 총, \(n-1\) 개의 근들이 반환되며, 크기 \(2(n-1)\) 의 복소수 배열 z 에 실수부-허수부 순서로 반복되어 저장됩니다.

모든 근들을 찾으면, GSL_SUCCESS 값으 반환합니다. 만약 QR 차원 감소가 수렴하지 않으면, 오류 관리자가 호출되고 GSL_EFAILED 오류 값를 전달합니다. 유의할 점은 유한한 정확도로 인해, 높은 차수의 중첩근은 저하된 정확도를 가지는 여러개의 단일 근들로 반환됩니다. 이러한 고-중첩근은 다중 구조를 고려하는 특별한 알고리즘이 필요합니다 2 .

예시

다음의 예시는 \(P(x) = x^5 -1\) 다항식을 이용해 주어진 일반적인 다항식의 해를 찾는 기능을 보여줍니다.

이 다항식은 다음의 5개 근을 가집니다.

\[1, e^{2 \pi i /5}, e^{4 \pi i /5}, e^{6 \pi i /5}, e^{8 \pi i /5}\]

다음 코드는 이 근들을 찾아줍니다.

#include <stdio.h>
#include <gsl/gsl_poly.h>

int
main (void)
{
  int i;
  /* coefficients of P(x) =  -1 +    x^5  */
  double a[6] = { -1, 0, 0, 0, 0, 1 }    ;
  double z[10];

  gsl_poly_complex_workspace * w
      =   gsl_poly_complex_workspace_allo  c (6);

  gsl_poly_complex_solve (a, 6, w, z)    ;

  gsl_poly_complex_workspace_free (w)    ;

  for (i = 0; i < 5; i++)
    {
      printf ("z%d = %+.18f %+.  18f\n",
              i, z[2*i], z[2*i+1]);
    }

  return 0;
}

프로그램의 결과 값은 다음과 같습니다.

z0 = -0.809016994374947673 +0.587785252292473359
z1 = -0.809016994374947673 -0.587785252292473359
z2 = +0.309016994374947507 +0.951056516295152976
z3 = +0.309016994374947507 -0.951056516295152976
z4 = +0.999999999999999889 +0.000000000000000000

이 결과는 \(z_n = e^{2\pi n i /5}\) 의 해석적 값들과 일치합니다.

참고 문헌과 추가 자료

균형 QR 방법과 이 방법의 오차 분석은 다음의 논문들에 기술되어 있습니다.

R.S. Martin, G. Peters and J.H. Wilkinson, “The QR Algorithm for Real Hessenberg Matrices”, Numerische Mathematik, 14 (1970), 219-231.
B.N. Parlett and C. Reinsch, “Balancing a Matrix for Calculation of Eigenvalues and Eigenvectors”, Numerische Mathematik, 13 (1969), 293-304.
A. Edelman and H. Murakami, “Polynomial roots from companion matrix eigenvalues”, Mathematics of Computation, Vol.: 64, No.: 210 (1995), 763-776.

분할 차분법 식들은 다음 문헌들에 기반합니다.

Abramowitz and Stegun, Handbook of Mathematical Functions, Sections 25.1.4 and 25.2.26.
R. L. Burden and J. D. Faires, Numerical Analysis, 9th edition, ISBN 0-538-73351-9, 2011.

1

balanced-QR reduction

2

Zeng, Algorithm 835, ACM Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp 218-236 을 참고할 수 있습니다.