package Algorism;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
use Math::Matrix;
use Math::MatrixReal;

use Sci qw($pi tanh cosh sinh);

my $pi          = 3.1415926535;

#================================================
# 静的メンバー関数
#================================================
sub InterpolateValuesInCube
{
	my ($nx, $ny, $nz, $pValues, $x, $y, $z) = @_;

	$x = Utils::Reduce01($x);
	$y = Utils::Reduce01($y);
	$z = Utils::Reduce01($z);
	my $xx = $x * $nx;
	my $yy = $y * $ny;
	my $zz = $z * $nz;
	my $ix = int( $xx + 0.0 );
	my $iy = int( $yy + 0.0 );
	my $iz = int( $zz + 0.0 );
	my ($ix1, $iy1, $iz1) =  ($ix+1, $iy+1, $iz+1);
$ix1 -= $nx if($ix1 >= $nx);
$iy1 -= $ny if($iy1 >= $ny);
$iz1 -= $nz if($iz1 >= $nz);

	my $a0   = $pValues->[$ix][$iy][$iz];
	my $ax   = $pValues->[$ix1][$iy][$iz]   - $a0;
	my $ay   = $pValues->[$ix][$iy1][$iz]   - $a0;
	my $az   = $pValues->[$ix][$iy][$iz1]   - $a0;
	my $axy  = $pValues->[$ix1][$iy1][$iz]  - $a0 - $ax - $ay;
	my $ayz  = $pValues->[$ix][$iy1][$iz1]  - $a0 - $ay - $az;
	my $azx  = $pValues->[$ix1][$iy][$iz1]  - $a0 - $az - $ax;
	my $axyz = $pValues->[$ix1][$iy1][$iz1] - $a0 - $ax - $ay - $az - $axy - $ayz - $azx;

	return $a0 + $ax * ($xx-$ix) + $ay * ($yy-$iy) + $az * ($zz-$iz)
	     + $axy * ($xx-$ix)*($yy-$iy) + $ayz * ($yy-$iy)*($zz-$iz) + $azx * ($zz-$iz)*($xx-$ix)
	     + $axyz * ($xx-$ix)*($yy-$iy)*($zz-$iz);
#	return $pValues->[$ix][$iy][$iz];
}

sub Integrate3DSimple
{
	my ($nx, $ny, $nz, $pValues) = @_;

	my $S = 0.0;
	for(my $ix = 0 ; $ix < $nx ; $ix++) {
	for(my $iy = 0 ; $iy < $ny ; $iy++) {
	for(my $iz = 0 ; $iz < $nz ; $iz++) {
		$S += $pValues->[$ix][$iy][$iz];
	}
	}
	}

	return $S;
}

#周期境界条件を考えるときは、単純和に一致する
sub Integrate3D
{
	my ($nx, $ny, $nz, $pValues) = @_;

	my $S = 0.0;
#	for(my $ix = 0 ; $ix < $nx ; $ix++) {
	for(my $ix = 0 ; $ix < $nx-1 ; $ix++) {
#	for(my $iy = 0 ; $iy < $ny ; $iy++) {
	for(my $iy = 0 ; $iy < $ny-1 ; $iy++) {
#	for(my $iz = 0 ; $iz < $nz ; $iz++) {
	for(my $iz = 0 ; $iz < $nz-1 ; $iz++) {
		my ($ix1, $iy1, $iz1) =  ($ix+1, $iy+1, $iz+1);
#		$ix1 = 0 if($ix1 == $nx);
#		$iy1 = 0 if($iy1 == $ny);
#		$iz1 = 0 if($iz1 == $nz);

		my $a0   = $pValues->[$ix][$iy][$iz];
		my $ax   = $pValues->[$ix1][$iy][$iz]   - $a0;
		my $ay   = $pValues->[$ix][$iy1][$iz]   - $a0;
		my $az   = $pValues->[$ix][$iy][$iz1]   - $a0;
		my $axy  = $pValues->[$ix1][$iy1][$iz]  - $a0 - $ax - $ay;
		my $ayz  = $pValues->[$ix][$iy1][$iz1]  - $a0 - $ay - $az;
		my $azx  = $pValues->[$ix1][$iy][$iz1]  - $a0 - $az - $ax;
		my $axyz = $pValues->[$ix1][$iy1][$iz1] - $a0 - $ax - $ay - $az - $axy - $ayz - $azx;
#print "a=($a0, $ax, $ay, $az)\n";
		$S += $a0 + 0.5 * ($ax + $ay + $az) + 0.25 * ($axy + $ayz + $azx) + 0.125 * $axyz;
#print "S=$S\n";
#exit;
	}
	}
	}

	return $S;
}

sub ConvolutionByArray
{
	my ($pX, $pY, $pfArray, $w) = @_;

	my $nData = @$pX;
	my $ix0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		if($pX->[$i] >= 0.0) {
			$ix0 = $i;
			last;
		}
	}

	my $x0 = $pX->[0];
	my $x1 = $pX->[$nData - 1];
	my $dx = $pX->[1] - $pX->[0];
	my $dimax = (defined $w)? int($w / $dx + 0.99999) : undef;
#print "ix0=$ix0, dimax=$dimax, w=$w, dx=$dx\n";
#exit;
	my @y;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x0 = $pX->[$i];
		my $y0 = $pY->[$i];
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $di = $i - $j;
			next if($w > 0 and abs($di) > $dimax);
			my $ia = $ix0 + $di;
			next if($ia < 0 and $nData <= $ia);
			my $f = $pfArray->[$ia];
			$y[$j] += $y0 * $f * $dx;
		}
	}

	return \@y;
}

sub ConvolutionByFunc
{
	my ($pX, $pY, $pfunc, $w) = @_;

	my $nData = @$pX;
	my $x0 = $pX->[0];
	my $x1 = $pX->[$nData - 1];
	my $dx = $pX->[1] - $pX->[0];
	my @y;
#print "dx0,f\n";
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x0 = $pX->[$i];
		my $y0 = $pY->[$i];
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $dx0 = ($i-$j) * $dx;
			next if($w > 0.0 and abs($dx0) > $w);
			my $f = &$pfunc($dx0);
#print "$dx0,$f\n";
			$y[$j] += $y0 * $f * $dx;
		}
#exit;
	}

	return \@y;
}

sub Convolution
{
	my ($pX, $pY, $Width, $Function, $IgnoreZeroValue, $Normalize) = @_;
	$Function = "Gaussian" if(!defined $Function);
	$IgnoreZeroValue = 0 if(!defined $IgnoreZeroValue);
#	my $pi = Sci::pi();

	my $nData = @$pX;
	my $x0 = $pX->[0];
	my $x1 = $pX->[$nData - 1];
	my $dx = $pX->[1] - $pX->[0];
	my $di = int(($Width * 5.0) / $dx + 1.1);
	my $coeff = 1/sqrt($pi)/($Width/$dx);
	my @y;
	for(my $j = 0 ; $j < $nData ; $j++) {
		my $x0 = $pX->[$j];
		my $y0 = $pY->[$j];
		for(my $k = -$di ; $k <= $di ; $k++) {
			next if($j+$k < 0 or $j+$k >= $nData);
			next if($IgnoreZeroValue and $pY->[$j+$k] == 0.0);

			my $dvx = $dx*$k;
			$y[$j+$k] += $coeff * $y0 * exp(-$dvx*$dvx/$Width/$Width);
		}
	}
	if($Normalize) {
		my @S;
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $x0 = $pX->[$j];
			my $y0 = $pY->[$j];
			$S[$j] = 0.0;
			for(my $k = -$di ; $k <= $di ; $k++) {
				next if($j+$k < 0 or $j+$k >= $nData);
				next if($IgnoreZeroValue and $pY->[$j+$k] == 0.0);

				my $dvx = $dx*$k;
				$S[$j] += $coeff * exp(-$dvx*$dvx/$Width/$Width);
			}
		}
		for(my $j = 0 ; $j < $nData ; $j++) {
			next if($S[$j] == 0.0);
			$y[$j] /= $S[$j];
		}
	}

	return \@y;
}

sub minver
{
	my ($A) = @_;

#	my $A = new Math::MatrixReal(2, 2);
#	$A->assign(1, 1, $this->{"a11"});
#	$A->assign(1, 2, $this->{"a21"});
#	$A->assign(2, 1, $this->{"a12"});
#	$A->assign(2, 2, $this->{"a22"});
#	print "A:<br>\n";print $A;print "<br>\n";
	return $A->inverse();
#	print "RA:<br>\n";print $RA;print "<br>\n";
#	my $X = new Math::MatrixReal(2,1);
#	$X->assign(1, 1, $xc);
#	$X->assign(2, 1, $yc);
#	print "X:<br>\n";print $X;print "<br>\n";
#	my $XX = $RA * $X;
#	print "XX:<br>\n";print $XX;print "<br>\n";
}

#入力
# $pX->[$i], $pY->[$i], $n : データ
# $m : 近似多項式の次数
#出力
# \@c : 多項式のi次の係数へのリファレンス
sub lstsq
{
	my ($pX, $pY, $n, $m, $c) = @_;
	return 999 if($m >= $n or $m < 1);

	my $nMtx = $m + 1;
	my $W = new Math::MatrixReal($nMtx, $nMtx);
	my $V = new Math::MatrixReal($nMtx, 1);
	
	for(my $j = 0 ; $j < $m+1 ; $j++) {
		for(my $k = $j ; $k < $m+1 ; $k++) {
			my $w = 0.0;
			for(my $i = 0 ; $i < $n ; $i++) {
				my $x = $pX->[$i];
				my $y = $pY->[$i];
				$w += $x ** ($j+$k);
			}
			$W->assign($j+1, $k+1, $w);
			$W->assign($k+1, $j+1, $w);
		}
		my $v = 0.0;
		for(my $i = 0 ; $i < $n ; $i++) {
			my $x = $pX->[$i];
			my $y = $pY->[$i];
			$v += $y * $x ** $j;
		}
		$V->assign($j+1, 1, $v);
	}

#print "W:\n";print $W;print "\n";
#print "V:\n";print $V;print "\n";
	my $RW = $W->inverse();
#print "RW:\n";print $RW;print "\n";
	my $XX = $RW * $V;
#print "XX:\n";print $XX;print "\n";

	for(my $i = 0 ; $i < $m+1 ; $i++)
	{
		$c->[$i] = $XX->element($i+1, 1);
	}
	return 0;
}

