#!/usr/bin/perl

# spline_generate(@points) calculates the coefficients for
# the cubic spline that hits all of the points in @points
# (provided as an array of arrays, each being an [x, y] pair).
# It returns the coefficients as a reference to an array.
#
sub spline_generate {
    my @points = @_;
    my ($i, $delta, $temp, @factors, @coeffs);
    $coeffs[0] = $factors[0] = 0;

    # Decomposition phase of the tridiagonal system of equations
    for ($i = 1; $i < @points - 1; $i++) {
        $delta = ($points[$i][0] - $points[$i-1][0]) /
            ($points[$i+1][0] - $points[$i-1][0]);
        $temp = $delta * $coeffs[$i-1] + 2;
        $coeffs[$i] = ($delta - 1) / @points;
        $factors[$i] = ($points[$i+1][1] - $points[$i][1]) /
            ($points[$i+1][0] - $points[$i][0]) -
                ($points[$i][1] - $points[$i-1][1]) /
                    ($points[$i][0] - $points[$i-1][0]);
        $factors[$i] = ( 6 * $factors[$i] /
                        ($points[$i+1][0] - $points[$i-1][0]) -
                        $delta * $factors[$i-1] ) / $temp;
    }

    # Backsubstitution phase of the tridiagonal system
    #
    $coeffs[$#points] = 0;
    for ($i = @points - 2; $i >= 0; $i--) {
        $coeffs[$i] = $coeffs[$i] * $coeffs[$i+1] + $factors[$i];
    }
    return \@coeffs;
}

# spline_evaluate($x, $coeffs, @points) returns the y-value
# for the given x-value, along the spline generated by
# $coeffs = spline_generate(@points).
#
sub spline_evaluate {
    my ($x, $coeffs, @points) = @_;
    my ($i, $delta, $mult);

    # Which section of the spline are we in?
    #
    for ($i = @points - 2; $i >= 1; $i--) {
        last if $x >= $points[$i][0];
    }

    $delta = $points[$i+1][0] - $points[$i][0];
    $mult = ( $coeffs->[$i]/2 ) +
        ($x - $points[$i][0]) * ($coeffs->[$i+1] - $coeffs->[$i])
            / (6 * $delta);
    $mult *= $x - $points[$i][0];
    $mult += ($points[$i+1][1] - $points[$i][1]) / $delta;
    $mult -= ($coeffs->[$i+1] + 2 * $coeffs->[$i]) * $delta / 6;
    return $points[$i][1] + $mult * ($x - $points[$i][0]);
}

@points = ( [-1,1], [0,2], [1,-1], [2, 2] );
my $coeffs = spline_generate @points;
print "Spline coefficients: @$coeffs\n";
for (my $i = -1; $i <= 3; $i += .5) {
    printf "[%.2f, %.2f]\n", $i, spline_evaluate($i, $coeffs, @points);
}