#!/usr/bin/perl

# jacobian($func_array, $point) calculates the Jacobian matrix at
# $point for the array of functions referred to by $func_array.
# $point is a reference to an array of coordinates.
#
sub jacobian {
    my ($func_array, $point) = @_;
    my ($delta, $i, $j, $k, $coord, @values, @func, @jacobian);
    my $epsilon = 1e-8;

    # Feed the point into each function.
    #
    for ($i = 0; $i < @$func_array; $i++) {
        $values[$i] = &{$func_array->[$i]}( @$point );
    }

    for ($i = 0; $i < @$point; $i++) {
        $coord = $point->[$i];
        $delta = $epsilon * abs($coord) || $epsilon;
        $point->[$i] = $delta + $coord;
        $delta = $point->[$i] - $coord;
        for ($k = 0; $k < @$func_array; $k++) {
            $func[$k] = &{$func_array->[$k]}( @$point );
        }
        $point->[$i] = $coord;
        for ($j = 0; $j < @$func_array; $j++) {
            $jacobian[$j][$i] = ($func[$j] - $values[$j]) / $delta;
        }
    }
    return @jacobian;
}

sub f {
    my ($x, $y, $z) = @_;
    return 3*$x + 2*$x*$y + 4*$y*$z;
}

sub g {
    my ($x, $y, $z) = @_;
    return unless $x;
    return 4/$x + log($x*$y) + $z**$x;
}

sub h {
    my ($x, $y, $z) = @_;
    return $x * $y * $z;
}

@jacobian = jacobian( [\&f, \&g, \&h], [3, 4, -2] );
foreach $row (@jacobian) {
    for ($column = 0; $column < @$row; $column++) {
        print $row->[$column], " ";
    }
    print "\n";
}