sub SolveByNewtonRapson
{
	my ($pFunc, $pDiffFunc, $iniX, $dX, $ramda, $Xth, $fth, $nMaxIter, $IsPrint) = @_;
	
	my $x = $iniX;
	my ($Xdiff, $fdiff, $nIter, $mess);

	my $delta;
 	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
 		my $f1;
 		if($pDiffFunc) {
			$f1 = &$pDiffFunc($x);
#print "f1=$f1\n";
 		}
 		else {
			my $f  = &$pFunc($x-$dX);
			my $f2 = &$pFunc($x+$dX);
#print "dX=$dX: $f2 - $f = $df\n";
			if(defined $f and defined $f2) {
				my $df = $f2 - $f;
				$f1 = $df / 2.0 / $dX;
			}
 		}
		if(!defined $f1) {
			return (undef, undef, undef, undef, "function did not return a value");
		}

		my $f  = &$pFunc($x);
		if($f1 < 0.0) {
 			$f1 -= $ramda;
 		}
		else {
 			$f1 += $ramda;
		}
		$delta = -$f / $f1;
print "$i: x = $x => " if($IsPrint);
		$x += $delta;
print "$x ($delta)\n" if($IsPrint);

		if($Xth > 0 and abs($delta) < $Xth) {
			$nIter = $i+1;
			$fdiff = &$pFunc($x);
			$mess = "converged for X";
			return ($x, $delta, $fdiff, $nIter, $mess);
		}
		if($fth > 0 and abs($f) < $fth) {
			$nIter = $i+1;
			$fdiff = &$pFunc($x);
			$mess = "converged for f";
			return ($x, $delta, $fdiff, $nIter, $mess);
		}
	}

	$mess  = "Error: not converged";
	return (undef, $delta, undef, $nMaxIter, $mess);
}

sub SolveBySecant
{
	my ($pFunc, $iniX, $initialdX, $ramda, $Xth, $fth, $nMaxIter, $IsPrint) = @_;
	
	my $x = $iniX;
	my ($Xdiff, $fdiff, $nIter, $mess);

	my $fprev = &$pFunc($x);
	my $dX    = $initialdX;
	$x += $dX;
 	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
		my $f2 = &$pFunc($x);
#print "dX=$dX: $f2 - $f = $df\n";
		if(!defined $f2) {
			return (undef, undef, undef, undef, "function did not return a value");
		}

		my $df = $f2 - $fprev;
		my $f1 = $df / $dX;
		if($f1 < 0.0) {
 			$f1 -= $ramda;
 		}
		else {
 			$f1 += $ramda;
		}
		$dX = -$f2 / $f1;
print "$i: x = $x => " if($IsPrint);
		$x += $dX;
print "$x ($dX)\n" if($IsPrint);

		if($Xth > 0 and abs($dX) < $Xth) {
			$nIter = $i+1;
			$fdiff = &$pFunc($x);
			$mess = "converged for X";
			return ($x, $dX, $fdiff, $nIter, $mess);
		}
		if($fth > 0 and abs($f2) < $fth) {
			$nIter = $i+1;
			$fdiff = &$pFunc($x);
			$mess = "converged for f";
			return ($x, $dX, $fdiff, $nIter, $mess);
		}

		$fprev = $f2;
	}
	
	$mess  = "Error: not converged";
	return (undef, $dX, undef, $nMaxIter, $mess);
}

sub SolveByBisection
{
	my ($pFunc, $x0, $x1, $Xth, $fth, $nMaxIter, $IsPrint) = @_;

	my $mess = 'converged';
	my $fmin = &$pFunc($x0);
	my $fmax = &$pFunc($x1);
	if($fmin * $fmax > 0.0) {
#print "f(x0)=$fmin and f(x1)=$fmax have same sign.\n" if($IsPrint);
		$mess = "Error: f(x0) and f(x1) have the same sign\n";
		return (undef, undef, undef, undef, $mess);
	}

	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
		my $dx = abs($x1 - $x0);
#print "$i: $dx\n" if($IsPrint);
		return (($x0+$x1)/2.0, $dx, ($fmin+$fmax)/2.0, $i, $mess) if(abs($dx) < $Xth);

		my $xhalf  = ($x0 + $x1) / 2.0;
		my $fhalf = &$pFunc($xhalf);
#print "$i: E: $x0 - $xhalf - $x1: $fmin - $fhalf - $fmax\n" if($IsPrint);
print "$i: $x0 - $xhalf - $x1: f=$fhalf\n" if($IsPrint);
		return (($x0+$x1)/2.0, $dx, $fhalf, $i, $mess) if(abs($fhalf) < $fth);

		if($fmin * $fhalf < 0.0) {
			$x1   = $xhalf;
			$fmax = $fhalf;
		}
		else {
			$x0   = $xhalf;
			$fmin = $fhalf;
		}
	}

	$mess  = "Error: not converged";
	return (undef, undef, undef, $nMaxIter, $mess);
}

sub SolveByNewton
{
	my ($func, $iniX, $dX, $Xth, $fth, $nMaxIter) = @_;
	
	my $x = $iniX;
	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
		my $curx = $x;
		my $f = eval($func);
#print "$i: x=$x: y=$f\n";
		$x += $dX;
		my $f2 = eval($func);
		my $df = $f2-$f;
		my $f1 = $df / $dX;

		my $delta = -$f / $f1;
#print "  $f-$f2=$f1: dX=$delta\n";
#print "  x:$curx  d:$delta\n";
		$x = $curx + $delta;
#print "  x:$x (new)\n";
#print "  dx:$delta ($Xth)  df:$df ($fth)\n";
		last if($Xth > 0 and abs($delta) < $Xth);
		last if($fth > 0 and abs($df) < $fth);
	}
	return $x;
}

sub IsFactor
{
	my ($I, $idiv) = @_;
	return 0 if($idiv <= 0);
	my $res = $I - int($I / $idiv) * $idiv;
	return 1 if(abs($res) < 1.0e-10);
	return 0;
}

sub Factorize
{
	my ($I) = @_;
	my $iMax = $I / 2;
	my @check;
	for(my $i = 0 ; $i <= $iMax ; $i++) {
		$check[$i] = 0;
	}
	my @e;
	my $ic = 2;
	while(1) {
		last if($ic > $iMax);
		if($check[$ic]) {
			$ic++;
			next;
		}
		if(IsFactor($I, $ic)) {
			push(@e, $ic);
			$I = $I / $ic;
			next;
		}
		for(my $j = 1 ; ; $j++) {
			last if($ic*$j > $I/2);
			$check[$ic*$j] = 1;
		}
		$ic++;
	}

	push(@e, $I) if(@e == 0);
	return @e;
}

sub SmoothingBySimpleAverage
{
	my ($pY, $n) = @_;

	my $n2 = int($n / 2);
	my $nData = @$pY;
	my @y;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $c = 0;
		$y[$i] = 0.0;
		for(my $k = $i - $n2 ; $k <= $i + $n2 ; $k++) {
			next if($k < 0 or $k >= $nData);
			$y[$i] += $pY->[$k];
			$c++;
		}
		if($c > 0) {
			$y[$i] /= $c;
		}
		else {
			$y[$i] = $pY->[$i];
		}
	}
	return \@y;
}


sub AdaptiveSmoothing
{
	my ($pY, $n) = @_;

	my $m = int($n / 2);
	my $W23 = (4.0 * $m * $m - 1.0) * (2.0 * $m + 3.0) / 3.0;
	my @w23j;
	for(my $j = -$m ; $j <= $m ; $j++) {
		$w23j[$j + $m] = (3.0 * $m * ($m+1.0) - 1.0 - 5.0 * $j * $j) / $W23;
	}

	my $nData = @$pY;
	my @y;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $c = 0.0;
		$y[$i] = 0.0;
		for(my $j = -$m ; $j <= $m ; $j++) {
			my $k = $i + $j;
			next if($k < 0 or $k >= $nData);
			$y[$i] += $w23j[$j+$m] * $pY->[$k];
			$c += $w23j[$j+$m];
		}
		if($c > 0) {
			$y[$i] /= $c;
		}
		else {
			$y[$i] = $pY->[$i];
		}
	}
	return \@y;
}

sub SmoothingBymLSQ
{
	my ($pX, $pY, $nSmooth, $nOrder, $iPrintLevel) = @_;

	my  (@wx, @wy, @x, @y);
	my $nData = @$pX;
	my $nhalf = int($nSmooth/2);

	my $c = 0;
	for(my $i = $nhalf ; $i < $nData-$nhalf ; $i++) {
		my $Opt = new Optimize;
		@wx = ();
		@wy = ();
		for(my $j = 0 ; $j < $nSmooth ; $j++) {
			$wx[$j] = $pX->[$i-$nhalf + $j];
			$wy[$j] = $pY->[$i-$nhalf + $j];
		}
		my ($pCi, $S2) = $Opt->Optimize("mlsq", \@wx, \@wy, $nOrder, $iPrintLevel);
#print "Optimized at (", join(',',@{$pCi}), ") with S2=$S2\n";

		$x[$c] = $pX->[$i];
		$y[$c] = $Opt->YCal($x[$c]);
		$c++;
#print "$x\t$y[$i]\t$ycal\n";
	}

	return (\@x, \@y);
}

