zeta_2n.pl 979 B

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849
  1. #!/usr/bin/perl
  2. # Author: Daniel "Trizen" Șuteu
  3. # License: GPLv3
  4. # Date: 06 September 2015
  5. # Website: https://github.com/trizen
  6. # Calculate zeta(2n) using a closed-form formula.
  7. # See: https://en.wikipedia.org/wiki/Riemann_zeta_function
  8. use 5.010;
  9. use strict;
  10. use warnings;
  11. use Memoize qw(memoize);
  12. use Math::AnyNum qw(:overload pi);
  13. sub bernoulli_number {
  14. my ($n) = @_;
  15. return 0 if $n > 1 && $n % 2; # Bn = 0 for all odd n > 1
  16. my @A;
  17. for my $m (0 .. $n) {
  18. $A[$m] = 1 / ($m + 1);
  19. for (my $j = $m ; $j > 0 ; $j--) {
  20. $A[$j - 1] = $j * ($A[$j - 1] - $A[$j]);
  21. }
  22. }
  23. return $A[0]; # which is Bn
  24. }
  25. sub factorial {
  26. $_[0] < 2 ? 1 : factorial($_[0] - 1) * $_[0];
  27. }
  28. memoize('factorial');
  29. sub zeta_2n {
  30. my ($n2) = 2 * $_[0];
  31. ((-1)**($_[0] + 1) * 2**($n2 - 1) * pi**$n2 * bernoulli_number($n2)) / factorial($n2);
  32. }
  33. for my $i (1 .. 10) {
  34. say "zeta(", 2 * $i, ") = ", zeta_2n($i);
  35. }