#!/usr/bin/perl -w

use lib "d:/Programs/Perl/lib";

use strict;
use Sci::Science;

my $n = 20;
my $pi = Sci::pi();

my @x;
my @y;
my $xx = 1.0;
for(my $i = 0 ; $i <= $n ; $i++) {
	$x[$i] = $pi * $i / $n;
	$y[$i] = sin($x[$i]);
}

my $v = InterpolateBy3rdSpline(\@x, \@y, $xx);
my $z = sin($xx);
print "x=$xx, v=$v, z=$z\n";

exit;

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

	my $nData = @$pX;
	if($nData < 3 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 $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] = $x[$i] - $x[$i-1];
	}
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my $hip1 = $x[$i+1] - $x[$i];
		$sp[$i] = 6.0 * (($y[$i+1] - $y[$i]) / ($h[$i] * $hip1)
			    - ($y[$i] - $y[$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);
}