# 第１種変形ベッセル関数(漸化式) Modified Bessel Function of the First Kind
# 引数 次数 変数 ($N, $X)
# 戻り値 第１種変形ベッセル関数 ($ModifiedBesselFunction)
sub MODIFIEDBESSELFUNCTIONFIRST{
	my ($N, $X) = @_;
	my $ModifiedBesselFunction = 0;
	my $NextModifiedBessel = 0;
	my $ModifiedBessel = 0;
	my $PrevModifiedBessel = 0;
	my $Diff = 0;
	my $Count = abs($N);

	if((-2 < $N) && ($N < 2)){
		# 第１種変形ベッセル関数(展開) Modified Bessel Function of the First Kind 
		$ModifiedBessel = &MODIFIEDBESSELFUNCTIONEXPANSION(abs($N), $X);
	}else {
		if(2 <= $N){
			$Diff = $N - int($N);
		}else {
			$Diff = abs($N - int($N));
		}

		# 第１種変形ベッセル関数(展開) Modified Bessel Function of the First Kind 
		$PrevModifiedBessel = &MODIFIEDBESSELFUNCTIONEXPANSION($Diff, $X);
		$ModifiedBessel = &MODIFIEDBESSELFUNCTIONEXPANSION($Diff + 1, $X);

		for(my $i = $Diff + 1; $i < $Count; $i++){
			$NextModifiedBessel = $PrevModifiedBessel - (((2 * $i) / $X) * $ModifiedBessel);

			# 漸化式
			$PrevModifiedBessel = $ModifiedBessel;
			$ModifiedBessel = $NextModifiedBessel;
		}
	}

	# 第１種変形ベッセル関数 Modified Bessel Function of the First Kind
	$ModifiedBesselFunction = $ModifiedBessel;

	return $ModifiedBesselFunction;
}

# 第１種変形ベッセル関数(展開) Modified Bessel Function of the First Kind 
# 引数 次数 変数 ($N, $X)
# 戻り値 第１種変形ベッセル関数(展開) ($ModifiedBesselFunction)
sub MODIFIEDBESSELFUNCTIONEXPANSION{
	my ($N, $X) = @_;
	my $ModifiedBesselFunction = 0;
	my $ModifiedBessel = 0;
	my $PrevModifiedBessel = 0;
	my $Num1 = 0;
	my $Den1 = 0;
	my $Den2 = 0;
	my $AbsN = abs($N);
	my $Limit = 100;
	my $Epsilon = 1.0e-10;
	
	for(my $i = 0; $i < $Limit; $i++){
		# 分子
		$Num1 = ($X / 2) ** ((2 * $i) + $AbsN);
		# 分母
		$Den1 = ($i == 0 ? 1 : ($Den1 * $i));
		$Den2 = &GAMMAFUNCTION($i + $AbsN + 1);

		# 一つ前
		$PrevModifiedBessel = $ModifiedBessel;
		# 第１種変形ベッセル関数 Modified Bessel Function of the First Kind 
		$ModifiedBessel += ($Num1 / ($Den1 * $Den2));

		# 収束判定
		last if(abs($ModifiedBessel - $PrevModifiedBessel) < $Epsilon);
	}

	# 第１種変形ベッセル関数 Modified Bessel Function of the First Kind 
	$ModifiedBesselFunction = $ModifiedBessel;

	return $ModifiedBesselFunction;
}
	
# ガンマ関数(近似) 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;
}