sub SmoothingByFFT
{
	my ($pY, $imin, $imax, $IsPrint) = @_;

	my $fft = new MyFFT;

	my ($n, $pX, $pData) = $fft->SetData(undef, $pY, undef, "forward:real", "truncate");
	my ($pFreq, $pReal, $pImag) = $fft->FFT();
#print "n(forward)=$n\n";

	my (@real, @image);
	for(my $i = 0 ; $i < @$pFreq ; $i++) {
		my $x = $pFreq->[$i];
		if($i < $imin or $imax < $i) {
#		if($x < $k0 or $k1 < $x) {
			$real[$i]  = 0.0;
			$image[$i] = 0.0;
		}
		else {
			$real[$i]  = $pReal->[$i];
			$image[$i] = $pImag->[$i];
		}
	}

	($n, $pX, $pData) = $fft->SetData(undef, \@real, \@image, "inverse:real", "truncate");
	my ($pTime, $pOrgReal) = $fft->FFT();
#print "n(inverse)=$n\n";

	return ($pFreq, $pReal, $pImag, $pTime, $pOrgReal);
}

sub IntegrateByRectangle
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	return undef if($nData < 2);

	my @y;
	$y[0] = 0.0; #0.5 * ($pX->[1] - $pX->[0]) * $pY->[0];
	for(my $i = 1 ; $i < $nData ; $i++) {
		$y[$i] = $y[$i-1] + ($pX->[$i] - $pX->[$i-1]) * $pY->[$i-1];
	}
	$y[$nData-1] -= 0.5 * ($pX->[$nData-1] - $pX->[$nData-2]) * $pY->[$nData-1];
	return \@y;
}

sub IntegrateByRectangleWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 1);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-1 ; $i++) {
		my $x = $x0 + $i * $h;
		my $y = &$pFunc($x);
#print "y'[$i]($x)=$y\n";
		$S += $y;
	}
	return $S * $h;
}

sub IntegrateByTrapezoid
{
	my ($pX, $pY, $idx0, $idx1) = @_;

	my $nData = @$pX;
my $nYData = @$pY;
	return undef if($nData < 2);

	$idx0 = 0      if(!defined $idx0);
	$idx1 = $nData if(!defined $idx1);

	my @y;
	for(my $i = 0 ; $i < $idx0 ; $i++) {
		$y[$i] = 0.0;
	}
	
	$y[$idx0]   = 0.0; #0.5 * ($pX->[$idx0+1] - $pX->[$idx0]) * $pY->[$idx0];
	$y[$idx0+1] = 0.5 * ($pX->[$idx0+1] - $pX->[$idx0]) * ($pY->[$idx0+1] + $pY->[$idx0]);
	for(my $i = $idx0+2 ; $i < $idx1 ; $i++) {
#print "i=$idx0, $idx1, $i (n=$nData, $nYData,  p=$pX, $pY)\n";
#print "eq: [y[i] = $y[$i-1] + 0.5 * ($pX->[$i] - $pX->[$i-1]) * ($pY->[$i] + $pY->[$i-1]);]\n";
		$y[$i] = $y[$i-1] + 0.5 * ($pX->[$i] - $pX->[$i-1]) * ($pY->[$i] + $pY->[$i-1]);
#exit if(!defined $y[$i]);
	}

	my $last = $y[$idx1-1];
	for(my $i = $idx1 ; $i < $nData ; $i++) {
		$y[$i] = $last;
	}

	return \@y;
}

sub IntegrateByTrapezoidWithFunction_NotUsed
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 1);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-1 ; $i++) {
		my $xi  = $x0 + $i * $h;
		my $xi1 = $xi + $h;
		$S += &$pFunc($xi) + &$pFunc($xi1);
	}

	return $S / 2.0 * $h;
}

sub IntegrateByTrapezoidWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 1);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.5 * (&$pFunc($x0) + &$pFunc($x1));
	for(my $i = 1 ; $i <= $nBlock-1 ; $i++) {
		my $x  = $x0 + $i * $h;
		my $y = &$pFunc($x);
		$S += $y;
	}

	return $S * $h;
}

sub IntegrateBySimpsonWithFunction_NotUsed
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 2 or $nBlock %2 != 0);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-2 ; $i += 2) {
		my $xi  = $x0  + $i * $h;
		my $xi1 = $xi  + $h;
		my $xi2 = $xi1 + $h;
		$S += &$pFunc($xi) + 4.0 * &$pFunc($xi1) + &$pFunc($xi2);;
	}

	return $S / 3.0 * $h;
}

sub IntegrateBySimpsonWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 2 or $nBlock %2 != 0);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = &$pFunc($x0) + &$pFunc($x1);
	my $S4 = 0.0;
	for(my $i = 1 ; $i <= $nBlock-1 ; $i += 2) {
		my $x  = $x0 + $i * $h;
		my $y = &$pFunc($x);
		$S4 += $y;
	}
	my $S2 = 0.0;
	for(my $i = 2 ; $i <= $nBlock-2 ; $i += 2) {
		my $x  = $x0 + $i * $h;
		my $y = &$pFunc($x);
		$S2 += $y;
	}

	return ($S + 4.0 * $S4 + 2.0 * $S2) / 3.0 * $h;
}

sub SinglePointIntegrateByConstantStepSimpson
{
	my ($nData, $h, $pY) = @_;

	if($nData < 3 or $nData % 2 == 0 or $h <= 0.0) {
		return undef;
	}

	my $w = 0.0;
	for(my $i = 1 ; $i < $nData-1 ; $i += 2) {
		$w += $pY->[$i-1] + 4.0 * $pY->[$i] + $pY->[$i+1];
	}
	my $S = $w * $h / 3.0;
	return $S;
}

sub SinglePointIntegrateBySimpson
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	if($nData < 3 or $nData % 2 == 0) {
		return undef;
	}

	my $w = 0.0;
	for(my $i = 1 ; $i < $nData-1 ; $i += 2) {
		my $h   = 0.5 * ($pX->[$i+1] - $pX->[$i-1]);
		my $hpd = 1.0 / ($pX->[$i]   - $pX->[$i-1]);
		my $hmd = 1.0 / ($pX->[$i+1] - $pX->[$i]);
		my $d   = $pX->[$i] - $pX->[$i-1] - $h;
		$w += ( (1.0 + 2.0*$d*$hpd) * $pY->[$i-1]
			  +  2.0*$h*($hpd+$hmd) * $pY->[$i]
			  + (1.0 - 2.0*$d*$hmd) * $pY->[$i+1] ) * $h;
	}
	my $S = $w / 3.0;
	return $S;
}

sub IntegrateByConstantStepSimpson
{
	my ($nData, $h, $pY) = @_;
	return undef if($nData < 2);

	my @y;
	$y[0] = 0.0;
	$y[1] = $h * $pY->[0];
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my @yy = ($pY->[$i-1], $pY->[$i], $pY->[$i+1]);
		my $S = SinglePointIntegrateByConstantStepSimpson(3, $h, \@yy);
		$y[$i+1] = $y[$i-1] + $S;
	}
	return \@y;
}

sub IntegrateBySimpson
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	return undef if($nData < 2);

	my @y;
	$y[0] = 0.0;
	$y[1] = ($pX->[1] - $pX->[0]) * $pY->[0];
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my @xx = ($pX->[$i-1], $pX->[$i], $pX->[$i+1]);
		my @yy = ($pY->[$i-1], $pY->[$i], $pY->[$i+1]);
		my $S = SinglePointIntegrateBySimpson(\@xx, \@yy);
		$y[$i+1] = $y[$i-1] + $S;
	}
	return \@y;
}

sub IntegrateBySimpson38WithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 3 or $nBlock % 3 != 0);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-3 ; $i += 3) {
		my $xi  = $x0  + $i * $h;
		my $xi1 = $xi  + $h;
		my $xi2 = $xi1 + $h;
		my $xi3 = $xi2 + $h;
		$S += 3.0 * &$pFunc($xi) + 9.0 * &$pFunc($xi1) + 9.0 * &$pFunc($xi2) + 3.0 * &$pFunc($xi3);
#		$S += 3.0 * &$pFunc($xi) + 9.0 * &$pFunc($xi) + 9.0 * &$pFunc($xi) + 3.0 * &$pFunc($xi);
#		$S += 24.0 * &$pFunc($xi);# + 9.0 * &$pFunc($xi) + 9.0 * &$pFunc($xi) + 3.0 * &$pFunc($xi);
	}

	return $S / 8.0 * $h;
}

sub IntegrateByBodeWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 4 or $nBlock % 4 != 0);

	my $h = ($x1 - $x0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-4 ; $i += 4) {
		my $xi  = $x0  + $i * $h;
		my $xi1 = $xi  + $h;
		my $xi2 = $xi1 + $h;
		my $xi3 = $xi2 + $h;
		my $xi4 = $xi3 + $h;
		$S += 14.0 * &$pFunc($xi) + 64.0 * &$pFunc($xi1) + 24.0 * &$pFunc($xi2)
				+ 64.0 * &$pFunc($xi3) + 14.0 * &$pFunc($xi4);
	}

	return $S / 45.0 * $h;
}

sub IntegrateByConstantStepBode
{
	my ($nData, $h, $pY) = @_;
	return undef if($nData < 2);

	my $py03 = IntegrateBySimpson([0.0, $h, 2.0*$h, 3.0*$h], [@$pY[0..4]]);
	my @y = @$py03;
	$h /= 45.0;
	for(my $i = 0 ; $i < $nData-4 ; $i++) {
		my $sum = 14.0 * $pY->[$i]   + 64.0 * $pY->[$i+1] + 24.0 * $pY->[$i+2]
				+ 64.0 * $pY->[$i+3] + 14.0 * $pY->[$i+4];
		$y[$i+4] = $y[$i] + $sum * $h;
	}
	return \@y;
}

