sub integrate(&@)
{
    my $f = shift;
    ## Set up some sensical defaults for ranges and vars:
    my ($from, $to, $by, $vars) = do {
        my %o = @_;
        for (qw(from to by vars)) {
            $o{$_} = [$o{$_}] if exists $o{$_} && ! ref $o{$_};
        }
        my $dim = @{$o{from}};
        $o{vars} = [qw(x y z w)[0..$dim-1]] unless $o{vars};
        $o{from} = [(0) x $dim] unless $o{from};
        $o{by} = [(0.01) x $dim] unless $o{by};
        @o{qw(from to by vars)};
    };
    my $vol = 1;
    $vol *= $_ for @$by;
    ## Generate nested evaluation loops:
    local *intgen = sub {
        my ($n, $body) = @_;
        if ($n < 0) {
            '$sum += &$f * ' . $vol;
        } else {
            my ($by, $v) = ($by->[$n], $vars->[$n]);
            my ($lo,$hi) = ($from->[$n] / $by, $to->[$n] / $by);
            "for \$$v ($lo .. $hi) { \$$v *= $by ;\n"
                . intgen($n-1)
                . "\n}\n";
        }
    };

    ## Do it:
    (eval 'sub { my $f = shift; my $sum = 0;'.intgen($#{$from}).' $sum }')
        ->($f);
}

#> integrate { $x * $y } from => [0,0], to => [1,1]
## 0.255025