#!/usr/bin/perl -w

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

use strict;
use Sci::Science;

my $n = 10;
my $pi = Sci::pi();
my $xx = 0.5;

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

my $diff = DifferentiateByLagrange(\@x, \@y, $xx);

print "Differentiate by Lagrange\n";
print "  dy/dx=$diff at x=$xx\n";

exit;

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

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