lerch_zeta_function.pl 1.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152
  1. #!/usr/bin/perl
  2. # Daniel "Trizen" Șuteu
  3. # Date: 03 January 2018
  4. # https://github.com/trizen
  5. # Simple implementation of the Lerch zeta function Φ(z, s, t), for real(z) < 1/2.
  6. # Formula due to Guillera and Sondow (2005).
  7. # See also:
  8. # https://mathworld.wolfram.com/LerchTranscendent.html
  9. # https://en.wikipedia.org/wiki/Lerch_zeta_function
  10. use 5.020;
  11. use warnings;
  12. use experimental qw(signatures);
  13. use Math::AnyNum qw(:overload pi binomial factorial);
  14. sub lerch ($z, $s, $t, $reps = 100) {
  15. my $sum = 0.0;
  16. my $r = (-$z) / (1 - $z);
  17. foreach my $n (0 .. $reps) {
  18. my $temp = 0.0;
  19. foreach my $k (0 .. $n) {
  20. $temp += (-1)**$k * binomial($n, $k) * ($t + $k)**(-$s);
  21. }
  22. $sum += $r**$n * $temp;
  23. }
  24. $sum / (1 - $z);
  25. }
  26. say "zeta(2)/2 =~ ", lerch(-1, 2, 1); # 0.822467033424113...
  27. say "4*catalan =~ ", lerch(-1, 2, 1 / 2); # 3.663862376708876...
  28. say '';
  29. sub A281964 ($n) {
  30. (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;
  31. }
  32. foreach my $n (1 .. 10) {
  33. printf("a(%2d) = %s\n", $n, A281964($n));
  34. }