12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152 |
- #!/usr/bin/perl
- # A very high-level computation of the nth-Bernoulli number, using prime numbers.
- # Algorithm due to Kevin J. McGown (December 8, 2005)
- # See his paper: "Computing Bernoulli Numbers Quickly"
- use 5.010;
- use strict;
- use warnings;
- use lib qw(../lib);
- #use Math::AnyNum qw(:overload); # can be uncommented
- use Math::AnyNum qw(factorial next_prime ceil float);
- sub bernoulli_from_primes {
- my ($n) = @_;
- $n == 0 and return Math::AnyNum->one;
- $n == 1 and return Math::AnyNum->new('1/2');
- $n < 0 and return Math::AnyNum->nan;
- $n % 2 and return Math::AnyNum->zero;
- my $tau = 6.28318530717958647692528676655900576839433879875;
- my $log2B = (log(4 * $tau * $n) / 2 + $n * log($n) - $n * log($tau) - $n) / log(2);
- local $Math::AnyNum::PREC = int($n + $log2B) + ($n <= 90 ? 18 : 0);
- my $K = factorial($n) * 2 / Math::AnyNum->tau**$n;
- my $d = 1;
- for (my $p = 2 ; $p <= $n + 1 ; $p = next_prime($p)) {
- if ($n % ($p - 1) == 0) {
- $d *= $p;
- }
- }
- my $N = ceil(($K * $d)->root($n - 1));
- my $z = 1.0;
- for (my $p = 2 ; $p <= $N ; $p = next_prime($p)) {
- my $u = float($p)**$n;
- $z *= $u / ($u-1);
- }
- (-1)**($n / 2 + 1) * int(ceil($d * $K * $z)) / $d;
- }
- foreach my $n (0 .. 100) {
- printf "B%-3d = %s\n", 2 * $n, bernoulli_from_primes(2 * $n);
- }
|