sub IntegrateByRombergFunc
{
	my ($x) = @_;

	return exp(-($x**2));
#	return sqrt(1.0 - $x*$x);
}

sub IntegrateByRomberg
{
	my ($pFunc, $x0, $x1, $nMax, $EPS, $IsPrint) = @_;

	my $y = &$pFunc($x0);
print "y($x0)=$y\n" if($IsPrint);

	my @S;
	$S[0] = [];

	my $n = 1;
	my $h   = $x1 - $x0;
	my $y0  = &$pFunc($x0);
	my $y1  = &$pFunc($x1);
	$S[0][0] = ($y0 + $y1) / 2.0 * $h;
print "S00 = $S[0][0]\n" if($IsPrint);

	my $kfin;
	for(my $k = 1 ; $k < $nMax ; $k++) {
		$n *= 2;
print "k=$k (n=$n)\n" if($IsPrint);
		$h = ($x1 - $x0) / $n;

# 台形則の計算
		$S[$k][0] = 0.0;
		for(my $i = 0 ; $i < $n ; $i++) {
			my $xi  = $x0 + $i * $h;
#print "i=$i, xi=$xi\n";
			my $xi1 = $xi + $h;
			my $yi  = &$pFunc($xi);
			my $yi1 = &$pFunc($xi1);
			$S[$k][0] += ($yi + $yi1) / 2.0 * $h;
		}
print "  S[$k][0] = $S[$k][0]\n" if($IsPrint);

		my $d4 = 4.0;
		for(my $d = 1 ; $d <= $k ; $d++) {
#print "    d=$d\n";
			$S[$k][$d] = ($d4 * $S[$k][$d-1] - $S[$k-1][$d-1]) / ($d4 - 1.0);
print "    S[$k][$d] = $S[$k][$d]\n" if($IsPrint);

			$d4 *= 4.0;
		}
		
		my $delta = abs($S[$k][$k] - $S[$k-1][$k-1]);
		if($delta < $EPS) {
			$kfin = $k;
			last;
		}
	}

	my $mess;
	if(!defined $kfin) {
		$mess = "Not converged";
		$kfin = $nMax-1;
	}
	my $Sfin = $S[$kfin][$kfin];
	my $diff = ($kfin >= 1)? $Sfin - $S[$kfin-1][$kfin-1] : undef;

	return ($Sfin, $n, $kfin, $diff, $mess);
}

sub IntegrateByIMTWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock, $nGL, $nGLOrder) = @_;
	$nGL      = 2 if(!defined $nGL);
	$nGLOrder = 5 if(!defined $nGLOrder);

	return undef if($nBlock < 2);

# x = ($x1 - $x0) * x' + $x0)]
# [x = $x0, $x1] => [x' = 0, 1]
# S= Int(x, $x0, $x1) = ($x1 - $x0) * Int(x', 0, 1)

	my $Q = 0.00702985840661; # = int(exp(-1/t - 1/(1-t))

	my $h = 1.0 / $nBlock;
	my $S = 0.0;
	for(my $i = 1 ; $i <= $nBlock-1 ; $i++) {
#		my $ui   = $i * $h;
#		my $xi   = ($x1 - $x0) * $ui + $x0;
#		my $yi   = &$pFunc($xi);
## phi  = exp(-1.0 / $ui - 1.0 / (1.0 - $ui))
#		my $phi  = exp(-1.0 / $ui - 1.0 / (1.0 - $ui));
## phi' = [1/$ui/$ui - 1/(1-$ui)^2] * phi
#		my $ui1m = 1.0 - $ui;
#		my $phi1 = (1.0/$ui/$ui - 1.0/$ui1m/$ui1m) * $phi;
#
#		$S += $yi * $phi1;

		my $ui     = $i * $h;
		my $phiui = Algorism::IntegrateByGaussLegendreWithFunction(
			sub {
				my ($x) = @_;
				return exp(-1.0 / $x - 1.0 / (1.0-$x)) / $Q;
			},
			0.0, $ui, $nGL, $nGLOrder);
		my $xi  = ($x1 - $x0) * $phiui + $x0;
		my $yi   = &$pFunc($xi);
#print "i=$i: ui=$ui => phi(u)=$phiui => xi=$xi\n";
		$S += $yi * exp(-1.0 / $ui - 1.0 / (1.0 - $ui));
	}

	return $S / $Q * $h * ($x1 - $x0);
}

sub IntegrateByDoubleIntegralWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock, $u0, $u1) = @_;
	$u0 = -4.0 if(!defined $u0);
	$u1 =  4.0 if(!defined $u1);

	return undef if($nBlock < 2);

# x = [($x1 - $x0) * x'+ ($x0 + $x1)] / 2.0
# [x = $x0, $x1] => [x' = -1, 1]
# S= Int(x, $x0, $x1) = ($x1 - $x0) / 2.0 * Int(x', -1, 1)

	my $h = ($u1 - $u0) / $nBlock;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock ; $i++) {
		my $ui    = $u0 + $i * $h;
		my $phiui = tanh($pi / 2.0 * sinh($ui));
		my $xi    = (($x1 - $x0) * $phiui + $x0 + $x1) / 2.0;
#print "i=$i: ui=$ui => phi(u)=$phiui => xi=$xi\n";
		my $cosh = cosh($ui);
		my $sinh = sinh($ui);
		my $coshsinh = cosh($pi / 2.0 * sinh($ui));
		$S += &$pFunc($xi) * $cosh / ($coshsinh * $coshsinh);
	}

	return ($pi / 2.0) * $S * $h * ($x1 - $x0) / 2.0;
}

sub IntegrateByCubicSplineWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock) = @_;

	return undef if($nBlock < 4);

	my (@x, @y);
	my $h = ($x1 - $x0) / $nBlock;
	for(my $i = 0 ; $i <= $nBlock ; $i++) {
		$x[$i]  = $x0 + $i * $h;
		$y[$i] = &$pFunc($x[$i]);
	}

	return Algorism::IntegrateByCubicSpline(\@x, \@y);
}

sub IntegrateByCubicSpline
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	if($nData < 4) {
		return undef;
	}
	for(my $i = 0 ; $i < $nData-1 ; $i++) {
		if($pX->[$i+1] <= $pX->[$i]) {
			return undef;
		}
	}

	my ($ph, $psp) = subspl($pX, $pY);

	my $w = 0.0;
	for(my $i = 1 ; $i < $nData ; $i++) {
		my $hi = $ph->[$i];
		my $hi2 = $hi * $hi;
		$w += $hi * ( $pY->[$i] + $pY->[$i-1] 
					- $hi2 * ($psp->[$i-1] + $psp->[$i]) * 8.3333333333333e-2);
	}
	return $w * 0.5;
}

my %ZPGaussLegendreIntegration = (
	'2' => [-1.0 / sqrt(3), + 1.0 / sqrt(3)],
	'3' => [-sqrt(3.0/5.0), 0.0, +sqrt(3.0/5.0)],
	'4' => [-0.861136311594052, -0.339981043584856,      +0.339981043584856, +0.861136311594052],
#	'5' => [-0.906179845938663, -0.538469310105683, 0.0, +0.538469310105683, +0.906179845938663],
	'5' => [-0.906179845938663992797626878299, -0.538469310105683091036314420700, 0.0, 
			+0.538469310105683091036314420700, +0.906179845938663992797626878299],

	);
my %WGaussLegendreIntegration  = (
	'2' => [1.0, 1.0],
	'3' => [5.0/9.0, 8.0/9.0, 5.0/9.0],
	'4' => [0.347854845137453, 0.652145154862546, 0.652145154862546, 0.347854845137453],
#	'5' => [0.236226885056189, 0.478628670499366, 0.568888888888889, 0.478628670499366, 0.236226885056189],
	'5' => [0.236926885056189087514264040719, 0.478628670499366468041291514835, 0.568888888888888888888888888888, 
			0.478628670499366468041291514835, 0.236926885056189087514264040719]
	);
#4 点公式
#?0.861136311594052575223946488892 0.347854845137453857373063949221
#?0.339981043584856264802665759103 0.652145154862546142626936050778
#+0.339981043584856264802665759103 0.652145154862546142626936050778
#+0.861136311594052575223946488892 0.347854845137453857373063949221
#5 点公式
#?0.906179845938663992797626878299 0.236926885056189087514264040719
#?0.538469310105683091036314420700 0.478628670499366468041291514835
#0 0.568888888888888888888888888888
#+0.538469310105683091036314420700 0.478628670499366468041291514835
#+0.906179845938663992797626878299 0.236926885056189087514264040719
#6 点公式
#?0.932469514203152027812301554493 0.171324492379170345040296142172
#?0.661209386466264513661399595019 0.360761573048438607569833513837
#?0.238619186093196908630501721680 0.467913934572691047389870343989
#+0.238619186093196908630501721680 0.467913934572691047389870343989
#+0.661209386466264513661399595019 0.360761573048438607569833513837
#+0.932469514203152027812301554493 0.171324492379170345040296142172
#7 点公式
#?0.949107912342758524526189684047 0.129484966168869693270611432679
#?0.741531185599394439863864773280 0.279705391489276667901467771423
#?0.405845151377397166906606412076 0.381830050505118944950369775488
#0 0.417959183673469387755102040816
#+0.405845151377397166906606412076 0.381830050505118944950369775488
#+0.741531185599394439863864773280 0.279705391489276667901467771423
#+0.949107912342758524526189684047 0.129484966168869693270611432679

