아래의 소스는
Coefficients for the Lanczos Approximation to the Gamma Function
의 것을 Microsoft Visual C/C++ 명령줄 컴파일러 cl 로 컴파일되도록 수정한 것이다.
// Filename: test_gamma_lanczos.c
//
// Compile: cl test_gamma_lanczos.c /EHsc /utf-8
// Execute: test_gamma_lanczos
#include <stdio.h>
#define _USE_MATH_DEFINES // for C
#include <math.h>
#define LG_g 5.0 // Lanczos parameter "g"
#define LG_N 6 // Range of coefficients i=[0..N]
const double lct[LG_N+1] = {
1.000000000190015,
76.18009172947146,
-86.50532032941677,
24.01409824083091,
-1.231739572450155,
0.1208650973866179e-2,
-0.5395239384953e-5
};
const double ln_sqrt_2_pi = 0.91893853320467274178;
const double g_pi = 3.14159265358979323846;
double lanczos_ln_gamma(double z)
{
double sum;
double base;
double rv;
int i;
if (z < 0.5) {
return log(g_pi / sin(g_pi * z)) - lanczos_ln_gamma(1.0 - z);
}
z = z - 1.0;
base = z + LG_g + 0.5; // Base of the Lanczos exponential
sum = 0;
for(i=LG_N; i>=1; i--) {
sum += lct[i] / (z + ((double) i));
}
sum += lct[0];
return ((ln_sqrt_2_pi + log(sum)) - base) + log(base)*(z+0.5);
}
double lanczos_gamma(double z)
{
return(exp(lanczos_ln_gamma(z)));
}
int main()
{
printf(" Gamma(0.5) Log(Gamma(0.5))\n");
printf(" ---------------------------------------------------\n");
printf(" Lanczos Gamma: %22.18f %22.18f\n", lanczos_gamma(0.5), lanczos_ln_gamma(0.5));
printf(" Built-in: %22.18f %22.18f\n", exp(lgamma(0.5)), lgamma(0.5));
printf(" Use sqrt(PI): %22.18f %22.18f\n", sqrt(M_PI), log(sqrt(M_PI)));
return 0;
}
빌트인 Gamma 함수와 비교하여 12자리 까지만 맞는데,
이는 허용오차가 약 sqrt(PI)*N 이기 때문이다.
여기서 N은 Lanczos 계수의 개수이다, (소스에서는 LG_N+1)
Python 인터프러터를 실행하여 잠깐 그 허용오차를 알아본다.
Python 3.7.5
>>> import math
>>> math.gamma(0.5)
1.7724538509055159
>>> import math
>>> math.sqrt(math.pi)*7
12.407176956338612
N=7 인 경우, 유효수자의 개수는 약 12개이다.
'프로그래밍 > C' 카테고리의 다른 글
C 언어에서 큰 부동소수점수(native double) 의 정확도 (0) | 2023.03.19 |
---|---|
MSYS2 에서 C 언어의 printf(...) 가 제대로 출력하지 못할 때... (0) | 2022.01.08 |
MPFR 라이브러리를 이용하여 Gamma 함수값 계산하기 (0) | 2021.01.28 |
cygwin 의 gcc 로 UTF-8 한글 처리하는 간단한 예제 (0) | 2014.04.13 |
cygwin/mingw 의 gcc 로 utf-8 한글 처리하기 (0) | 2014.04.02 |