# 第1種ベッセル関数(漸化式) Bessel Function Of The First Kind
# 引数 次数 変数 ($N, $X)
# 戻り値 第1種ベッセル関数 ($BesselFunction)
sub BESSELFUNCTIONFIRST{
	my ($N, $X) = @_;
	my $BesselFunction = 0;
	my $NextBessel = 0;
	my $Bessel = 0;
	my $PrevBessel = 0;
	my $Diff = $N - int($N);
	my $Add = 0;
	my $Count = abs($N);

	if((-2 < $N) && ($N < 2)){
		# 第1種ベッセル関数(展開) Bessel Function Of The First Kind
		$Bessel = &BESSELFUNCTIONEXPANSION($N, $X);
	}else {
		if(2 <= $N){
			$Add = 1;
		}else {
			$Add = -1;
		}

		# 第1種ベッセル関数(展開) Bessel Function Of The First Kind
		$PrevBessel = &BESSELFUNCTIONEXPANSION($Diff, $X);
		$Bessel = &BESSELFUNCTIONEXPANSION($Diff + $Add, $X);

		for(my $i = $Diff + $Add; abs($i) < $Count; $i += $Add){
			$NextBessel = (((2 * $i) / $X) * $Bessel) - $PrevBessel;

			# 漸化式
			$PrevBessel = $Bessel;
			$Bessel = $NextBessel;
		}
	}

	# 第1種ベッセル関数 Bessel Function Of The First Kind
	$BesselFunction = $Bessel;

	return $BesselFunction;
}

# 第1種ベッセル関数(展開) Bessel Function Of The First Kind
# 引数 次数 変数 ($N, $X)
# 戻り値 第1種ベッセル関数(展開) ($BesselFunction)
sub BESSELFUNCTIONEXPANSION{
	my ($N, $X) = @_;
	my $BesselFunction = 0;
	my $Bessel = 0;
	my $PrevBessel = 0;
	my $Num1 = 0;
	my $Num2 = 0;
	my $Den1 = 0;
	my $Den2 = 0;
	my $Diff = $N - int($N);
	my $AbsN = ($Diff == 0 ? abs($N) : $N);
	my $Limit = 100;
	my $Epsilon = 1.0e-10;
	
	for(my $i = 0; $i < $Limit; $i++){
		# 分子
		$Num1 = (($i % 2) == 0 ? 1: -1);
		$Num2 = ($X / 2) ** ((2 * $i) + $AbsN);
		# 分母
		$Den1 = ($i == 0 ? 1 : ($Den1 * $i));
		$Den2 = &GAMMAFUNCTION($AbsN + $i + 1);
	
		# 一つ前
		$PrevBessel = $Bessel;
		# 第1種ベッセル関数 Bessel Function
		$Bessel += (($Num1 * $Num2) / ($Den1 * $Den2));
	
		# 収束判定
		last if(abs($Bessel - $PrevBessel) < $Epsilon);
	}
	
	# 負の整数の場合
	if(($N < 0) && ($Diff == 0)){
		$Bessel = $Bessel * ((-1) ** $AbsN);
	}
	
	# 第1種ベッセル関数 Bessel Function Of The First Kind
	$BesselFunction = $Bessel;
	
	return $BesselFunction;
}
	
# ガンマ関数(近似) Gamma Function Approximation
# 引数 値 ($X)
# 戻り値 ガンマ関数(近似) ($GammaFunction)
sub GAMMAFUNCTION{
	my ($X) = @_;
	my $GammaFunction = 0;
	my $Factorial = 1;
	my $Diff = abs($X - int($X));
	my $Pi = atan2(1, 1) * 4;
	my $Temp = 0;
	
	# 値の確認
	if(($X <= 0) && ($Diff == 0)){
		return 0;
	}
	
	if($Diff == 0){
		for(my $i = $X - 1; $i >= 2; $i--){
			$Factorial *= $i;
		}
	
		# X-1の階乗 ガンマ関数 Gamma Function
		$GammaFunction = $Factorial;
	}elsif($Diff == 0.5){
		my $Start = ($X > 0 ? ((2 * int($X)) - 1) : ((2 * int(abs($X) + 1)) - 1));
	
		for(my $i = $Start; $i > 1; $i -= 2){
			$Factorial *= $i;
		}
	
		if($X > 0){
			$Temp = ($Factorial / (2 ** int($X))) * sqrt($Pi);
		}else {
			$Temp = ((((-1) ** int(abs($X) + 1)) * (2 ** int(abs($X) + 1))) / $Factorial) * sqrt($Pi);
		}
	
		# 半整数 ガンマ関数 Gamma Function
		$GammaFunction = $Temp
	}else {
		if($X > 0){
			$Temp = exp(&LOGGAMMAFUNCTION($X));
		}else {
			$Temp = $Pi / (sin($Pi * $X) * exp(&LOGGAMMAFUNCTION(1 - $X)));
		}
	
		# ガンマ関数(近似) Gamma Function Approximation
		$GammaFunction = $Temp
	}
	
	return $GammaFunction;
}
	
# ログガンマ関数(近似) Log Gamma Function Approximation
# 引数 値 ($x)
# 戻り値 ログガンマ関数(近似) ($LogGammaFunction)
sub LOGGAMMAFUNCTION{
	my ($x) = @_;
	my $LogGammaFunction = 0;
	my @Bernoulli = ((1 / 12), (1 / 360), (1 / 1260), (1 / 1680), (1 / 1188), (691 / 360360), (1 / 156), (3617 / 122400), (43867 / 244188), (174611 / 125400), (77683 / 5796), (236364091 / 1506960));
	my $Pi = atan2(1, 1) * 4;
	my $X = $x;
	my $Sum = 0;
	my $Sign = 0;
	my $Power1 = 1;
	my $Power2 = 1;
	my $Count = @Bernoulli - 1;
	
	for($X = $x; $X <= $Count; $X++){
		$Power1 *= $X;
	}
	$Power2 = 1 / ($X * $X);
	
	for(my $i = $Count; $i >= 0; $i--){
		# 符号
		$Sign = (($i % 2) == 0 ? 1: -1);
	
		$Sum = $Sum + ($Sign * $Bernoulli[$i]);
		$Sum = $Sum * $Power2 if($i != 0);
	}
	$Sum = $Sum / $X;
	
	# ログガンマ関数(近似) Log Gamma Function Approximation
	$LogGammaFunction = ((1 / 2) * log(2 * $Pi)) - log($Power1) - $X + (($X - (1 / 2)) * log($X)) + $Sum;
	
	return $LogGammaFunction;
}