sub IntegrateByGaussLegendreWithFunction
{
	my ($pFunc, $x0, $x1, $nBlock, $nOrder, $pZP, $pW) = @_;

	return undef if($nBlock < 1);

	if(!defined $pZP) {
		$pZP = $ZPGaussLegendreIntegration{"$nOrder"};
		$pW  = $WGaussLegendreIntegration{"$nOrder"};
	}
	if(!defined $pZP or !defined $pW) {
		print "Error in Algorism::IntegrateByGaussLegendreWithFunction: Can not find zero points and weights for order [$nOrder].\n";
		return undef;
	}
#	if($nOrder == 5) {
#		print "\n\nWarning: Algorism::IntegrateByGaussLegendreWithFunction: G-L integration with nOrder=5 is not guaranteed.\n";
#		print "         The acuracy of the zero points and/or weight must be re-examined.\n\n\n";
#		
#	}
#for(my $i = 0 ; $i < $nOrder ; $i++) {
#	print "$i: Zero at $pZP->[$i]\tw=$pW->[$i]\n";
#}

	my $h  = ($x1 - $x0) / $nBlock;
	my $h2 = $h / 2.0;

	my $S = 0.0;
	for(my $i = 0 ; $i <= $nBlock-1 ; $i++) {
		my $xc = $x0 + ($i + 0.5) * $h;
		for(my $j = 0 ; $j < $nOrder ; $j++) {
			my $xj = $xc + $pZP->[$j] * $h2;
#print "i=$i: j=$j: xc=$xc: xj=$xj (ZP=$pZP->[$j])\n";
			$S += $pW->[$j] * &$pFunc($xj);
		}
	}
#print "S=$S\n";
#exit;
	return $S * $h2;
}

sub Differentiate2ndByIndex
{
	my ($i, $py, $h) = @_;
	return ($py->[$i+1] + $py->[$i-1] - 2.0 * $py->[$i]) / $h / $h;
}

sub Differentiate2nd
{
	my ($pFunc, $x, $h) = @_;
	return (&$pFunc($x + $h) - &$pFunc($x)) / $h;
	return (&$pFunc($x + $h) + &$pFunc($x - $h) - 2.0 * &$pFunc($x)) / $h / $h;
}

sub Differentiate2
{
	my ($pFunc, $x, $h) = @_;
	return (&$pFunc($x + $h) - &$pFunc($x)) / $h;
}

sub Differentiate3
{
	my ($pFunc, $x, $h) = @_;
	return (&$pFunc($x + $h) - &$pFunc($x - $h)) / 2.0 / $h;
}

sub Differentiate5
{
	my ($pFunc, $x, $h) = @_;
	return (-&$pFunc($x + 2*$h) + 8*&$pFunc($x + $h) - 8*&$pFunc($x - $h) + &$pFunc($x - 2*$h)) / 12.0 / $h;
}

sub Differentiate7
{
	my ($pFunc, $x, $h) = @_;
	return (&$pFunc($x + 3*$h) - 9*&$pFunc($x + 2*$h) + 45*&$pFunc($x + $h) 
			- 45*&$pFunc($x - $h) + 9*&$pFunc($x - 2*$h) - &$pFunc($x - 3*$h)) / 60.0 / $h;
}

sub Differentiate7WithIndex
{
	my ($pY, $i0, $h) = @_;
	return ($pY->[$i0 + 3] - 9*$pY->[$i0 + 2] + 45*$pY->[$i0 + 1] 
		- 45*$pY->[$i0 - 1] + 9*$pY->[$i0 - 2] - $pY->[$i0 - 3]) / 60.0 / $h;
}

sub DifferentiateByLinearAverage
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	return undef if($nData < 2);

	my @y;
	$y[0] = ($pY->[1] - $pY->[0]) / ($pX->[1] - $pX->[0]);
	$y[$nData-1] = ($pY->[$nData-1] - $pY->[$nData-2]) / ($pX->[$nData-1] - $pX->[$nData-2]);
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my $d1 = ($pY->[$i]   - $pY->[$i-1]) / ($pX->[$i]   - $pX->[$i-1]);
		my $d2 = ($pY->[$i+1] - $pY->[$i]  ) / ($pX->[$i+1] - $pX->[$i]  );
		$y[$i] = ($d1 + $d2) * 0.5;
	}
	return \@y;
}

sub DifferentiateByLagrange
{
	my ($pX, $pY, $x) = @_;

	my $nData = @$pX;
	if($nData < 2 or $x < $pX->[0] or $x >= $pX->[$nData-1]) {
		return undef;
	}
	for(my $i = 0 ; $i < $nData-1 ; $i++) {
		if($pX->[$i+1] <= $pX->[$i]) {
			return undef;
		}
	}

	my $w = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $sm0 = 0.0;
		my $sm2 = 1.0;
		my $xi = $pX->[$i];
		for(my $j = 0 ; $j < $nData ; $j++) {
			if($i != $j) {
				$sm2 *= ($xi - $pX->[$j]);
				my $sm1 = 1.0;
				for(my $k = 0 ; $k < $nData ; $k++) {
					if($k != $i and $k != $j) {
						$sm1 *= ($x - $pX->[$k]);
					}
				}
				$sm0 += $sm1;
			}
		}
		$w += $pY->[$i] * $sm0 / $sm2;
	}
	return $w;
}

sub DifferentiateByCubicSpline
{
	my ($pX, $pY, $x) = @_;

	my $nData = @$pX;
	if($nData < 2 or $x < $pX->[0] or $x >= $pX->[$nData-1]) {
		return undef;
	}
	for(my $i = 0 ; $i < $nData-1 ; $i++) {
		if($pX->[$i+1] <= $pX->[$i]) {
			return undef;
		}
	}

	my ($ph, $psp) = subspl($pX, $pY);

	for(my $i = 1 ; $i < $nData ; $i++) {
		if($pX->[$i-1] <= $x and $x <= $pX->[$i]) {
			my $hi  = $ph->[$i];
			my $hi2 = $hi * $hi;
			my $dxp2 = ($pX->[$i] - $x) * ($pX->[$i] - $x);
			my $dxm2 = ($x - $pX->[$i-1]) * ($x - $pX->[$i-1]);
			my $diff = ($psp->[$i-1] * ($hi2 - 3.0 * $dxp2)
					  + $psp->[$i] * (3.0 * $dxm2 - $hi2)
					  + 6.0 * ($pY->[$i] - $pY->[$i-1]) ) / (6.0 * $hi);
			return $diff;
		}
	}
	return undef;
}

sub DifferentiateByRichardsonExtrapolation
{
	my ($pFunc, $x, $h0, $nMax, $EPS, $IsPrint) = @_;

	my $h = $h0;
	my @D;
	my $delta;
	my $mfin;
	$D[0][0] = (&$pFunc($x + $h) - &$pFunc($x - $h)) / 2.0 / $h;
print "D00: $D[0][0]\n" if($IsPrint);
	for(my $m = 1 ; $m <= $nMax ; $m++) {
		$h *= 0.5;
		$D[0][$m] = (&$pFunc($x + $h) - &$pFunc($x - $h)) / 2.0 / $h;
print "D0$m: $D[0][$m]\n" if($IsPrint);
		for(my $k = $m-1 ; $k >= 0 ; $k--) {
			my $m1 = $m - $k;
			$D[$m1][$k] = (4**$m1 * $D[$m1-1][$k+1] - $D[$m1-1][$k]) / (4**$m1 - 1);
print "D$m1$k: $D[$m1][$k]\n" if($IsPrint);
		}

		$delta = abs($D[$m][0] - $D[$m-1][0]);
		if($delta < $EPS) {
			$mfin = $m;
			last;
		}

	}
	
	if(!$mfin) {
		$mfin = $nMax;
	}

	return ($D[$nMax][0], $mfin, $h0 * (0.5**$mfin), $delta, \@D);
}

sub DifferentiateByRichardsonExtrapolation_old
{
	my ($pFunc, $x, $h0, $nMax, $EPS, $IsPrint) = @_;

	my $h = $h0;
	my @D;
	my $delta;
	for(my $k = 0 ; $k <= $nMax ; $k++) {
		$D[0][$k] = (&$pFunc($x + $h) - &$pFunc($x - $h)) / 2.0 / $h;
print "D0$k: $D[0][$k]\n" if($IsPrint);
		$h /= 2.0;
	}
	my $mfin;
print "D00: $D[0][0]\n" if($IsPrint);
	for(my $m = 1 ; $m <= $nMax ; $m++) {
		for(my $k = 0 ; $k <= $nMax-$m ; $k++) {
			$D[$m][$k] = (4**$m * $D[$m-1][$k+1] - $D[$m-1][$k]) / (4**$m - 1);
print "D$m$k: $D[$m][$k]\n" if($IsPrint);
		}

		$delta = abs($D[$m][0] - $D[$m-1][0]);
		if($delta < $EPS) {
			$mfin = $m;
			last;
		}

	}
	
	if(!$mfin) {
		$mfin = $nMax;
	}

	return ($D[$nMax][0], $mfin, $h0 * (0.5**$mfin), $delta, \@D);
}

sub FindIndexByBinarySearch
{
	my ($pX, $x) = @_;
	return 0 if($x <= $pX->[0]);
	my $nData = @$pX;
	return $nData - 1 if($pX->[$nData-1] <= $x);

	my $idx0 = 0;
	my $idx1 = $nData - 1;
	my $x0 = $pX->[$idx0];
	my $x1 = $pX->[$idx1];
	while(1) {
		my $idxh = int( ($idx0 + $idx1) /2 );
		my $xh   = $pX->[$idxh];
#print "x: $x ($x0,$xh,$x1) => ";
#		return $idxh if(($x1 - $x) * ($xh - $x) == 0.0);
		if(($x1 - $x) * ($xh - $x) <= 0.0) {
			$idx0 = $idxh;
			$x0 = $pX->[$idx0];
		}
		else {
#		elsif(($xh - $x) * ($x0 - $x) < 0.0) {
			$idx1 = $idxh;
			$x1 = $pX->[$idx1];
		}
#		else {
#			return undef;
#		}
#print "($x0 [$idx0],$x1 [$idx1])\n";
		last if($idx1 - $idx0 <= 1);
	}
#exit if($x > 1140);
	return $idx0;
}

