#!/usr/local/bin/perl

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2005-01-11
#
# This is an example of Lagrange interpolation, using degree 3.
#
# Given x_0 .. x_3, y_0 .. y_3, compute the polynomial
#
#     (x - x_1)(x - x_2)(x - x_3)
# y_0 ---------------------------- +
#     (x_0-x_1)(x_0-x_2)(x_0-x_3)
#
#     (x - x_0)(x - x_2)(x - x_3)
# y_1 ---------------------------- +
#     (x_1-x_0)(x_1-x_2)(x_1-x_3)
#
#     (x - x_0)(x - x_1)(x - x_3)
# y_2 ---------------------------- +
#     (x_2-x_0)(x_2-x_1)(x_2-x_3)
#
#     (x - x_0)(x - x_1)(x - x_2)
# y_3 ----------------------------
#     (x_3-x_0)(x_3-x_1)(x_3-x_2)
#
# Program output:
#
#   0.00000 -0.34510  1.00000 -1.3451e+00
#   0.00391 -0.32857  1.00391 -1.3325e+00
#   0.00781 -0.31211  1.00784 -1.3200e+00
#   ...
#   0.98828  2.68467  2.68661 -1.9441e-03
#   0.99219  2.69585  2.69713 -1.2788e-03
#   0.99609  2.70705  2.70768 -6.3083e-04
#   1.00000  2.71828  2.71828  0.0000e+00
#   1.00391  2.72953  2.72892  6.1388e-04
#   1.00781  2.74081  2.73960  1.2110e-03
#   1.01172  2.75212  2.75032  1.7915e-03
#   ...
#   1.23828  3.46298  3.44968  1.3299e-02
#   1.24219  3.47640  3.46318  1.3215e-02
#   1.24609  3.48986  3.47674  1.3125e-02
#   1.25000  3.50337  3.49034  1.3029e-02
#   1.25391  3.51693  3.50400  1.2927e-02
#   1.25781  3.53054  3.51772  1.2819e-02
#   1.26172  3.54419  3.53149  1.2705e-02
#   ...
#   1.48828  4.43016  4.42948  6.8910e-04
#   1.49219  4.44727  4.44681  4.5811e-04
#   1.49609  4.46445  4.46422  2.2838e-04
#   1.50000  4.48169  4.48169  0.0000e+00
#   1.50391  4.49900  4.49923 -2.2694e-04
#   1.50781  4.51639  4.51684 -4.5236e-04
#   1.51172  4.53384  4.53452 -6.7618e-04
#   ...
#   1.98828  7.30223  7.30297 -7.4339e-04
#   1.99219  7.33106  7.33155 -4.9811e-04
#   1.99609  7.36000  7.36025 -2.5028e-04
#   2.00000  7.38906  7.38906  0.0000e+00
#   2.00391  7.41823  7.41798  2.5265e-04
#   2.00781  7.44752  7.44701  5.0758e-04
#   2.01562  7.50644  7.50542  1.0239e-03
#   ...
#   2.48828 12.04297 12.04056  2.4043e-03
#   2.49219 12.08932 12.08769  1.6277e-03
#   2.49609 12.13583 12.13500  8.2641e-04
#   2.50000 12.18249 12.18249  0.0000e+00
#   2.50391 12.22932 12.23017 -8.5189e-04
#   2.50781 12.27631 12.27804 -1.7296e-03
#   2.51172 12.32346 12.32610 -2.6336e-03
#   ...
#   2.98828 19.39450 19.85153 -4.5704e-01
#   2.99219 19.46416 19.92923 -4.6507e-01
#   2.99609 19.53404 20.00723 -4.7319e-01
#
# Note that the fit is good at the interpolation points, but worse in
# between them and even worse at the edges.
#
# When the data is plotted, examine the error curve and note that it has
# appoximately the classic W shape of a quartic.  This makes sense, since we've
# used a cubic fit; the next-order term is quartic.
# ----------------------------------------------------------------

# Input mesh
$x0 = 0.20;
$x1 = 0.40;
$x2 = 0.60;
$x3 = 0.80;

# Function values to be estimated
$y0 = f($x0);
$y1 = f($x1);
$y2 = f($x2);
$y3 = f($x3);

# Interpolation mesh
$xlo = 0.0;
$xhi = 1.0;

#$nx  = 400;
#$nx  = $ARGV[0] if (@ARGV);
#$dx  = ($xhi - $xlo) / $nx;

# Alternatively:
$dx  = 1.0 / 64;

$c0 = $y0 / (          ($x0-$x1)*($x0-$x2)*($x0-$x3));
$c1 = $y1 / (($x1-$x0)          *($x1-$x2)*($x1-$x3));
$c2 = $y2 / (($x2-$x0)*($x2-$x1)          *($x2-$x3));
$c3 = $y3 / (($x3-$x0)*($x3-$x1)*($x3-$x2));

for ($x = $xlo; $x < $xhi; $x += $dx) {
	$actual = f($x);
	$estimated =
		$c0 *          ($x-$x1)*($x-$x2)*($x-$x3) +
		$c1 * ($x-$x0)         *($x-$x2)*($x-$x3) +
		$c2 * ($x-$x0)*($x-$x1)         *($x-$x3) +
		$c3 * ($x-$x0)*($x-$x1)*($x-$x2) ;

	printf "%10.5f %10.5f %10.5f %9.4e\n",
		$x, $estimated, $actual, $estimated - $actual;
}

# ----------------------------------------------------------------
sub f {
	my ($x) = @_;
	return atan2($x, 1.0);
}