Java實作伽馬函數gamma()

實習的項目裡遇到了第一類貝塞爾函數，C++裡有直接的貝塞爾函數可以使用，但是是實數域的而且僅C++17支援。項目裡是要用到複數域的半整階的第一類貝塞爾函數和第二類漢克爾函數，是以就直接根據原始定義實作了半整階的第一類貝塞爾函數。貝塞爾函數定義公式可以用伽馬函數求解，就順便把伽馬函數也實作了。寫完了發現C++的cmath裡有現成的伽馬函數tgamma()有float、double、long double三種形式。

思路分析

幾個特殊的伽馬函數值

Γ(0) = 1

Γ(1) = 1

Γ(1/2) = √π

将這些特殊值作為遞歸結束的條件。

在x > 1時，使用遞推公式Γ(x + 1) = x * Γ(x)推至0 < x < 1

在x < 0時，使用遞推公式Γ(x) = Γ(x + 1) / x推至0 < x < 1

在(0,1)區間内的伽馬值若無法遞推到特殊值，則使用伽馬函數原始定義公式實作。

注*：因為程式在計算時，會損失精度，導緻計算的結果不精确。故在運作時設定setAbsRelaErr=0.00001相對誤差絕對值來保護精度。

遞歸實作伽馬函數

public static double gamma(double x, double setAbsRelaErr) {
	//setAbsRelaErr 相對誤差絕對值
	//遞歸結束條件
	if(x < 0) {
		return gamma(x + 1, setAbsRelaErr) / x;
	}
	if(Math.abs(1.0 - x) < 0.00001) {
		return 1;
	}
	if(Math.abs(0.5 - x) < 0.00001) {
		return Math.sqrt(3.1415926);
	}
	
	if(x > 1.0) {
		return (x - 1) * gamma(x - 1, setAbsRelaErr);
	}
	
	double res = 0.0;
	double temp = 1.0;
	double check = 0.0;
	int i = 1;
	while(Math.abs((check - temp) / temp) > setAbsRelaErr){
		check = temp;
		temp *= i / (x - 1 + i);
		i++;
	}
	res = temp * Math.pow(i, x - 1);
	return res;
}

public static void main(String[] args) {
	System.out.println(gamma(-0.5, 0.00001));
}
           

循環實作伽馬函數

使用伽馬函數的遞推公式Γ(x + 1) = x * Γ(x)實作，傳進去的形參是x，根據公式在函數裡先減1，再根據公式計算。while循環的判斷條件是前後兩個數的內插補點不超過0.00001，以此來限制精度。跟遞歸實作相比，時間複雜度相當高。

public static double gamma(double x, double setAbsRelaErr) {
	x -= 1;
	double res = 0.0;
	int i = 1;
	double temp = 1.0;
	double temp0 = 1.0;
	double temp1 = 1.0;
	temp = 1 / (x + 1);
	do {
		temp0 = temp * Math.pow(i, x);
		i++;
		temp *= i / (x + i);
		temp1 = temp * Math.pow(i, x);
	}while(Math.abs(temp1 - temp0) > 0.00001);
	res = temp1;
	return res;
}