sub Interpolate
{
	my ($pX, $pY, $x) = @_;
	return $pY->[0] if($x <= $pX->[0]);
	my $nData = @$pX;
	return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

	my $idx  = FindIndexByBinarySearch($pX, $x);
	my $idx1 = $idx+1;
	my $k = ($pY->[$idx1] - $pY->[$idx]) / ($pX->[$idx1] - $pX->[$idx]);
	return $pY->[$idx] + $k * ($x - $pX->[$idx]);
}

sub InterpolateForConstantStepX
{
	my ($nData, $x0, $dx, $pY, $x) = @_;
	my $idx = int(($x - $x0) / $dx);
	return $pY->[0] if($idx < 0);
	return $pY->[$nData-1] if($idx >= $nData);
#return $pY->[$idx];
	$idx-- if($idx == $nData-1);

	my $idx1 = $idx+1;
#print "idx: $idx, $idx1 [$nData]\n";
	my $xa = $x0 + $idx*$dx;
	my $xb = $x0 + $idx1*$dx;
#print "x: $xa < $x < $xb\n";
#print " y[$idx1]=$pY->[$idx1], $pY->[$idx]\n";
	my $k = ($pY->[$idx1] - $pY->[$idx]) / ($xb - $xa);
#print "k=$k\n";
	return $pY->[$idx] + $k * ($x - $xa);
}

sub InterpolateByLagrange
{
	my ($pX, $pY, $x) = @_;

	my $nData = @$pX;
	if($nData < 3) {
#	if($nData < 3 or $x < $pX->[0] or $x >= $pX->[$nData-1]) {
		return undef;
	}

	my $w = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $w1 = 1.0;
		my $w2 = $pX->[$i];
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $w3 = $pX->[$j];
			if($i != $j) {
				$w1 *= ($x - $w3) / ($w2 - $w3);
			}
		}
		$w += $w1 * $pY->[$i];
	}
#print ("y=$w for $x [$pX->[0],$pX->[1],$pX->[2]] [$pY->[0],$pY->[1],$pY->[2]]\n");
#last if($x > 1150);
	return $w;
}

sub InterpolateBySimpson
{
	my ($pX, $pY, $x) = @_;
	return $pY->[0] if($x <= $pX->[0]);
	my $nData = @$pX;
	return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

	my $idx  = FindIndexByBinarySearch($pX, $x);
	$idx = 1        if($idx <= 1);
	$idx = $nData-2 if($idx >= $nData-3);
#my $idxm1=$idx-1;
#my $idxp1=$idx+1;
#print "idx=$idx\tx=$x ($pX->[$idxm1],$pX->[$idx],$pX->[$idxp1])\n";
	my $y = InterpolateByLagrange( [@$pX[$idx-1..$idx+1]], [@$pY[$idx-1..$idx+1]], $x);
#print "\t\ty=$y ($pY->[$idxm1],$pY->[$idx],$pY->[$idxp1])\n";
	return $y;
}

sub InterpolateBy4thPolynomial
{
	my ($pX, $pY, $x) = @_;
	return $pY->[0] if($x <= $pX->[0]);
	my $nData = @$pX;
	return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

	my $idx  = FindIndexByBinarySearch($pX, $x);
	$idx = 2        if($idx <= 2);
	$idx = $nData-3 if($idx >= $nData-4);
#my $idxm1=$idx-1;
#my $idxp1=$idx+1;
#print "idx=$idx\tx=$x ($pX->[$idxm1],$pX->[$idx],$pX->[$idxp1])\n";
	my $y = InterpolateByLagrange( [@$pX[$idx-2..$idx+2]], [@$pY[$idx-2..$idx+2]], $x);
#print "\t\ty=$y ($pY->[$idxm1],$pY->[$idx],$pY->[$idxp1])\n";
	return $y;
}

sub InterpolateByCubicSpline
{
	my ($pX, $pY, $x) = @_;

	my $nData = @$pX;
	return undef           if($nData < 3);
	return $pY->[0]        if($x <= $pX->[0]);
	return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

	for(my $i = 0 ; $i < $nData-1 ; $i++) {
		if($pX->[$i+1] <= $pX->[$i]) {
my $i1 = $i+1;
print "Warning in Algorism::InterpolateByCubicSpline: $i: d[i+1] must be > d[i] ($pX->[$i1] < $pX->[$i1])\n";
			return undef;
		}
	}

	my ($ph, $psp) = subspl($pX, $pY);

	for(my $i = 1 ; $i < $nData ; $i++) {
		if($pX->[$i-1] <= $x and $x < $pX->[$i]) {
			my $sm  = $psp->[$i-1];
			my $si  = $psp->[$i];
			my $hi  = $ph->[$i];
			my $hi2 = $hi * $hi;
			my $dxp = $pX->[$i] - $x;
			my $dxm = $x - $pX->[$i-1];
			my $v = $sm * $dxp * $dxp * $dxp
				  + $si * $dxm * $dxm * $dxm
				  + (6.0 * $pY->[$i-1] - $hi2 * $sm) * $dxp
				  + (6.0 * $pY->[$i]   - $hi2 * $si) * $dxm;
			return $v / $hi / 6.0;
		}
	}
	return undef;
}

sub subspl
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	my (@h, @sp, @du, @dc);
	$du[$nData-1] = 0.0;
	$dc[$nData-1] = 0.0;
	$sp[$nData-1] = 0.0;
	for(my $i = 1 ; $i < $nData ; $i++) {
		$h[$i] = $pX->[$i] - $pX->[$i-1];
	}
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my $hip1 = $pX->[$i+1] - $pX->[$i];
		$sp[$i] = 6.0 * (($pY->[$i+1] - $pY->[$i]) / ($h[$i] * $hip1)
			    - ($pY->[$i] - $pY->[$i-1]) / ($h[$i] * $h[$i]));
		$du[$i] = $h[$i+1] / $h[$i];
		$dc[$i] = 2.0 * (1.0 + $du[$i]);
	}
	$dc[$nData-1] += 1.0;
	$dc[$nData-2] += $h[$nData-1] / $h[$nData-2];
	$du[1] /= $dc[1];
	$sp[1] /= $dc[1];
	for(my $i = 2 ; $i < $nData ; $i++) {
#print "$i: $dc[$i], $du[$i-1]\n";
		my $g = 1.0 / ($dc[$i] - $du[$i-1]);
		$du[$i] *= $g;
		$sp[$i] = ($sp[$i] - $sp[$i-1]) * $g;
	}
	for(my $i = $nData-2 ; $i >= 1 ; $i--) {
#print "$i: $sp[$i+1], $du[$i]\n";
		$sp[$i] -= $sp[$i+1] * $du[$i];
	}
	$sp[0] = $sp[1];
	$sp[$nData-1] = $sp[$nData-2];
	return (\@h, \@sp);
}

sub dydxFunc
{
	my ($x, $y) = @_;

# dy / dx = y

	return $y;
}

sub SolveByEuler
{
	my ($pdydxFunc, $x0, $x1, $h, $pYinitial) = @_;
	
	my (@x, @y);

	$x[0] = $x0;
	$y[0] = $pYinitial->[0];
	for(my $i = 1, my $x = $x0 + $h ; $x <= $x1 ; $i++, $x = $x0 + $i * $h) {
		$x[$i] = $x;
		my $dydx = &$pdydxFunc($x[$i-1], $y[$i-1]);
		$y[$i] = $y[$i-1] + $h * $dydx;
#print "$i: ($x[$i], $y[$i])\n";
	}
	return (\@x, \@y);
}

sub dydxMultiFunc
{
	my ($idx, $x, $py) = @_;

# dy / dx = y
	if($idx == 0) {
		return $py->[0];
	}
	else {
		print "error: Function index idx=$idx is not found.\n";
		exit;
	}
}

sub SolveSimultaneousDifferentialEquationsByEuler
{
	my ($nEq, $pdydxMultiFunc, $x0, $x1, $h, $ppYinitial) = @_;
	
	my (@x, @py);

	$x[0] = $x0;
	for(my $ie = 0 ; $ie < $nEq ; $ie++) {
		$py[$ie] = [];
		$py[$ie][0] = $ppYinitial->[$ie][0];
#print "ie=$ie: y(x0) = $py[$ie][0]\n";
	}

	my $x;
	for(my $i = 1 ; $x <= $x1 ; $i++) {
		$x[$i] = $x = $x0 + $i * $h;
		my @YArray;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1];
		}
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $dydx = &$pdydxMultiFunc($ie, $x[$i-1], \@YArray);
			$py[$ie][$i] = $py[$ie][$i-1] + $h * $dydx;
#print "i=$i, ie=$ie: ($x[$i], $py[$ie][$i]) (dy/dx=$dydx)\n";
		}
	}
	return (\@x, \@py);
}

sub SolveByHeun
{
	my ($pdydxFunc, $x0, $x1, $h, $pYinitial) = @_;
	
	my (@x, @y);

	$x[0] = $x0;
	$y[0] = $pYinitial->[0];
	for(my $i = 1, my $x = $x0 + $h ; $x <= $x1 ; $i++, $x = $x0 + $i * $h) {
		$x[$i] = $x;
		my $dydx = &$pdydxFunc($x[$i-1], $y[$i-1]);
		my $k1 = $h * $dydx;
		my $k2 = $h * &$pdydxFunc($x[$i], $y[$i-1] + $k1);
		$y[$i] = $y[$i-1] + ($k1 + $k2) / 2.0;
	}
	return (\@x, \@y);
}

