#!/usr/bin/perl

sub poly_fit {
    my @points = @_;
    my ($i, $j, $y, @solution);
    for ($i = 0; $i < @points; $i++) {
        ($x, $y) = ( @ {$points[$i]} );

        for ($j = 0; $j < @points-1; $j++) {
            $points[$i][$j] = $x ** (@points - $j - 1);
        }
        $points[$i][@points-1] = 1;
        $points[$i][@points]   = $y;
    }
    return linear_solve( @points );
}

@solution = poly_fit( [2,7000], [3,6000], [4,9000] );
print "@solution\n";                         # 2000 -11000 21000

use Math::MatrixReal;  # http://www.perl.com/CPAN/modules/by-module/Math

sub linear_solve {
    my @equations = @_;
    my ($i, $j, $solution, @solution, $dimension, $base_matrix);

    # Create $matrix, representing the lefthand side of our equations.
    #
    my $matrix = new Math::MatrixReal( scalar @equations,
                                       scalar @equations );

    # Create $vector, representing the y values.
    my $vector = new Math::MatrixReal( scalar @equations, 1 );

    # Fill $matrix and $vector.
    #
    for ($i = 0; $i < @equations; $i++) {
        for ($j = 0; $j < @equations; $j++) {
            assign $matrix ( $i+1, $j+1, $equations[$i][$j] );
        }
        assign $vector ( $i+1, 1, $equations[$i][-1] );
    }

    # Transform $matrix into an LR matrix.
    #
    my $LR = decompose_LR $matrix;

    # Solve the LR matrix for $vector.
    #
    ($dimension, $solution, $base_matrix) = $LR->solve_LR( $vector );

    for ($i = 0; $i < @equations; $i++) {
        $solution[$i] = element $solution( $i+1, 1 );
    }
    return @solution;
}