use warnings;
use strict;

# 例 x^2 定積分の近似値
my $a = 0;
my $b = 1;
my $N = 100;
my $Integrand = '$x*$x';
	
my $TrapezoidalRule = &TRAPEZOIDALRULE($a, $b, $N, $Integrand);
print "$TrapezoidalRule\n";

# 台形公式 Trapezoidal Rule
# 引数 左端a 右端b 分割数 被積分関数 ($a, $b, $N, \@Integrand)
# 戻り値 台形公式 ($TrapezoidalRule)
sub TRAPEZOIDALRULE{
	my ($a, $b, $N, $Integrand) = @_;
	   ($a, $b) = ($b, $a) if($a > $b);
	my $TrapezoidalRule = 0;
	my $X = 0;
	my $Sum = 0;
	my $TempA = 0;
	my $TempB = 0;
	my $Pi = atan2(1, 1) * 4;
	my $h = (2 * $Pi) / $N;
	my $HyperbolicSine = 0;
	my $HyperbolicCosine = 0;
	my $Pi_2 = $Pi / 2;
	my $BA_2 = (($b - $a) / 2);
	my $N_2 = int(($N / 2) + 0.5);
	
	# 配列数と分割数の確認
	if($N < 2){
		return 0;
	}
	
	# 計算	
	for(my $i = -$N_2; $i <= $N_2; $i++){
		$HyperbolicSine = &HYPERBOLICSINE($i * $h);
		$HyperbolicCosine = &HYPERBOLICCOSINE($Pi_2 * $HyperbolicSine);
		$X = (($BA_2 * &HYPERBOLICTANGENT($Pi_2 * $HyperbolicSine)) + (($b + $a) / 2));
		$TempA = &INTEGRAND($X, $Integrand);
		$TempB = &HYPERBOLICCOSINE($i * $h) / ($HyperbolicCosine * $HyperbolicCosine);
	
		$Sum += ($TempA * $TempB);
	}
	
	# 台形公式 Trapezoidal Rule
	$TrapezoidalRule = ((($b - $a) * $Pi * $h) / 4) * $Sum;
	
	return $TrapezoidalRule;
}

# 被積分関数 Integrand
# 引数 積分変数 被積分関数 (\@Integrand)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
	my ($x, $Integrand) = @_;
	my $tmp = eval($Integrand);

	return $tmp;
}

# 双曲線正弦 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;
}