sub SolveSimultaneousDifferentialEquationsByHeun
{
	my ($nEq, $pdydxMultiFunc, $x0, $x1, $h, $ppYinitial) = @_;
	
	my (@x, @py);

	$x[0] = $x0;
	for(my $ie = 0 ; $ie < $nEq ; $ie++) {
		$py[$ie] = [];
		$py[$ie][0] = $ppYinitial->[$ie][0];
#print "ie=$ie: y(x0) = $py[$ie][0]\n";
	}

	my $x;
	for(my $i = 1 ; $x <= $x1 ; $i++) {
		$x[$i] = $x = $x0 + $i * $h;
		my @YArray;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1];
		}
		my @k1;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $dydx = &$pdydxMultiFunc($ie, $x[$i-1], \@YArray);
			$k1[$ie] = $h * $dydx;
		}
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k1[$ie];
		}
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $k2 = $h * &$pdydxMultiFunc($ie, $x[$i], \@YArray);
			$py[$ie][$i] = $py[$ie][$i-1] + ($k1[$ie] + $k2) / 2.0;
#print "i=$i, ie=$ie: ($x[$i], $py[$ie][$i]) (dy/dx=$dydx)\n";
		}
	}
	return (\@x, \@py);
}

sub SolveByRungekutta3
{
	my ($pdydxFunc, $x0, $x1, $h, $pYinitial) = @_;
	
	my (@x, @y);

	$x[0] = $x0;
	$y[0] = $pYinitial->[0];
	for(my $i = 1, my $x = $x0 + $h ; $x <= $x1 ; $i++, $x = $x0 + $i * $h) {
		$x[$i] = $x;
		my $dydx = &$pdydxFunc($x[$i-1], $y[$i-1]);
		my $k1 = $h * $dydx;
		my $k2 = $h * &$pdydxFunc($x[$i-1] + $h / 2.0, $y[$i-1] + $k1 / 2.0);
		my $k3 = $h * &$pdydxFunc($x[$i],            , $y[$i-1] + $k2 + $k2 - $k1);
		$y[$i] = $y[$i-1] + ($k1 + 4.0 * $k2 + $k3) / 6.0;
	}
	return (\@x, \@y);
}

sub SolveSimultaneousDifferentialEquationsByRungekutta3
{
	my ($nEq, $pdydxMultiFunc, $x0, $x1, $h, $ppYinitial) = @_;
	
	my (@x, @py);

	$x[0] = $x0;
	for(my $ie = 0 ; $ie < $nEq ; $ie++) {
		$py[$ie] = [];
		$py[$ie][0] = $ppYinitial->[$ie][0];
#print "ie=$ie: y(x0) = $py[$ie][0]\n";
	}

	my $x;
	for(my $i = 1 ; $x <= $x1 ; $i++) {
		$x[$i] = $x = $x0 + $i * $h;

		my @YArray;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1];
		}
		my @k1;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $dydx = &$pdydxMultiFunc($ie, $x[$i-1], \@YArray);
			$k1[$ie] = $h * $dydx;
		}

		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k1[$ie] / 2.0;
		}
		my @k2;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$k2[$ie] = $h * &$pdydxMultiFunc($ie, $x[$i-1] + $h / 2.0, \@YArray);
		}

		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k2[$ie] + $k2[$ie] - $k1[$ie];
		}
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $k3 = $h * &$pdydxMultiFunc($ie, $x[$i], \@YArray);
			$py[$ie][$i] = $py[$ie][$i-1] + ($k1[$ie] + 4.0 * $k2[$ie] + $k3) / 6.0;
#print "i=$i, ie=$ie: ($x[$i], $py[$ie][$i]) (dy/dx=$dydx)\n";
		}


	}
	return (\@x, \@py);
}

sub SolveByRungekutta4
{
	my ($pdydxFunc, $x0, $x1, $h, $pYinitial) = @_;
	
	my (@x, @y);

	$x[0] = $x0;
	$y[0] = $pYinitial->[0];
	for(my $i = 1, my $x = $x0 + $h ; $x <= $x1 ; $i++, $x = $x0 + $i * $h) {
		$x[$i] = $x;
		my $dydx = &$pdydxFunc($x[$i-1], $y[$i-1]);
		my $k1 = $h * $dydx;
		my $k2 = $h * &$pdydxFunc($x[$i-1] + $h / 2.0, $y[$i-1] + $k1 / 2.0);
		my $k3 = $h * &$pdydxFunc($x[$i-1] + $h / 2.0, $y[$i-1] + $k2 / 2.0);
		my $k4 = $h * &$pdydxFunc($x[$i]             , $y[$i-1] + $k3);
		$y[$i] = $y[$i-1] + ($k1 + 2.0 * $k2 + 2.0 * $k3 + $k4) / 6.0;
	}
	return (\@x, \@y);
}

sub SolveSimultaneousDifferentialEquationsByRungekutta4
{
	my ($nEq, $pdydxMultiFunc, $x0, $x1, $h, $ppYinitial) = @_;
	
	my (@x, @py);

	$x[0] = $x0;
	for(my $ie = 0 ; $ie < $nEq ; $ie++) {
		$py[$ie] = [];
		$py[$ie][0] = $ppYinitial->[$ie][0];
#print "ie=$ie: y(x0) = $py[$ie][0]\n";
	}

	my $x;
	for(my $i = 1 ; $x <= $x1 ; $i++) {
		$x[$i] = $x = $x0 + $i * $h;

		my @YArray;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1];
		}
		my @k1;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $dydx = &$pdydxMultiFunc($ie, $x[$i-1], \@YArray);
			$k1[$ie] = $h * $dydx;
		}

		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k1[$ie] / 2.0;
		}
		my @k2;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$k2[$ie] = $h * &$pdydxMultiFunc($ie, $x[$i-1] + $h / 2.0, \@YArray);
		}

		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k2[$ie] / 2.0;
		}
		my @k3;
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$k3[$ie] = $h * &$pdydxMultiFunc($ie, $x[$i-1] + $h / 2.0, \@YArray);
		}

		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			$YArray[$ie] = $py[$ie][$i-1] + $k3[$ie];
		}
		for(my $ie = 0 ; $ie < $nEq ; $ie++) {
			my $k4 = $h * &$pdydxMultiFunc($ie, $x[$i], \@YArray);
			$py[$ie][$i] = $py[$ie][$i-1] + ($k1[$ie] + 2.0 * $k2[$ie] + 2.0 * $k3[$ie] + $k4) / 6.0;
#print "i=$i, ie=$ie: ($x[$i], $py[$ie][$i]) (dy/dx=$dydx)\n";
		}


	}
	return (\@x, \@py);
}

sub SolveByVerlet
{
	my ($pdydxFunc, $x0, $x1, $h, $ppYinitial) = @_;
	
	my (@x, @y, @y2);
	
	my $kh2 = 1.0 / 2.0 / $h;
	my $h2  = $h * $h;

	$x[0]  = $x0;
	$y[0]  = $ppYinitial->[0][0];
	$y2[0] = $ppYinitial->[1][0];
	$x[1]  = $x0 + $h;
	$y[1]  = $ppYinitial->[0][1];
	$y2[1] = $ppYinitial->[1][1];
#	$y2[0] = ($y[1] - $y[0]) / $h;
	my $i = 2;
	for(my $x = $x0 + 2.0 * $h ; $x <= $x1 + $h ; $i++, $x = $x0 + $i * $h) {
		$x[$i]    = $x;
		$y[$i]    = 2.0 * $y[$i-1] - $y[$i-2] + $h2 * &$pdydxFunc($x[$i-1], [$y[$i-1], $y2[$i-1]]);
		$y2[$i-1] = ($y[$i] - $y[$i-2]) * $kh2;
	}
	$y2[$i] = ($y[$i] - $y[$i-1]) / $h;

	return (\@x, [\@y, \@y2]);
}

