#!/usr/bin/perl

use constant pi => 3.14159265358979;

# bernoulli($n) returns the $nth Bernoulli number when called
#   in a scalar context, and all of the Bernoulli numbers up
#   to $n when called in a list context.
#
sub bernoulli {
    my ($n) = shift;

    # Take care of the easy cases.
    return 1    if $n == 0;
    return -0.5 if $n == 1;
    return 0    if $n % 2;

    my (@bernoulli) = (1, -0.5);
    my ($i, $j);

    # Stop when our array of Bernoulli numbers has
    # grown to have $n members.
    #
    while ($#bernoulli < $n) {

        # If we're looking at an odd index, fill in zero and move on.
        $#bernoulli % 2 || (push(@bernoulli, 0), next);

        # Otherwise, compute the next Bernoulli number,
        # using (as we must) the numbers that we just computed.
        # choose() is defined in Probability, (Chapter 14).
        #
        for ($i = 0, $j = 0; $i <= $#bernoulli; $i++) {
            $j += choose($#bernoulli + 2, $i) * $bernoulli[$i];
        }
        $j /= - ($#bernoulli + 2);
        push @bernoulli, $j;
    }

    # Return all the numbers if called in list context,
    # or the last number if called in scalar context.
    #
    return wantarray ? @bernoulli : $bernoulli[$#bernoulli];
}

@bernoullis = bernoulli(4);   # list context

print "Bernoulli numbers: @bernoullis\n";

# zeta_iterative($n, $terms) computes the first $terms terms
#   of Riemann's zeta function.  This will work for either even
#   or odd $n, but should only be used for odd $n since a
#   closed-form solution exists for even $n.
#
sub zeta_iterative {
    my ($n, $terms) = @_;
    my ($i, $result);
    for ($i = 1; $i <= $terms; $i++) {
        $result += 1 / ($i ** $n);
    }
    return $result;
}

sub zeta {
    my ($n) = shift;
    if ($n % 2) {                         # if $n is odd
        return zeta_iterative($n, 10000);
    } else {
        return .5 * abs(bernoulli($n)) * ((2*pi) ** $n) / factorial($n);
    }
}

foreach $n (2, 4, 6, 8) {
    print zeta($n), "\n";
}

# choose($n, $k) is the number of ways to choose $k elements from a set
# of $n elements, when the order of selection is irrelevant.
#
sub choose {
    my ($n, $k) = @_;
    my ($result, $j) = (1, 1);

    return 0 if $k > $n || $k < 0;
    $k = ($n - $k) if ($n - $k) < $k;

    while ( $j <= $k ) {
        $result *= $n--;
        $result /= $j++;
    }
    return $result;
}

sub factorial {
    my ($n, $res) = (shift, 1);

    # Nonintegers require the gamma function,
    # discussed later in the chapter.
    return undef unless $n >= 0 and $n == int($n);

    $res *= $n-- while $n > 1;
    return $res;
}