bernoulli_numbers_formulas.pl

#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# License: GPLv3
# Date: 21 September 2015
# Website: https://github.com/trizen

use 5.014;
use strict;
use warnings;

use Math::AnyNum qw(:overload);

# Translation of:
#   https://en.wikipedia.org/wiki/Bernoulli_number#Algorithmic_description

sub bernoulli_number {
    my ($n) = @_;

    # return 0 if $n > 1 && $n % 2;    # Bn = 0 for all odd n > 1

    my @A;
    for my $m (0 .. $n) {
        $A[$m] = 1 / ($m + 1);

        for (my $j = $m ; $j > 0 ; $j--) {
            $A[$j - 1] = ((($j == 1 ? '' : "$j*") . '(' . join('-', ($A[$j - 1], $A[$j])) . ')') =~ s/^\((.*?)\)\z/$1/r);
        }
    }

    return $A[0];    # which is Bn
}

foreach my $i (0 .. 6) {
    printf("B(%d) = %s\n", $i, bernoulli_number($i));
}

__END__
B(0) = 1
B(1) = 1-1/2
B(2) = 1-1/2-2*(1/2-1/3)
B(3) = 1-1/2-2*(1/2-1/3)-2*(2*(1/2-1/3)-3*(1/3-1/4))
B(4) = 1-1/2-2*(1/2-1/3)-2*(2*(1/2-1/3)-3*(1/3-1/4))-2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))
B(5) = 1-1/2-2*(1/2-1/3)-2*(2*(1/2-1/3)-3*(1/3-1/4))-2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))-2*(2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))-3*(3*(3*(1/3-1/4)-4*(1/4-1/5))-4*(4*(1/4-1/5)-5*(1/5-1/6))))
B(6) = 1-1/2-2*(1/2-1/3)-2*(2*(1/2-1/3)-3*(1/3-1/4))-2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))-2*(2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))-3*(3*(3*(1/3-1/4)-4*(1/4-1/5))-4*(4*(1/4-1/5)-5*(1/5-1/6))))-2*(2*(2*(2*(2*(1/2-1/3)-3*(1/3-1/4))-3*(3*(1/3-1/4)-4*(1/4-1/5)))-3*(3*(3*(1/3-1/4)-4*(1/4-1/5))-4*(4*(1/4-1/5)-5*(1/5-1/6))))-3*(3*(3*(3*(1/3-1/4)-4*(1/4-1/5))-4*(4*(1/4-1/5)-5*(1/5-1/6)))-4*(4*(4*(1/4-1/5)-5*(1/5-1/6))-5*(5*(1/5-1/6)-6*(1/6-1/7)))))