# Sande-Tukey Algorism (科学計測のための波形データ処理)
# ar[n], ai[n]: Real and imaginary parts of input data.
# n: datum number
# $mode: 'forward': Fourier transform. 'reverse': Reverse Fourier transform
sub FFT_ST # proper operation not confirmed
{
	my ($n, $pX, $ar_org, $ai_org, $mode) = @_;
	$ar_org = [] if(!defined $ar_org);
	$ai_org = [] if(!defined $ai_org);
	my $ar = [@$ar_org];
	my $ai = [@$ai_org];

	if($n < 2) {
print "Error in Algorism::FFT: n ($n) is smaller than 2.\n";
		return undef;
	}

	my $nf = $n / 2;
	my $m = int(log($n) / log(2.0) + 0.001);
#print "iter=$m (n=$n) (", 2**$m, ")\n";
	if(2**$m != $n) {
print "Error in Algorism::FFT: n ($n) must be 2^k.\n";
		return undef;
	}

	my (@S, @C);
	my $A = 0;
	my $B = $pi * 2.0 / $n;
	for(my $i= 0 ; $i <= $nf ; $i++) {
		if($mode eq 'reverse') {
			$S[$i] = sin($A);
			$C[$i] = cos($A);
		}
		else {
			$S[$i] = -sin($A);
			$C[$i] = cos($A);
		}
		$A += $B;
	}

# FFT
	my $j = 0;
	my $p = 0;

	my $l = $n;
	my $k = 0;
	my $h = 1;
	for(my $g = 1 ; $g <= $m ; $g++) {
		$l = $l / 2;
		$k = 0;
		for(my $q = 1 ; $q <= $h ; $q++) {
			$p = 0;
			for(my $i = $k ; $i <= $l + $k - 1 ; $i++) {
				$j = $i + $l;
				$A = $ar->[$i] - $ar->[$j];
				$B = $ai->[$i] - $ai->[$j];
				$ar->[$i] += $ar->[$j];
				$ai->[$i] += $ai->[$j];
				if($p == 0) {
					$ar->[$j] = $A;
					$ai->[$j] = $B;
				}
				else {
					$ar->[$j] = $A * $C[$p] + $B * $S[$p];
					$ai->[$j] = $B * $C[$p] - $A * $S[$p];
				}
				$p += $h;
			}
			$k += $l + $l;
		}
		$h += $h;
	}

# bit reversal
	$j = $nf;
	for(my $i = 1; $i <= $n-1; $i++) {
		$k = $n;
		if($j < $i) {
			my $r = $ar->[$i];
			$ar->[$i] = $ar->[$j];
			$ar->[$j] = $r;
			$r = $ai->[$i];
			$ai->[$i] = $ai->[$j];
			$ai->[$j] = $r;
		}
		$k /= 2;
		while($j >= $k) {
			$j -= $k;
			$k /= 2;
		}
		$j += $k;
	}

	my $df = 1.0 / ($pX->[1] + $pX->[$n-1] - 2.0 * $pX->[0]);
	my @f;
	for(my $i = 0; $i < $nf; $i++) {
		$f[$i] = $df * ($i+1);
	}

	if($mode eq 'reverse') {
		my $w = -1.0 / $n;
		for(my $i = 0; $i < $nf; $i++) {
			$ar->[$i] = $w * $ar->[$i];
			$ai->[$i] = $w * $ai->[$i];
		}
	}

#print "a=[$ar, $ai]\n";
	return($nf, \@f, $ar, $ai);
}

# ar[n], ai[n]: Real and imaginary parts of input data.
# n: datum number
# $mode: 'forward': Fourier transform. 'reverse': Reverse Fourier transform
sub FFT
{
	my ($n, $pX, $ar_org, $ai_org, $mode) = @_;
	$ar_org = [] if(!defined $ar_org);
	$ai_org = [] if(!defined $ai_org);
	my $ar = [@$ar_org];
	my $ai = [@$ai_org];

	if($n < 2) {
print "Error in Algorism::FFT: n ($n) is smaller than 2.\n";
		return undef;
	}

	my $nf = $n / 2;
	my $iter = int(log($n) / log(2.0) + 0.001);
#print "iter=$iter (n=$n)\n";

	my $j = 1;
	for(my $i = 0; $i < $iter; $i++) {
		$j *= 2;
	}
	if($n != $j) {
print "Error in Algorism::FFT: n ($n) must be 2^k.\n";
		return undef;
	}

	my $sign = 1.0;
	if($mode eq 'reverse') {
		$sign = -1.0;
	}

	my $df = 1.0 / ($pX->[1] + $pX->[$n-1] - 2.0 * $pX->[0]);

	my ($j1, $j2, $dr1, $dr2, $di1, $di2, $tr, $ti);
	my $xp2 = $n;
	for(my $it = 0; $it < $iter; $it++) {
		my $xp = $xp2;
		$xp2 = $xp / 2;
		my $w = $pi / $xp2;
		for(my $k = 0; $k < $xp2; $k++) {
			my $arg = $k * $w;
			my $wr = cos($arg);
			my $wi = $sign * sin($arg);
			my $i = $k - $xp;
			for($j = $xp; $j <= $n; $j += $xp) {
				$j1 = $j  + $i;
				$j2 = $j1 + $xp2;
				$dr1 = $ar->[$j1];
				$dr2 = $ar->[$j2];
				$di1 = $ai->[$j1];
				$di2 = $ai->[$j2];
				$tr  = $dr1 - $dr2;
				$ti  = $di1 - $di2;
				$ar->[$j1] = $dr1 + $dr2;
				$ai->[$j1] = $di1 + $di2;
				$ar->[$j2] = $tr * $wr - $ti * $wi;
				$ai->[$j2] = $ti * $wr + $tr * $wi;
			}
		}
	}
	$j1 = $nf;
	$j2 = $n - 1;
	$j = 1;
	for(my $i = 1; $i <= $j2; $i++) {
		if($i < $j) {
			my $im1 = $i - 1;
			my $jm1 = $j - 1;
			my $tr  = $ar->[$jm1];
			my $ti  = $ai->[$jm1];
			$ar->[$jm1] = $ar->[$im1];
			$ai->[$jm1] = $ai->[$im1];
			$ar->[$im1] = $tr;
			$ai->[$im1] = $ti;
		}
		my $k = $j1;
		while($k < $j) {
			$j -= $k;
			$k /= 2;
		}
		$j += $k;
	}

	my @f;
	for(my $i = 0; $i < $nf; $i++) {
		$f[$i] = $df * ($i+1);
	}

	if($mode eq 'reverse') {
		my $w = -1.0 / $n;
		for(my $i = 0; $i < $nf; $i++) {
			$ar->[$i] = $w * $ar->[$i];
			$ai->[$i] = $w * $ai->[$i];
		}
	}

	return($nf, \@f, $ar, $ai);
}

sub DFT
{
	my ($n, $pX, $pR, $pI, $mode) = @_;

#	my $df = 1.0 / ($pX->[$n-1] - $pX->[0]);
	my $df = 1.0 / ($pX->[1] + $pX->[$n-1] - 2.0 * $pX->[0]);

	my (@f, @yr, @yi);
	my $C = 1.0;
	$C = -1.0 if($mode eq 'inverse');
	
	my (@wr, @wi);
	$wr[0] = 1.0;
	$wi[0] = 0.0;
	my $phi = $C * 2.0*$pi / $n;
	my $wr = cos($phi);
	my $wi = sin($phi);
	$wr[1] = $wr;
	$wi[1] = $wi;
	for(my $k = 2 ; $k < $n ; $k++) {
#	for(my $k = 0 ; $k < $n ; $k++) {
#		my $phi = $C * 2.0*$pi / $n * $k;
#		$wr[$k] = cos($phi);
#		$wi[$k] = sin($phi);
# w^k = (wr(k-1) + i * wi(k-1)) * (wr + i * wi)
#     = (wr(k-1) * wr - wi(k-1) * wi) + i * (wi(k-1) * wr + wr(k-1) * wi)
		$wr[$k] = $wr[$k-1] * $wr - $wi[$k-1] * $wi;
		$wi[$k] = $wr[$k-1] * $wi + $wi[$k-1] * $wr;
	}

	for(my $k = 0 ; $k < $n/2 ; $k++) {
		$f[$k]  = $df * ($k+1);
		$yr[$k] = 0.0;
		$yi[$k] = 0.0;
		for(my $i = 0 ; $i < $n ; $i++) {
			my $ik = ($i * $k) % $n;
			my $wr = $wr[$ik];
			my $wi = $wi[$ik];
# x*exp(2pik) = (R + i I)(cos(phi) + i sin(phi))
#      Rcos(phi) - Isin(phi) + i(Rsin(phi) + Icos(phi))
			my $R = ($pR)? $pR->[$i] : 0.0;
			my $I = ($pI)? $pI->[$i] : 0.0;
			$yr[$k] += $R * $wr - $I * $wi;
			$yi[$k] += $R * $wi + $I * $wr;
		}
	}
	return ($n/2, \@f, \@yr, \@yi);
}

sub DFT2
{
	my ($n, $pX, $pR, $pI, $mode) = @_;

#	my $df = 1.0 / ($pX->[$n-1] - $pX->[0]);
	my $df = 1.0 / ($pX->[1] + $pX->[$n-1] - 2.0 * $pX->[0]);

	my (@f, @yr, @yi);
	my $C = 1.0;
	$C = -1.0 if($mode eq 'inverse');

	for(my $k = 0 ; $k < $n/2 ; $k++) {
		$f[$k]  = $df * ($k+1);
		$yr[$k] = 0.0;
		$yi[$k] = 0.0;
		for(my $i = 0 ; $i < $n ; $i++) {
			my $phi  = $C * 2.0*$pi / $n * $i*$k;
			my $cosp = cos($phi);
			my $sinp = sin($phi);
# x*exp(2pik) = (R + i I)(cos(phi) + i sin(phi))
#      Rcos(phi) - Isin(phi) + i(Rsin(phi) + Icos(phi))
			my $R = ($pR)? $pR->[$i] : 0.0;
			my $I = ($pI)? $pI->[$i] : 0.0;
			$yr[$k] += $R * $cosp - $I * $sinp;
			$yi[$k] += $R * $sinp + $I * $cosp;
		}
	}
	return ($n/2, \@f, \@yr, \@yi);
}

sub randomGauss
{
	my ($avg, $sigma, $seed) = @_;
	$avg   = 0.0 if(!defined $avg);
	$sigma = 1.0 if(!defined $sigma);

	srand($seed) if(defined $seed);
	my $x = rand(1.0);
 	$x = 1.0e-100 if($x == 0.0);
	my $y = rand(1.0);
	my $z = sqrt(-2.0*log($x)) * sin(2.0 * $pi * $y);
	return $avg + $sigma * $z;
}

sub randomExp
{
	my ($rambda, $seed) = @_;
	$rambda = 1 if(!defined $rambda);

	srand($seed) if(defined $seed);
	my $x = rand(1.0);
	$x = 1.0e-100 if($x == 0.0);
	return -log($x) / $rambda;
}

#NTSCで使われているYIQカラーモデルに基づく方法
#  松田晃一、HTML5+JavaScriptによる画像・動画像処理入門、、カットシステム (2014)
sub RGBtoGray
{
	my ($R, $G, $B) = @_;
	return int(0.3 * $R + 0.59 * $G + 0.11 * $B);
}

1;