-
Notifications
You must be signed in to change notification settings - Fork 34
/
zeta_2n.pl
executable file
·48 lines (35 loc) · 979 Bytes
/
zeta_2n.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
#!/usr/bin/perl
# Author: Daniel "Trizen" Șuteu
# License: GPLv3
# Date: 06 September 2015
# Website: https://github.com/trizen
# Calculate zeta(2n) using a closed-form formula.
# See: https://en.wikipedia.org/wiki/Riemann_zeta_function
use 5.010;
use strict;
use warnings;
use Memoize qw(memoize);
use Math::AnyNum qw(:overload pi);
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 * ($A[$j - 1] - $A[$j]);
}
}
return $A[0]; # which is Bn
}
sub factorial {
$_[0] < 2 ? 1 : factorial($_[0] - 1) * $_[0];
}
memoize('factorial');
sub zeta_2n {
my ($n2) = 2 * $_[0];
((-1)**($_[0] + 1) * 2**($n2 - 1) * pi**$n2 * bernoulli_number($n2)) / factorial($n2);
}
for my $i (1 .. 10) {
say "zeta(", 2 * $i, ") = ", zeta_2n($i);
}