123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146 |
- #!/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
|