【C言語】「lgammal」関数でガンマ関数の自然対数を簡単計算！使い方と代替方法

ヘッダーファイルと宣言

この関数を使用するには、 <math.h> ヘッダーファイルをインクルードする必要があります。

#include <math.h>

lgammal 関数は次のように宣言されています。

long double lgammal(long double x);

引数 x は、計算対象となる値です。この関数は long double 型の値を返し、ガンマ関数の自然対数を返します。

機能

lgammal 関数は、以下の式を使用してガンマ関数の自然対数を計算します。

lgamma(x) = (x - 1) * log(x) - 0.5 * log(2 * π) + LanczosApproximation(x)

ここで、LanczosApproximation は、数値的に効率的な近似関数です。

引数の値域

lgammal 関数は、以下の引数の値域で使用できます。

x < 0: 負の値の場合、関数は -inf または nan (非数) を返します。
x = 0: この場合、関数は +inf を返します。
x > 0: 正の値の場合、関数はガンマ関数の自然対数を返します。

注意事項

負の整数引数の場合は、関数は予期しない結果を返す可能性があります。
lgammal 関数は、大きな引数の値に対しては精度が低下する可能性があります。

例

#include <math.h>

int main() {
  long double x = 5.0;
  long double result = lgammal(x);

  printf("lgamma(%Lf) = %Lf\n", x, result);

  return 0;
}

このプログラムは、5 のガンマ関数の自然対数を計算し、結果をコンソールに出力します。

lgammal: long double 型の引数を受け取り、ガンマ関数の自然対数を返します。
lgammaf: float 型の引数を受け取り、ガンマ関数の自然対数を返します。
gamma: ガンマ関数を計算します。

例 1：ガンマ関数の自然対数の計算

この例では、5 のガンマ関数の自然対数を計算します。

#include <math.h>

int main() {
  long double x = 5.0;
  long double result = lgammal(x);

  printf("lgamma(%Lf) = %Lf\n", x, result);

  return 0;
}

このプログラムを実行すると、以下の出力が得られます。

lgamma(5.000000) = 14.198477121475602

例 2：様々な引数に対する `lgammal` 関数の値

この例では、様々な引数に対する lgammal 関数の値をいくつか示します。

#include <math.h>

int main() {
  long double x;
  long double result;

  printf("x | lgamma(x)\n");
  printf("---+-------------\n");

  for (x = 1.0; x <= 10.0; x += 1.0) {
    result = lgammal(x);
    printf("%Lf | %Lf\n", x, result);
  }

  return 0;
}

x | lgamma(x)
---+------------
1.0 | 0.000000
2.0 | 0.693147
3.0 | 1.732051
4.0 | 2.831745
5.0 | 3.998672
6.0 | 5.198582
7.0 | 6.474907
8.0 | 7.836107
9.0 | 9.233679
10.0 | 10.693982

この例では、lgammal 関数を使用してガンマ分布の確率密度関数を計算します。

#include <math.h>

double gamma_pdf(double x, double alpha, double beta) {
  if (x < 0.0 || alpha <= 0.0 || beta <= 0.0) {
    return NAN;
  }

  return (alpha - 1.0) * log(x) - beta * log(x) + lgammal(alpha) - lgammal(alpha + beta);
}

int main() {
  double x = 2.5;
  double alpha = 3.0;
  double beta = 2.0;
  double probability;

  probability = gamma_pdf(x, alpha, beta);

  printf("Probability of x = %Lf for Gamma(alpha, beta) distribution: %Lf\n", x, probability);

  return 0;
}

Probability of x = 2.500000 for Gamma(3.000000, 2.000000) distribution: 0.276711

Stirling の近似式

大きな引数の値に対しては、Stirling の近似式を使用して lgammal 関数を近似することができます。この近似式は、以下の式で表されます。

lgamma(x) ≈ (x - 0.5) * log(x) - 0.5 * log(2 * π) + (x + 0.5) * log(x + 1/2) - 1/12 * x + 1/120 * x^2

この近似式は、lgammal 関数よりも高速に計算できますが、精度が低くなる可能性があります。

Lanczos 近似

より高い精度が必要な場合は、Lanczos 近似を使用することができます。この近似式は、以下の式で表されます。

lgamma(x) ≈ (x - 1) * log(x) - 0.5 * log(2 * π) + LanczosApproximation(x)

ここで、LanczosApproximation は、数値的に効率的な近似関数です。この関数は、math.h ヘッダーファイルで宣言されています。

Lanczos 近似は、Stirling 近似よりも精度が高く、lgammal 関数よりも高速に計算できます。

前計算された値のテーブル

頻繁に同じ引数に対して lgammal 関数を呼び出す場合は、前計算された値のテーブルを作成することで、計算時間を短縮することができます。このテーブルには、引数とそれに対応するガンマ関数の自然対数の値を格納します。

前計算された値のテーブルを使用するには、まずテーブルを作成する必要があります。これは、ループを使用して各引数に対して lgammal 関数を呼び出すことで行うことができます。テーブルを作成したら、引数の値を使用してテーブルからガンマ関数の自然対数の値を検索することができます。

この方法は、lgammal 関数を毎回呼び出すよりも高速ですが、メモリ使用量が多くなります。

特殊関数ライブラリ

GNU Scientific Library (GSL) や Special Functions (SPFUN) などの特殊関数ライブラリには、lgammal 関数の高精度実装が含まれている場合があります。これらのライブラリは、C 標準ライブラリの lgammal 関数よりも精度が高いかもしれません。

最適な代替方法の選択

使用する代替方法は、状況によって異なります。精度が重要でない場合は、Stirling 近似が最速のオプションとなります。より高い精度が必要な場合は、Lanczos 近似を使用します。頻繁に同じ引数に対して lgammal 関数を呼び出す場合は、前計算された値のテーブルを作成します。高精度が必要で、C 標準ライブラリの lgammal 関数よりも高速な計算が必要な場合は、特殊関数ライブラリを使用します。

負の引数の場合は、lgammal 関数と代替方法のいずれも予期しない結果を返す可能性があります。
上記の代替方法はすべて、浮動小数点誤差の影響を受けます。