perl-scripts/Math/bernoulli_numbers_from_factorials_visual.pl

#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# Date: 02 December 2017
# https://github.com/trizen

# A new algorithm for computing Bernoulli numbers (visualization).

# Inspired from Norman J. Wildberger video lecture:
#   https://www.youtube.com/watch?v=qmMs6tf8qZ8

# See also:
#   https://en.wikipedia.org/wiki/Bernoulli_number#Connection_with_Pascal’s_triangle

use 5.010;
use strict;
use warnings;

use Math::AnyNum qw(:overload factorial bernfrac);

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

    my @B = (1, (0) x $n);

    foreach my $i (1 .. $n) {

        if ($i % 2 != 0 and $i > 1) {
            ## next;
        }

        foreach my $k (0 .. $i - 1) {

            if ($k % 2 != 0 and $k > 1) {
                ## next;
            }

            my $f = factorial($i - $k + 1);
            my $d = $B[$i] - $B[$k] / $f;

            printf("[%2s, %s] -> %6s / %2s! - %6s / %s!  / %2s! = %6s  / %2s!\n",
                   $i, $k, $B[$i] * factorial($i),
                   $i, $B[$k] * factorial($k),
                   $k,
                   $i - $k + 1,
                   $d * factorial($i), $i);

            $B[$i] = $d;
        }

        say '';
    }

    map { $B[$_] * factorial($_) } 0 .. $#B;
}

my @B = bernoulli_numbers(10);    # first 10 Bernoulli numbers

foreach my $i (0 .. $#B) {

    # Verify the results
    if ($i > 1 and $B[$i] != bernfrac($i)) {
        die "error for i=$i";
    }

    say "B($i) = $B[$i]";
}

__END__
[ 1, 0] ->      0 /  1! -      1 / 0!  /  2! =   -1/2  /  1!

[ 2, 0] ->      0 /  2! -      1 / 0!  /  3! =   -1/3  /  2!
[ 2, 1] ->   -1/3 /  2! -   -1/2 / 1!  /  2! =    1/6  /  2!

[ 3, 0] ->      0 /  3! -      1 / 0!  /  4! =   -1/4  /  3!
[ 3, 1] ->   -1/4 /  3! -   -1/2 / 1!  /  3! =    1/4  /  3!
[ 3, 2] ->    1/4 /  3! -    1/6 / 2!  /  2! =      0  /  3!

[ 4, 0] ->      0 /  4! -      1 / 0!  /  5! =   -1/5  /  4!
[ 4, 1] ->   -1/5 /  4! -   -1/2 / 1!  /  4! =   3/10  /  4!
[ 4, 2] ->   3/10 /  4! -    1/6 / 2!  /  3! =  -1/30  /  4!
[ 4, 3] ->  -1/30 /  4! -      0 / 3!  /  2! =  -1/30  /  4!

[ 5, 0] ->      0 /  5! -      1 / 0!  /  6! =   -1/6  /  5!
[ 5, 1] ->   -1/6 /  5! -   -1/2 / 1!  /  5! =    1/3  /  5!
[ 5, 2] ->    1/3 /  5! -    1/6 / 2!  /  4! =  -1/12  /  5!
[ 5, 3] ->  -1/12 /  5! -      0 / 3!  /  3! =  -1/12  /  5!
[ 5, 4] ->  -1/12 /  5! -  -1/30 / 4!  /  2! =      0  /  5!

[ 6, 0] ->      0 /  6! -      1 / 0!  /  7! =   -1/7  /  6!
[ 6, 1] ->   -1/7 /  6! -   -1/2 / 1!  /  6! =   5/14  /  6!
[ 6, 2] ->   5/14 /  6! -    1/6 / 2!  /  5! =   -1/7  /  6!
[ 6, 3] ->   -1/7 /  6! -      0 / 3!  /  4! =   -1/7  /  6!
[ 6, 4] ->   -1/7 /  6! -  -1/30 / 4!  /  3! =   1/42  /  6!
[ 6, 5] ->   1/42 /  6! -      0 / 5!  /  2! =   1/42  /  6!

