use warnings;
use strict;

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

# ガウス-ルジャンドル数値積分 Legendre-Gauss Quadrature
# 引数 左端a 右端b 分割数 被積分関数 ($a, $b, $N, \@Integrand)
# 戻り値 ガウス-ルジャンドル数値積分 ($LegendreGaussQuadrature)
sub LEGENDREGAUSSQUADRATURE{
	my ($a, $b, $N, $Integrand) = @_;
	   ($a, $b) = ($b, $a) if($a > $b);
	my $LegendreGaussQuadrature = 0;
	my @NewtonsMethod = ();
	my $LegendrePolynomial = 0;
	my $X = 0;
	my $Sum = 0;
	my $AB2 = ($a + $b) / 2;
	my $BA2 = ($b - $a) / 2;

	# 分割数の確認
	if($N < 2){
		return 0;
	}
	
	# 計算
	for(my $i = 1; $i <= $N; $i++){
		# ニュートン法 Newton's Method
		@NewtonsMethod = &NEWTONSMETHOD($N, $i);
		$X = $NewtonsMethod[0];
		$LegendrePolynomial = $NewtonsMethod[1];
		# 重み
		my $Weight = (2 * (1 - ($X * $X))) / (($N * $LegendrePolynomial) * ($N * $LegendrePolynomial));

		$X = $AB2 + ($BA2 * $X);
		my $tmp = &INTEGRAND($X, $Integrand);
	
		$Sum += $Weight * $tmp;
	}

	# ガウス-ルジャンドル数値積分 Legendre-Gauss Quadrature
	$LegendreGaussQuadrature = $BA2 * $Sum;
	
	return $LegendreGaussQuadrature;
}

# ニュートン法 Newton's Method
# 引数 次数 位置 ($N, $K)
# 戻り値 ニュートン法 ($NewtonsMethod)
sub NEWTONSMETHOD{
	my ($N, $K) = @_;
	my @NewtonsMethod = ();
	my @LegendrePolynomial = ();
	my $Pi = atan2(1, 1) * 4;
	my $X = cos(($Pi * ($K - 0.25)) / ($N + 0.5));
	my $PrevX = 0;
	my $Limit = 100;
	my $Epsilon = 1.0e-20;

	# 計算
	for(my $i = 0; $i < $Limit; $i++){
		# ルジャンドル多項式 Legendre Polynomial
		@LegendrePolynomial = &LEGENDREPOLYNOMIAL($N, $X);
		my $Num = $LegendrePolynomial[0];
		my $Den = (($N * ($LegendrePolynomial[1] - ($X * $LegendrePolynomial[0])))  / (1 - ($X * $X)));

		# 近似解
		$PrevX = $X;
		$X = $X - ($Num / $Den);
		# ニュートン法 Newton's Method
		# [0] 近似解 [1] 一つ前のルジャンドル多項式の値
		$NewtonsMethod[0] = $X;
		$NewtonsMethod[1] = $LegendrePolynomial[1];
	
		# 収束判定
		last if(abs($X - $PrevX) < $Epsilon);
	}
	
	return @NewtonsMethod;
}

# ルジャンドル多項式 Legendre Polynomial
# 引数 次数 値 ($N, $x)
# 戻り値 ルジャンドル多項式 (@LegendrePolynomial)
sub LEGENDREPOLYNOMIAL{
	my ($N, $x) = @_;
	my @LegendrePolynomial = ($x , 1);

	# 計算
	for(my $i = 1; $i < $N; $i++){
		my $tmp = ((((2 * $i) + 1) * $x * $LegendrePolynomial[0]) - ($i * $LegendrePolynomial[1])) / ($i + 1);

		# [0] 最新 [1] 一つ前
		# ルジャンドル多項式 Legendre Polynomial
		$LegendrePolynomial[1] = $LegendrePolynomial[0];
		$LegendrePolynomial[0] = $tmp;
	}

	return @LegendrePolynomial;
}

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

	return $tmp;
}