# 不完全ベータ関数 Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 不完全ベータ関数 (@IncompleteBetaFunction)
sub INCOMPLETEBETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $IncompleteBetaFunction = 0;
	my $SimpsonsRule = 0;
	my $a = 0;
	my $b = $X;
	my $N = 100;
	my $h = ($b - $a) / $N;
	my $Sum = 0;
	my $SumA = 0;
	my $SumB = 0;
	
	# $AA,$BBが一未満だと機能しない?
	if(($X < 0) || (1 < $X) || ($AA <= 0) || ($BB <= 0)){
		return 0;
	}

	# シンプソンの公式 Simpson's Rule
	$SumA = ($a ** ($AA - 1)) * ((1 - $a) ** ($BB - 1));
	$SumB = ($b ** ($AA - 1)) * ((1 - $b) ** ($BB - 1));
	
	# 分割数は適当
	for(my $i = 1; $i < $N; $i++){
		my $t = $a + ($i * $h);
		my $num = ($i % 2 == 0 ? 2 : 4);
		my $tmp = ($t ** ($AA - 1)) * ((1 - $t) ** ($BB - 1));
	
		$Sum += $num * $tmp;
	}
	
	# 不完全ベータ関数 Incomplete Beta Function
	$IncompleteBetaFunction = (1 / &BETAFUNCTION($AA, $BB)) * (($SumA + $Sum + $SumB) * ($h / 3));
	
	return $IncompleteBetaFunction;
}

# ベータ関数 Beta Function
# 引数 変数 変数 ($AA, $BB)
# 戻り値 ベータ関数 (@BetaFunction)
sub BETAFUNCTION{
	my ($AA, $BB) = @_;
	my $BetaFunction = 0;

	# 一未満かつ0.5以外だと値が一致しない
	if((($AA < 1) && ($AA != 0.5)) || (($BB < 1) && ($BB != 0.5))){
		return 0;
	}

	# ベータ関数 Beta Function
	$BetaFunction = (&GAMMAFUNCTION($AA) * &GAMMAFUNCTION($BB)) / &GAMMAFUNCTION($AA + $BB);

	return $BetaFunction;
}

# ガンマ関数 Gamma Function
# 引数 変数 ($BB)
# 戻り値 ガンマ関数 (@GammaFunction)
sub GAMMAFUNCTION{
	my ($BB) = @_;
	my $GammaFunction = 0;
	my $X = $BB;
	my $Factorial = 1;
	my $Diff = $X - int($X);
	my $Pi = atan2(1, 1) * 4;

	if($Diff == 0){
		for(my $i = $X - 1; $i > 1; $i--){
			$Factorial *= $i;
		}

		# N-1の階乗
		$GammaFunction = $Factorial;
	}elsif($Diff == 0.5){
		for(my $i = int($X); $i > 0; $i--){
			$Factorial *= ((2 * $i) + 1) / 2;
		}

		# 半整数
		$GammaFunction = $Factorial * sqrt($Pi);
	}else {
		my $X1 = 1 / (12 * $X);
		my $X2 = 1 / (288 * ($X * $X));
		my $X3 = 139 / (51840 * ($X * $X * $X));
		my $X4 = 571 / (2488320 * ($X * $X * $X * $X));
		my $X5 = 163879 / (209018880 * ($X * $X * $X * $X * $X));
		my $X6 = 5246819 / (75246796800  * ($X * $X * $X * $X * $X * $X));
		my $X7 = 534703531 / (902961561600  * ($X * $X * $X * $X * $X * $X * $X));
		my $X8 = 4483131259 / (86684309913600 * ($X * $X * $X * $X * $X * $X * $X * $X));
		my $X9 = 432261921612371 / (514904800886784000 * ($X * $X * $X * $X * $X * $X * $X * $X));

		# スターリングの近似 Stirling's approximation
		$GammaFunction	 = sqrt(2 * $Pi) * ($X ** ($X - (1 / 2))) * exp(-$X) * (1 + $X1 + $X2 - $X3 - $X4 + $X5 + $X6 - $X7 - $X8 + $X9);
	}

	return $GammaFunction;
}