[ 7, 0] ->      0 /  7! -      1 / 0!  /  8! =   -1/8  /  7!
[ 7, 1] ->   -1/8 /  7! -   -1/2 / 1!  /  7! =    3/8  /  7!
[ 7, 2] ->    3/8 /  7! -    1/6 / 2!  /  6! =  -5/24  /  7!
[ 7, 3] ->  -5/24 /  7! -      0 / 3!  /  5! =  -5/24  /  7!
[ 7, 4] ->  -5/24 /  7! -  -1/30 / 4!  /  4! =   1/12  /  7!
[ 7, 5] ->   1/12 /  7! -      0 / 5!  /  3! =   1/12  /  7!
[ 7, 6] ->   1/12 /  7! -   1/42 / 6!  /  2! =      0  /  7!

[ 8, 0] ->      0 /  8! -      1 / 0!  /  9! =   -1/9  /  8!
[ 8, 1] ->   -1/9 /  8! -   -1/2 / 1!  /  8! =   7/18  /  8!
[ 8, 2] ->   7/18 /  8! -    1/6 / 2!  /  7! =  -5/18  /  8!
[ 8, 3] ->  -5/18 /  8! -      0 / 3!  /  6! =  -5/18  /  8!
[ 8, 4] ->  -5/18 /  8! -  -1/30 / 4!  /  5! =  17/90  /  8!
[ 8, 5] ->  17/90 /  8! -      0 / 5!  /  4! =  17/90  /  8!
[ 8, 6] ->  17/90 /  8! -   1/42 / 6!  /  3! =  -1/30  /  8!
[ 8, 7] ->  -1/30 /  8! -      0 / 7!  /  2! =  -1/30  /  8!

[ 9, 0] ->      0 /  9! -      1 / 0!  / 10! =  -1/10  /  9!
[ 9, 1] ->  -1/10 /  9! -   -1/2 / 1!  /  9! =    2/5  /  9!
[ 9, 2] ->    2/5 /  9! -    1/6 / 2!  /  8! =  -7/20  /  9!
[ 9, 3] ->  -7/20 /  9! -      0 / 3!  /  7! =  -7/20  /  9!
[ 9, 4] ->  -7/20 /  9! -  -1/30 / 4!  /  6! =   7/20  /  9!
[ 9, 5] ->   7/20 /  9! -      0 / 5!  /  5! =   7/20  /  9!
[ 9, 6] ->   7/20 /  9! -   1/42 / 6!  /  4! =  -3/20  /  9!
[ 9, 7] ->  -3/20 /  9! -      0 / 7!  /  3! =  -3/20  /  9!
[ 9, 8] ->  -3/20 /  9! -  -1/30 / 8!  /  2! =      0  /  9!

[10, 0] ->      0 / 10! -      1 / 0!  / 11! =  -1/11  / 10!
[10, 1] ->  -1/11 / 10! -   -1/2 / 1!  / 10! =   9/22  / 10!
[10, 2] ->   9/22 / 10! -    1/6 / 2!  /  9! = -14/33  / 10!
[10, 3] -> -14/33 / 10! -      0 / 3!  /  8! = -14/33  / 10!
[10, 4] -> -14/33 / 10! -  -1/30 / 4!  /  7! =  19/33  / 10!
[10, 5] ->  19/33 / 10! -      0 / 5!  /  6! =  19/33  / 10!
[10, 6] ->  19/33 / 10! -   1/42 / 6!  /  5! = -14/33  / 10!
[10, 7] -> -14/33 / 10! -      0 / 7!  /  4! = -14/33  / 10!
[10, 8] -> -14/33 / 10! -  -1/30 / 8!  /  3! =   5/66  / 10!
[10, 9] ->   5/66 / 10! -      0 / 9!  /  2! =   5/66  / 10!

B(0) = 1
B(1) = -1/2
B(2) = 1/6
B(3) = 0
B(4) = -1/30
B(5) = 0
B(6) = 1/42
B(7) = 0
B(8) = -1/30
B(9) = 0
B(10) = 5/66