Script Programming(스크립트 프로그래밍)

Visual Studio 2019 와 MSYS2 MinGW64 에서 테스트 된 소스입니다.

 혹시 MinGW 에서 컴파일되지 않으면

$ packman -S mpfr

명령으로 mpfr 라이브러리를 설치하고 컴파일하면 된다.

// Filename: calcGammaFn.c
//
// Compile: gcc -o calcGammaFn calcGammaFn.c -lmpfr -lgmp
// Execute: ./calcGammaFn
//     Or
// Compile: cl calcGammaFn.c /I. mpfr.lib
// Execute: calcGammaFn
//
// Date: 2021.01.28


#include <stdio.h>
#include <math.h>    // for log(10)
#include <mpfr.h>

int main()
{
  mpfr_t x; 
  int i, inex;

  mpfr_set_emin (-41);
  mpfr_init2 (x, 42);

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_gamma (x, x, MPFR_RNDZ);      
      mpfr_subnormalize (x, inex, MPFR_RNDZ);
      mpfr_dump (x);
  }

  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_gamma (x, x, MPFR_RNDZ);
      mpfr_printf("Gamma(%2d) = %2d! = %Rf\n", i, i - 1, x);
  }

  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_lngamma (x, x, MPFR_RNDZ);
      mpfr_printf("LogGamma(%2d) = Log(%2d!) = %Rf\n", i, i - 1, x);
  }

  printf("\n");

 double t10 = log(10.0);
 printf("log(10) = %f\n", t10);
 
  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_lngamma (x, x, MPFR_RNDZ);      
      inex = mpfr_div_d(x, x, t10,  MPFR_RNDZ);
      mpfr_printf("Log10Gamma(%2d) = Log10(%2d!) = %Rf\n", i, i - 1, x);
  }

  mpfr_clear (x);
  
  return 0;
}

/*
Output:
LogGamma( 1) = Log( 0!) = 0
LogGamma( 2) = Log( 1!) = 0
LogGamma( 3) = Log( 2!) = 0.693147
LogGamma( 4) = Log( 3!) = 1.791759
LogGamma( 5) = Log( 4!) = 3.178054
LogGamma( 6) = Log( 5!) = 4.787492
LogGamma( 7) = Log( 6!) = 6.579251
LogGamma( 8) = Log( 7!) = 8.525161
LogGamma( 9) = Log( 8!) = 10.604603
LogGamma(10) = Log( 9!) = 12.801827
LogGamma(11) = Log(10!) = 15.104413
LogGamma(12) = Log(11!) = 17.502308
LogGamma(13) = Log(12!) = 19.987214
LogGamma(14) = Log(13!) = 22.552164
LogGamma(15) = Log(14!) = 25.191221
LogGamma(16) = Log(15!) = 27.899271
LogGamma(17) = Log(16!) = 30.671860

log(10) = 2.302585

Log10Gamma( 1) = Log10( 0!) = 0
Log10Gamma( 2) = Log10( 1!) = 0
Log10Gamma( 3) = Log10( 2!) = 0.301030
Log10Gamma( 4) = Log10( 3!) = 0.778151
Log10Gamma( 5) = Log10( 4!) = 1.380211
Log10Gamma( 6) = Log10( 5!) = 2.079181
Log10Gamma( 7) = Log10( 6!) = 2.857332
Log10Gamma( 8) = Log10( 7!) = 3.702431
Log10Gamma( 9) = Log10( 8!) = 4.605521
Log10Gamma(10) = Log10( 9!) = 5.559763
Log10Gamma(11) = Log10(10!) = 6.559763
Log10Gamma(12) = Log10(11!) = 7.601156
Log10Gamma(13) = Log10(12!) = 8.680337
Log10Gamma(14) = Log10(13!) = 9.794280
Log10Gamma(15) = Log10(14!) = 10.940408
Log10Gamma(16) = Log10(15!) = 12.116500
Log10Gamma(17) = Log10(16!) = 13.320620
*/

Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다.  단지 exp 함수를 이용하는 부분에서는 exp 함수의 구현에 따라 오차가 더 있을 수 있다.

감마함수 Γ(x)의 정의

(i) \textrm{$\alpha > 0$일 때는}

         \Gamma (\alpha) = \int_0^\infty e^{-t} t^{\alpha - 1} \ dt

(ii) \textrm{$\alpha < 0$ 이고 $\alpha \ne$정수 일 때는}

   \textrm{$\alpha + k> 0$ 이 되는 최소의 양의 정수 $k$를 찾아서}

         \Gamma (\alpha) = \dfrac{\Gamma(\alpha + k)}{\alpha (\alpha + 1) \cdots (\alpha +k - 1)}

[감마함수 Γ(x)의 특징]

(1) \Gamma (\alpha + 1) = \alpha \cdot \Gamma(\alpha)

(2) \textrm{특히 $n$이 양의 정수일 때는 \ } \Gamma (n) = (n - 1)!

아래에서는 Ubuntu에서 gnuplot을 이용하여 감마함수 Γ(x)의 그래프를 단계적으로 완성해 나가는 과정을 보여준다. 과정은 모두 5단계로 나뉘어져 있다. 일단 그려주고 나서 옵션을 조금씩 변경 또는 추가 하면서 그림을 다시 그려주는 방식을 취한다. (그려진 그림을 다시 그리는 gnuplot 명령은 replot이다)

* 단계 1: 감마함수의 그래프를 일단 그린다. (옵션은 차후에 조절한다.)

* 단계 2: x좌표축, y좌표축에 레이블을 붙여주고, 그래프의 제목을 표시한다.

* 단계 3: 좌표를 읽기 쉽게 정수 간격으로 그리드를 만든다. (set grid 명령은 그리드를 ON 하는 명령)

* 단계 4: x축과 y축을 (0.8 두께의 청색 실선으로) 그려준다.

* 단계 5: 수직 방향의 점근선 부근에 잘못된 부분 처리하기  ('set samples 숫자' 명령 사용)

감마함수 Γ(x)의 정의

(i) \textrm{$\alpha > 0$일 때는}

         \Gamma (\alpha) = \int_0^\infty e^{-t} t^{\alpha - 1} \ dt

(ii) \textrm{$\alpha < 0$ 이고 $\alpha \ne$정수 일 때는}

   \textrm{$\alpha + k> 0$ 이 되는 최소의 양의 정수 $k$를 찾아서}

         \Gamma (\alpha) = \dfrac{\Gamma(\alpha + k)}{\alpha (\alpha + 1) \cdots (\alpha +k - 1)}

[감마함수 Γ(x)의 특징]

(1) \Gamma (\alpha + 1) = \alpha \cdot \Gamma(\alpha)

(2) \textrm{특히 $n$이 양의 정수일 때는 \ } \Gamma (n) = (n - 1)!

* Ubuntu의 KmPlot을 이용하여 감마함수 Γ(x)의 그래프를 그린 화면

    (함수식에 f(x) = gamma(x) 라고 적어주고, 'Create" 버튼을 클릭하면 그려진다.)

* KmPlot 의 메뉴에서 "File" -> "Export..." 를 선택하여 출력한 감마함수 Γ(x)의 그래프