#!/usr/bin/perl -l

# integrate() uses the Romberg algorithm to estimate the definite integral
# of the function $func (provided as a code reference) from $lo to $hi.
#
# The subroutine will compute roughly ($steps + 1) * ($steps + 2) / 2
# estimates for the integral, of which the last will be the most accurate.
#
# integrate() returns early if intermediate estimates change by less
# than $epsilon.
#
sub integrate {
    my ($func, $lo, $hi, $steps, $epsilon) = @_;
    my ($h) = $hi - $lo;
    my ($i, $j, @r, $sum);

    # Our initial estimate.
    $est[0][0] = ($h / 2) * ( &{$func}( $lo ) + &{$func}( $hi ) );

    # Compute each row of the Romberg array.
    for ($i = 1; $i <= $steps; $i++) {

        $h /= 2;
        $sum = 0;

        # Compute the first column of the current row.
        for ($j = 1; $j < 2 ** $i; $j += 2) {
            $sum += &{$func}( $lo + $j * $h );
        }
        $est[$i][0] = $est[$i-1][0] / 2 + $sum * $h;

        # Compute the rest of the columns in this row.
        for ($j = 1; $j <= $i; $j++) {
            $est[$i][$j] = ($est[$i][$j-1] - $est[$i-1][$j-1])
                / (4**$j - 1) + $est[$i][$j-1];
        }

        # Are we close enough?
        return $est[$i][$i] if $epsilon and
            abs($est[$i][$i] - $est[$i-1][$i-1]) <= $epsilon;
    }
    return $est[$steps][$steps];
}

# bandwidth() returns the bandwidth between Boston and Helsinki
# for a given hour, in kilobytes per second.
#
sub bandwidth {
    my $time = shift;
    $time = $time % 24 + $time - int($time);
    # Business hours in Boston or Helsinki
    if ($time >= 3 && $time < 17) {
        return 51 - 50 / (($time - 10) ** 2 + 1);
    } else {
        $time = 20 - $time if $time < 3;
        return 200 - 6 * (($time - 22) ** 2);
    }
}

use constant epsilon  => 1e-14;
use constant pi       => 3.14159265358979;
use constant infinity => 1000;

$five_to_eleven_am = integrate( \&bandwidth,  5, 11, 6, epsilon ) * 3.6e6;
print $five_to_eleven_am;                                          # 713 meg

$five_to_eleven_pm = integrate( \&bandwidth, 17, 23, 6, epsilon ) * 3.6e6;
print $five_to_eleven_pm;                                          # 3.4 gig

$nine_to_eleven_am = integrate( \&bandwidth,  9, 11, 6, epsilon )  * 3.6e6;
print $nine_to_eleven_am;                                          # 84 meg

sub area   { 2 * pi / $_[0] }    # Surface area of a slice of Gabriel's Horn
sub volume { pi / ($_[0] ** 2) } # Volume of a slice of Gabriel's Horn

$gabriel_area   = integrate(\&area,   1, infinity, 10, epsilon);
$gabriel_volume = integrate(\&volume, 1, infinity, 10, epsilon);

print "Volume is $gabriel_volume, but area is $gabriel_area.\n";