-
Notifications
You must be signed in to change notification settings - Fork 33
/
bernoulli_numbers_from_factorials_mpz.pl
executable file
·96 lines (68 loc) · 2.36 KB
/
bernoulli_numbers_from_factorials_mpz.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#!/usr/bin/perl
# Author: Daniel "Trizen" Șuteu
# Date: 07 July 2018
# https://github.com/trizen
# A new algorithm for computing Bernoulli numbers.
# 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::GMPq;
use Math::GMPz;
sub bernoulli_numbers {
my ($n) = @_;
my @A = (Math::GMPz::Rmpz_init_set_ui(1));
my @B = (Math::GMPz::Rmpz_init_set_ui(1));
my @F = (Math::GMPz::Rmpz_init_set_ui(1));
foreach my $k (1 .. $n) {
$F[$k] = Math::GMPz::Rmpz_init();
$A[$k] = Math::GMPz::Rmpz_init_set_ui(0);
$B[$k] = Math::GMPz::Rmpz_init_set_ui(1);
Math::GMPz::Rmpz_mul_ui($F[$k], $F[$k - 1], $k);
}
Math::GMPz::Rmpz_mul_ui($F[$n + 1] = Math::GMPz::Rmpz_init(), $F[$n], $n + 1);
my $t = Math::GMPz::Rmpz_init();
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 $r = $i - $k + 1;
Math::GMPz::Rmpz_mul($A[$i], $A[$i], $F[$r]);
Math::GMPz::Rmpz_mul($A[$i], $A[$i], $B[$k]);
Math::GMPz::Rmpz_submul($A[$i], $B[$i], $A[$k]);
Math::GMPz::Rmpz_mul($B[$i], $B[$i], $F[$r]);
Math::GMPz::Rmpz_mul($B[$i], $B[$i], $B[$k]);
Math::GMPz::Rmpz_gcd($t, $A[$i], $B[$i]);
Math::GMPz::Rmpz_divexact($A[$i], $A[$i], $t);
Math::GMPz::Rmpz_divexact($B[$i], $B[$i], $t);
}
}
my @R = @A;
for (my $k = 2 ; $k <= $#B ; $k += 2) {
Math::GMPz::Rmpz_mul($A[$k], $A[$k], $F[$k]);
my $bern = Math::GMPq::Rmpq_init();
Math::GMPq::Rmpq_set_num($bern, $A[$k]);
Math::GMPq::Rmpq_set_den($bern, $B[$k]);
Math::GMPq::Rmpq_canonicalize($bern);
$R[$k] = $bern;
}
if ($#R > 0) {
my $bern = Math::GMPq::Rmpq_init();
Math::GMPq::Rmpq_set_num($bern, $A[1]);
Math::GMPq::Rmpq_set_den($bern, $B[1]);
Math::GMPq::Rmpq_canonicalize($bern);
$R[1] = $bern;
}
return @R;
}
my @B = bernoulli_numbers(100); # first 100 Bernoulli numbers
foreach my $i (0 .. $#B) {
say "B($i) = $B[$i]";
}