use Math::Cephes qw(gamma);

# 第1種ベッセル関数(漸化式) Bessel Function Of The First Kind
# 引数 次数 変数 ($N, $X)
# 戻り値 第1種ベッセル関数 ($BesselFunction)
sub BESSELFUNCTION{
	my ($N, $X) = @_;
	my $BesselFunction = 0;
	my @PrevBesselFunction = ();
	my $Num = 0;
	my $Diff = 0;
	my $Limit = abs($N);

	if((-2 < $N) && ($N < 2)){
		# 展開式
		# 第1種ベッセル関数 Bessel Function Of The First Kind
		$BesselFunction = &BESSELFUNCTIONEXPANSION($N, $X);
	}else {
		if(2 <= $N){
			$Num = 1;
			$Diff = $N - int($N);
		}else {
			$Num = -1;
			$Diff = $N - int($N);
		}

		# 展開式
		$PrevBesselFunction[1] = &BESSELFUNCTIONEXPANSION($Diff, $X);
		$PrevBesselFunction[0] = &BESSELFUNCTIONEXPANSION($Diff + $Num, $X);

		# 漸化式
		for(my $i = $Diff + $Num; abs($i) < $Limit; $i += $Num){
			my $tmp = (((2 * $i) / $X) * $PrevBesselFunction[0]) - $PrevBesselFunction[1];

			# 第1種ベッセル関数 Bessel Function Of The First Kind
			$BesselFunction = $tmp;

			# [0] 一つ前 [1] 二つ前
			$PrevBesselFunction[1] = $PrevBesselFunction[0];
			$PrevBesselFunction[0] = $tmp;
		}
	}

	return $BesselFunction;
}

# 第1種ベッセル関数(展開) Bessel Function Of The First Kind
# 引数 次数 変数 ($N, $X)
# 戻り値 第1種ベッセル関数 ($BesselFunction)
sub BESSELFUNCTIONEXPANSION{
	my ($N, $X) = @_;
	my $BesselFunction = 0;
	my $PrevBesselFunction = 0;
	my $Num1 = 0;
	my $Den1 = 0;
	my $Den2 = 0;
	my $Limit = 100;
	my $Epsilon = 1.0e-10;
	
	for(my $i = 0; $i < $Limit; $i++){
		# 分子
		$Num1 = ((-1) ** $i);
		# 分母
		$Den1 = ($i == 0 ? 1 : ($Den1 * $i));
		$Den2 = &gamma($N + $i + 1);

		# 一つ前
		$PrevBesselFunction = $BesselFunction;
		# 第1種ベッセル関数(展開) Bessel Function Of The First Kind
		$BesselFunction += ($Num1 / ($Den1 * $Den2)) * (($X / 2) ** ((2 * $i) + $N));

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

	return $BesselFunction;
}