# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 正則化された不完全ベータ関数 ($RegularizedIncompleteBetaFunction)
sub REGULARIZEDINCOMPLETEBETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $RegularizedIncompleteBetaFunction = 0;
	
	# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
	$RegularizedIncompleteBetaFunction = &INCOMPLETEBETAFUNCTION($AA, $BB, $X) / &BETAFUNCTION($AA, $BB);
	
	return $RegularizedIncompleteBetaFunction;
}
	
# ベータ関数(近似) Beta Function Approximation
# 引数 値 ($X)
# 戻り値 ベータ関数(近似) ($BetaFunction)
sub BETAFUNCTION{
	my ($X, $Y) = @_;
	my $BetaFunction = 0;
	
	# 値の確認
	if(($X <= 0) && (abs($X - int($X)) == 0) || ($Y <= 0) && (abs($Y - int($Y)) == 0)){
		return 0;
	}
	
	# ベータ関数(近似) Beta Function Approximation
	$BetaFunction = (&GAMMAFUNCTION($X) * &GAMMAFUNCTION($Y)) / &GAMMAFUNCTION($X + $Y);
	
	return $BetaFunction;
}
	
# ガンマ関数(近似) 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($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;
}
	
# 不完全ベータ関数 Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 不完全ベータ関数 ($IncompleteBetaFunction)
sub INCOMPLETEBETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $IncompleteBetaFunction = 0;
	my $TrapezoidalRule = 0;
	my $T = 0;
	my $Sum = 0;
	my $HyperbolicSine = 0;
	my $HyperbolicCosine = 0;
	my $TempA = 0;
	my $TempB = 0;
	my $a = 0;
	my $b = $X;
	my $N = 100;
	my $Pi = atan2(1, 1) * 4;
	my $Pi_2 = $Pi / 2;
	my $BA_2 = (($b - $a) / 2);
	my $h = (2 * $Pi) / $N;
	my $Begin = -($N / 2);
	my $End = $N / 2;
	
	# 値の確認
	if(($AA <= 0) || ($BB <= 0) || ($X < 0) || (1 < $X)){
		return 0;
	}
	
	# 計算
	for(my $i = $Begin; $i <= $End; $i++){
		$HyperbolicSine = &HYPERBOLICSINE($i * $h);
		$HyperbolicCosine = &HYPERBOLICCOSINE($Pi_2 * $HyperbolicSine);
		$T = (($BA_2 * &HYPERBOLICTANGENT($Pi_2 * $HyperbolicSine)) + (($b + $a) / 2));
		$TempA = &FUNCTION($AA, $BB, $T);
		$TempB = &HYPERBOLICCOSINE($i * $h) / ($HyperbolicCosine * $HyperbolicCosine);
	
		$Sum += ($TempA * $TempB);
	}
	
	# 台形公式 Trapezoidal Rule
	$TrapezoidalRule = ((($b - $a) * $Pi * $h) / 4) * $Sum;
	
	# 不完全ベータ関数 Incomplete Beta Function
	$IncompleteBetaFunction = $TrapezoidalRule;
	
	return $IncompleteBetaFunction;
}
	
# 関数 Function
# 引数 積分変数 被積分関数 (\@Integrand)
# 戻り値 関数 ($Function)
sub FUNCTION{
	my ($AA, $BB, $X) = @_;
	my $Function = ($X ** ($AA - 1)) * ((1 - $X) ** ($BB - 1));
	
	return $Function;
}
	
# 双曲線正弦 Hyperbolic Sine
# 引数 値 ($X)
# 戻り値 双曲線正弦 ($HyperbolicSine)
sub HYPERBOLICSINE{
	my ($X) = @_;
	my $HyperbolicSine = 0;
	
	# 双曲線正弦 Hyperbolic Sine
	$HyperbolicSine = (exp($X) - exp(-$X)) / 2;
	
	return $HyperbolicSine;
}
	
# 双曲線余弦 Hyperbolic Cosine
# 引数 値 ($X)
# 戻り値 双曲線余弦 ($HyperbolicCosine)
sub HYPERBOLICCOSINE{
	my ($X) = @_;
	my $HyperbolicCosine = 0;
	
	# 双曲線余弦 Hyperbolic Cosine
	$HyperbolicCosine = (exp($X) + exp(-$X)) / 2;
	
	return $HyperbolicCosine;
}
	
# 双曲線正接 Hyperbolic Tangent
# 引数 値 ($X)
# 戻り値 双曲線正接 ($HyperbolicTangent)
sub HYPERBOLICTANGENT{
	my ($X) = @_;
	my $HyperbolicSine = (exp($X) - exp(-$X)) / 2;
	my $HyperbolicCosine = (exp($X) + exp(-$X)) / 2;
	my $HyperbolicTangent = 0;
	
	# 双曲線正接 Hyperbolic Tangent
	$HyperbolicTangent = $HyperbolicSine / $HyperbolicCosine;
	
	return $HyperbolicTangent;
}