-
Notifications
You must be signed in to change notification settings - Fork 34
/
lerch_zeta_function.pl
executable file
·51 lines (34 loc) · 1.14 KB
/
lerch_zeta_function.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
#!/usr/bin/perl
# Daniel "Trizen" Șuteu
# Date: 03 January 2018
# https://github.com/trizen
# Simple implementation of the Lerch zeta function Φ(z, s, t), for real(z) < 1/2.
# Formula due to Guillera and Sondow (2005).
# See also:
# https://mathworld.wolfram.com/LerchTranscendent.html
# https://en.wikipedia.org/wiki/Lerch_zeta_function
use 5.020;
use warnings;
use experimental qw(signatures);
use Math::AnyNum qw(:overload pi binomial factorial);
sub lerch ($z, $s, $t, $reps = 100) {
my $sum = 0.0;
my $r = (-$z) / (1 - $z);
foreach my $n (0 .. $reps) {
my $temp = 0.0;
foreach my $k (0 .. $n) {
$temp += (-1)**$k * binomial($n, $k) * ($t + $k)**(-$s);
}
$sum += $r**$n * $temp;
}
$sum / (1 - $z);
}
say "zeta(2)/2 =~ ", lerch(-1, 2, 1); # 0.822467033424113...
say "4*catalan =~ ", lerch(-1, 2, 1 / 2); # 3.663862376708876...
say '';
sub A281964 ($n) {
(factorial($n) * (-2 * i * i**$n * (lerch(-1, 1, $n / 2 + 1) - i * lerch(-1, 1, ($n + 1) / 2)) + pi + 2 * i * log(2)) / 4)->real->round;
}
foreach my $n (1 .. 10) {
printf("a(%2d) = %s\n", $n, A281964($n));
}