Bernoulli numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Bernoulli numbers are used in some series expansions of several functions (trigonometric, hyperbolic, gamma, etc.), and are extremely important in number theory and analysis. Note that there are two definitions of Bernoulli numbers; this task will be using the modern usage (as per the National Institute of Standards and Technology convention). The 'th Bernoulli number is expressed as Bn.
- Task
- show the Bernoulli numbers B0 through B60, suppressing all output of values which are equal to zero. (Other than B1, all odd Bernoulli numbers have a value of 0 (zero).
- express the numbers as fractions (most are improper fractions).
- fractions should be reduced.
- index each number in some way so that it can be discerned which number is being displayed.
- align the solidi (/) if used (extra credit).
- An algorithm
The numbers can be computed with the following algorithm (from wikipedia)
for m from 0 by 1 to n do
A[m] ← 1/(m+1)
for j from m by -1 to 1 do
A[j-1] ← j×(A[j-1] - A[j])
return A[0] (which is Bn)
- See also
- Sequence A027641 Numerator of Bernoulli number B_n on The On-Line Encyclopedia of Integer Sequences.
- Sequence A027642 Denominator of Bernoulli number B_n on The On-Line Encyclopedia of Integer Sequences.
- Entry Bernoulli number on The Eric Weisstein's World of Mathematics (TM).
- Luschny's The Bernoulli Manifesto for a discussion on = -½ vs. +½.
Common Lisp
An implementation of the simple algorithm.
Be advised that the pseudocode algorithm specifies (j * (a[j-1] - a[j])) in the inner loop; implementing that as-is gives the wrong value (1/2) where n = 1, whereas subtracting a[j]-a[j-1] yields the correct value (B[1]=-1/2). See the numerator list.
<lang lisp>(defun bernouilli (n)
(loop with a = (make-array (list (1+ n)))
for m from 0 to n do
(setf (aref a m) (/ 1 (+ m 1)))
(loop for j from m downto 1 do
(setf (aref a (- j 1))
(* j (- (aref a j) (aref a (- j 1))))))
finally (return (aref a 0))))
- Print outputs to stdout
(loop for n from 0 to 60 do
(let ((b (bernouilli n)))
(when (not (zerop b))
(format t "~a: ~a~%" n b))))
- For the "extra credit" challenge, we need to align the slashes.
(let (results)
;;collect the results
(loop for n from 0 to 60 do
(let ((b (bernouilli n)))
(when (not (zerop b)) (push (cons b n) results))))
;;parse the numerators into strings; save the greatest length in max-length
(let ((max-length (apply #'max (mapcar (lambda (r)
(length (format nil "~a" (numerator r))))
(mapcar #'car results)))))
;;Print the numbers with using the fixed-width formatter: ~Nd, where N is
;;the number of leading spaces. We can't just pass in the width variable
;;but we can splice together a formatting string that includes it.
;;We also can't use the fixed-width formatter on a ratio, so we have to split
;;the ratio and splice it back together like idiots.
(loop for n in (mapcar #'cdr (reverse results))
for r in (mapcar #'car (reverse results)) do
(format t (concatenate 'string
"B(~2d): ~"
(format nil "~a" max-length)
"d/~a~%")
n
(numerator r)
(denominator r)))))</lang>
- Output:
B( 0): 1/1 B( 1): -1/2 B( 2): 1/6 B( 4): -1/30 B( 6): 1/42 B( 8): -1/30 B(10): 5/66 B(12): -691/2730 B(14): 7/6 B(16): -3617/510 B(18): 43867/798 B(20): -174611/330 B(22): 854513/138 B(24): -236364091/2730 B(26): 8553103/6 B(28): -23749461029/870 B(30): 8615841276005/14322 B(32): -7709321041217/510 B(34): 2577687858367/6 B(36): -26315271553053477373/1919190 B(38): 2929993913841559/6 B(40): -261082718496449122051/13530 B(42): 1520097643918070802691/1806 B(44): -27833269579301024235023/690 B(46): 596451111593912163277961/282 B(48): -5609403368997817686249127547/46410 B(50): 495057205241079648212477525/66 B(52): -801165718135489957347924991853/1590 B(54): 29149963634884862421418123812691/798 B(56): -2479392929313226753685415739663229/870 B(58): 84483613348880041862046775994036021/354 B(60): -1215233140483755572040304994079820246041491/56786730
D
This uses the D module from the Arithmetic/Rational task.
<lang d>import std.stdio, std.range, std.algorithm, std.typecons, std.conv,
arithmetic_rational;
auto bernoulli(in uint n) pure nothrow {
auto A = new Rational[n + 1];
foreach (immutable m; 0 .. n + 1) {
A[m] = Rational(1, m + 1);
foreach_reverse (immutable j; 1 .. m + 1)
A[j - 1] = j * (A[j - 1] - A[j]);
}
return A[0];
}
void main() {
immutable bn = 61.iota.map!bernoulli.enumerate.filter!(t => t[1]).array;
immutable width = bn.map!(b => b[1].numerator.text.length).reduce!max;
foreach (immutable b; bn)
writefln("B(%2d) = %*d/%d",
b[0], width, b[1].numerator, b[1].denominator);
}</lang> The output is exactly the same as the Python entry.
FunL
FunL has pre-defined function B in module integers, which is defined as:
<lang funl>import integers.choose
def B( n ) = sum( 1/(k + 1)*sum((if 2|r then 1 else -1)*choose(k, r)*(r^n) | r <- 0..k) | k <- 0..n )
for i <- 0..60 if i == 1 or 2|i
printf( "B(%2d) = %s\n", i, B(i) )</lang>
- Output:
B( 0) = 1 B( 1) = -1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
Go
<lang go>package main
import (
"fmt" "math/big" "strings"
)
func b(n int) *big.Rat {
var f big.Rat
a := make([]big.Rat, n+1)
for m := range a {
a[m].SetFrac64(1, int64(m+1))
for j := m; j >= 1; j-- {
d := &a[j-1]
d.Mul(f.SetInt64(int64(j)), d.Sub(d, &a[j]))
}
}
return f.Set(&a[0])
}
func align(b *big.Rat, w int) string {
s := b.String()
return strings.Repeat(" ", w-strings.Index(s, "/")) + s
}
func main() {
for n := 0; n <= 60; n++ {
if b := b(n); b.Num().BitLen() > 0 {
fmt.Printf("B(%2d) =%s\n", n, align(b, 45))
}
}
}</lang>
- Output:
B( 0) = 1/1 B( 1) = 1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
Haskell
This program works as a command line utility, that reads from stdin the number of elements to compute (default 60) and prints them in stdout. The implementation of the algorithm is in the function bernoullis. The rest is for printing the results.
<lang Haskell>module Main where
import Data.Ratio import System.Environment
main = getArgs >>= printM . defaultArg where
defaultArg as = if null as then 60 else read (head as)
printM m = mapM_ (putStrLn . printP) . takeWhile ((<= m).fst)
. filter (\(_,b) -> b /= 0%1) . zip [0..] $ bernoullis
printP (i,r) = "B(" ++ show i ++ ")=" ++ show (numerator r) ++ "/" ++ show (denominator r)
bernoullis = map head . iterate (ulli 1) . map berno $ enumFrom 0 where
berno i = 1 % (i+1) ulli _ [_] = [] ulli i (x:y:xs) = (i%1)*(x-y) : ulli (i+1) (y:xs)
</lang>
- Output:
B(0)=1/1 B(1)=1/2 B(2)=1/6 B(4)=-1/30 B(6)=1/42 B(8)=-1/30 B(10)=5/66 B(12)=-691/2730 B(14)=7/6 B(16)=-3617/510 B(18)=43867/798 B(20)=-174611/330 B(22)=854513/138 B(24)=-236364091/2730 B(26)=8553103/6 B(28)=-23749461029/870 B(30)=8615841276005/14322 B(32)=-7709321041217/510 B(34)=2577687858367/6 B(36)=-26315271553053477373/1919190 B(38)=2929993913841559/6 B(40)=-261082718496449122051/13530 B(42)=1520097643918070802691/1806 B(44)=-27833269579301024235023/690 B(46)=596451111593912163277961/282 B(48)=-5609403368997817686249127547/46410 B(50)=495057205241079648212477525/66 B(52)=-801165718135489957347924991853/1590 B(54)=29149963634884862421418123812691/798 B(56)=-2479392929313226753685415739663229/870 B(58)=84483613348880041862046775994036021/354 B(60)=-1215233140483755572040304994079820246041491/56786730
Icon and Unicon
The following works in both languages: <lang unicon>link "rational"
procedure main(args)
limit := integer(!args) | 60
every b := bernoulli(i := 0 to limit) do
if b.numer > 0 then write(right(i,3),": ",align(rat2str(b),60))
end
procedure bernoulli(n)
(A := table(0))[0] := rational(1,1,1)
every m := 1 to n do {
A[m] := rational(1,m+1,1)
every j := m to 1 by -1 do A[j-1] := mpyrat(rational(j,1,1), subrat(A[j-1],A[j]))
}
return A[0]
end
procedure align(r,n)
return repl(" ",n-find("/",r))||r
end</lang>
Sample run:
->bernoulli 60 0: (1/1) 1: (1/2) 2: (1/6) 4: (-1/30) 6: (1/42) 8: (-1/30) 10: (5/66) 12: (-691/2730) 14: (7/6) 16: (-3617/510) 18: (43867/798) 20: (-174611/330) 22: (854513/138) 24: (-236364091/2730) 26: (8553103/6) 28: (-23749461029/870) 30: (8615841276005/14322) 32: (-7709321041217/510) 34: (2577687858367/6) 36: (-26315271553053477373/1919190) 38: (2929993913841559/6) 40: (-261082718496449122051/13530) 42: (1520097643918070802691/1806) 44: (-27833269579301024235023/690) 46: (596451111593912163277961/282) 48: (-5609403368997817686249127547/46410) 50: (495057205241079648212477525/66) 52: (-801165718135489957347924991853/1590) 54: (29149963634884862421418123812691/798) 56: (-2479392929313226753685415739663229/870) 58: (84483613348880041862046775994036021/354) 60: (-1215233140483755572040304994079820246041491/56786730) ->
J
Implementation:
<lang J>B=:3 :0"0
+/,(<:*(_1^[)*!*(y^~1+[)%1+])"0/~i.1x+y
)</lang>
Task:
<lang J> require'strings'
'B',.rplc&'r/_-'"1": (#~ 0 ~: {:"1)(,.B) i.61
B 0 1 B 1 1/2 B 2 1/6 B 4 -1/30 B 6 1/42 B 8 -1/30 B10 5/66 B12 -691/2730 B14 7/6 B16 -3617/510 B18 43867/798 B20 -174611/330 B22 854513/138 B24 -236364091/2730 B26 8553103/6 B28 -23749461029/870 B30 8615841276005/14322 B32 -7709321041217/510 B34 2577687858367/6 B36 -26315271553053477373/1919190 B38 2929993913841559/6 B40 -261082718496449122051/13530 B42 1520097643918070802691/1806 B44 -27833269579301024235023/690 B46 596451111593912163277961/282 B48 -5609403368997817686249127547/46410 B50 495057205241079648212477525/66 B52 -801165718135489957347924991853/1590 B54 29149963634884862421418123812691/798 B56 -2479392929313226753685415739663229/870 B58 84483613348880041862046775994036021/354 B60 -1215233140483755572040304994079820246041491/56786730</lang>
Maple
<lang Maple>k:=length(numer(bernoulli(60))):
G := proc(n) local b,i;
b := `if`(n=1/2,1/2,bernoulli(2*n));
printf("%a%s%a\n",'B'[2*n],
cat(" "$i=1..(5+k-length(numer(b))
+(1+signum(b))/2)-length(2*n)),b);
NULL;
end proc:
G(0), G(1/2), seq(G(i),i=1..30);</lang>
Mathematica
Mathematica has no native way for starting an array at index 0. I therefore had to build the array from 1 to n+1 instead of from 0 to n, adjusting the formula accordingly. <lang Mathematica>bernoulli[n_] := Module[{a = ConstantArray[0, n + 2]},
Do[ am = 1/m; If[m == 1 && a1 != 0, Print[{m - 1, a1}]]; Do[ aj - 1 = (j - 1)*(aj - 1 - aj); If[j == 2 && a1 != 0, Print[{m - 1, a1}]]; , {j, m, 2, -1}]; , {m, 1, n + 1}]; ]
bernoulli[60]</lang>
- Output:
{0,1}
{1,1/2}
{2,1/6}
{4,-(1/30)}
{6,1/42}
{8,-(1/30)}
{10,5/66}
{12,-(691/2730)}
{14,7/6}
{16,-(3617/510)}
{18,43867/798}
{20,-(174611/330)}
{22,854513/138}
{24,-(236364091/2730)}
{26,8553103/6}
{28,-(23749461029/870)}
{30,8615841276005/14322}
{32,-(7709321041217/510)}
{34,2577687858367/6}
{36,-(26315271553053477373/1919190)}
{38,2929993913841559/6}
{40,-(261082718496449122051/13530)}
{42,1520097643918070802691/1806}
{44,-(27833269579301024235023/690)}
{46,596451111593912163277961/282}
{48,-(5609403368997817686249127547/46410)}
{50,495057205241079648212477525/66}
{52,-(801165718135489957347924991853/1590)}
{54,29149963634884862421418123812691/798}
{56,-(2479392929313226753685415739663229/870)}
{58,84483613348880041862046775994036021/354}
{60,-(1215233140483755572040304994079820246041491/56786730)}
Or, it's permissible to use the native Bernoulli number function instead of being forced to use the specified algorithm, we very simply have:
(Note from task's author: nobody is forced to use any specific algorithm, the one shown is just a suggestion.)
<lang Mathematica>Table[{i, BernoulliB[i]}, {i, 0, 60}]; Select[%, #2 != 0 &] // TableForm</lang>
- Output:
0 1 1 -(1/2) 2 1/6 4 -(1/30) 6 1/42 8 -(1/30) 10 5/66 12 -(691/2730) 14 7/6 16 -(3617/510) 18 43867/798 20 -(174611/330) 22 854513/138 24 -(236364091/2730) 26 8553103/6 28 -(23749461029/870) 30 8615841276005/14322 32 -(7709321041217/510) 34 2577687858367/6 36 -(26315271553053477373/1919190) 38 2929993913841559/6 40 -(261082718496449122051/13530) 42 1520097643918070802691/1806 44 -(27833269579301024235023/690) 46 596451111593912163277961/282 48 -(5609403368997817686249127547/46410) 50 495057205241079648212477525/66 52 -(801165718135489957347924991853/1590) 54 29149963634884862421418123812691/798 56 -(2479392929313226753685415739663229/870) 58 84483613348880041862046775994036021/354 60 -(1215233140483755572040304994079820246041491/56786730)
PARI/GP
<lang parigp>for(n=0,60,t=bernfrac(n);if(t,print(n" "t)))</lang>
- Output:
0 1 1 -1/2 2 1/6 4 -1/30 6 1/42 8 -1/30 10 5/66 12 -691/2730 14 7/6 16 -3617/510 18 43867/798 20 -174611/330 22 854513/138 24 -236364091/2730 26 8553103/6 28 -23749461029/870 30 8615841276005/14322 32 -7709321041217/510 34 2577687858367/6 36 -26315271553053477373/1919190 38 2929993913841559/6 40 -261082718496449122051/13530 42 1520097643918070802691/1806 44 -27833269579301024235023/690 46 596451111593912163277961/282 48 -5609403368997817686249127547/46410 50 495057205241079648212477525/66 52 -801165718135489957347924991853/1590 54 29149963634884862421418123812691/798 56 -2479392929313226753685415739663229/870 58 84483613348880041862046775994036021/354 60 -1215233140483755572040304994079820246041491/56786730
Perl
The only thing in the suggested algorithm which depends on N is the number of times through the inner block. This means that all but the last iteration through the loop produce the exact same values of A.
Instead of doing the same calculations over and over again, I retain the A array until the final Bernoulli number is produced.
<lang perl>#!perl use strict; use warnings; use List::Util qw(max); use Math::BigRat;
my $one = Math::BigRat->new(1); sub bernoulli_print { my @a; for my $m ( 0 .. 60 ) { push @a, $one / ($m + 1); for my $j ( reverse 1 .. $m ) { # This line: ( $a[$j-1] -= $a[$j] ) *= $j; # is a faster version of the following line: # $a[$j-1] = $j * ($a[$j-1] - $a[$j]); # since it avoids unnecessary object creation. } next unless $a[0]; printf "B(%2d) = %44s/%s\n", $m, $a[0]->parts; } }
bernoulli_print(); </lang> The output is exactly the same as the Python entry.
We can also use modules for faster results. E.g.
<lang perl>use ntheory qw/bernfrac/;
for my $n (0 .. 60) {
my($num,$den) = bernfrac($n); printf "B(%2d) = %44s/%s\n", $n, $num, $den if $num != 0;
}</lang> with identical output. Or: <lang perl>use Math::Pari qw/bernfrac/;
for my $n (0 .. 60) {
my($num,$den) = split "/", bernfrac($n);
printf("B(%2d) = %44s/%s\n", $n, $num, $den||1) if $num != 0;
}</lang> with the difference being that Pari chooses = -½.
Perl 6
First, a straighforward implementation of the naïve algorithm in the task description. <lang perl6>sub bernoulli($n) {
my @a;
for 0..$n -> $m {
@a[$m] = FatRat.new(1, $m + 1);
for reverse 1..$m -> $j {
@a[$j - 1] = $j * (@a[$j - 1] - @a[$j]);
}
}
return @a[0];
}
constant @bpairs = grep *.value.so, ($_ => bernoulli($_) for 0..60);
my $width = [max] @bpairs.map: *.value.numerator.chars; my $form = "B(%2d) = \%{$width}d/%d\n";
printf $form, .key, .value.nude for @bpairs;</lang>
- Output:
B( 0) = 1/1 B( 1) = 1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
Here is a much faster way, following the Perl solution that avoids recalculating previous values each time through the function. We do this in Perl 6 by not defining it as a function at all, but by defining it as an infinite sequence that we can read however many values we like from (52, in this case, to get up to B(100)). In this solution we've also avoided subscripting operations; rather we use a sequence operator (...) iterated over the list of the previous solution to find the next solution. We reverse the array in this case to make reference to the previous value in the list more natural, which means we take the last value of the list rather than the first value, and do so conditionally to avoid 0 values. <lang perl6>constant bernoulli = gather {
my @a;
for 0..* -> $m {
@a = FatRat.new(1, $m + 1),
-> $prev {
my $j = @a.elems;
$j * (@a.shift - $prev);
} ... { not @a.elems }
take $m => @a[*-1] if @a[*-1];
}
}
constant @bpairs = bernoulli[^52];
my $width = [max] @bpairs.map: *.value.numerator.chars; my $form = "B(%d)\t= \%{$width}d/%d\n";
printf $form, .key, .value.nude for @bpairs;</lang>
- Output:
B(0) = 1/1 B(1) = 1/2 B(2) = 1/6 B(4) = -1/30 B(6) = 1/42 B(8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730 B(62) = 12300585434086858541953039857403386151/6 B(64) = -106783830147866529886385444979142647942017/510 B(66) = 1472600022126335654051619428551932342241899101/64722 B(68) = -78773130858718728141909149208474606244347001/30 B(70) = 1505381347333367003803076567377857208511438160235/4686 B(72) = -5827954961669944110438277244641067365282488301844260429/140100870 B(74) = 34152417289221168014330073731472635186688307783087/6 B(76) = -24655088825935372707687196040585199904365267828865801/30 B(78) = 414846365575400828295179035549542073492199375372400483487/3318 B(80) = -4603784299479457646935574969019046849794257872751288919656867/230010 B(82) = 1677014149185145836823154509786269900207736027570253414881613/498 B(84) = -2024576195935290360231131160111731009989917391198090877281083932477/3404310 B(86) = 660714619417678653573847847426261496277830686653388931761996983/6 B(88) = -1311426488674017507995511424019311843345750275572028644296919890574047/61410 B(90) = 1179057279021082799884123351249215083775254949669647116231545215727922535/272118 B(92) = -1295585948207537527989427828538576749659341483719435143023316326829946247/1410 B(94) = 1220813806579744469607301679413201203958508415202696621436215105284649447/6 B(96) = -211600449597266513097597728109824233673043954389060234150638733420050668349987259/4501770 B(98) = 67908260672905495624051117546403605607342195728504487509073961249992947058239/6 B(100) = -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
And if you're a pure enough FP programmer to dislike destroying and reconstructing the array each time, here's the same algorithm without side effects. We use zip with the pair constructor => to keep values associated with their indices. This provides sufficient local information that we can define our own binary operator "bop" to reduce between each two terms, using the "triangle" form (called "scan" in Haskell) to return the intermediate results that will be important to compute the next Bernoulli number. Output is identical to the previous solution, but runs about 3.5 times slower in rakudo, as of this writing (2014-03). <lang perl6>my sub infix:<bop>(\prev,\this) { this.key => this.key * (this.value - prev.value) }
constant bernoulli = grep *.value, map { (.key => .value.[*-1]) }, do
0 => [FatRat.new(1,1)],
-> (:key($pm),:value(@pa)) {
$pm + 1 => [ map *.value, [\bop] ($pm + 2 ... 1) Z=> FatRat.new(1, $pm + 2), @pa ];
} ... *;</lang>
PicoLisp
Brute force and method by Srinivasa Ramanujan. <lang PicoLisp>(load "@lib/frac.l")
(de fact (N)
(cache '(NIL) N
(if (=0 N) 1 (apply * (range 1 N))) ) )
(de binomial (N K)
(frac
(/
(fact N)
(* (fact (- N K)) (fact K)) )
1 ) )
(de A (N M)
(let Sum (0 . 1)
(for X M
(setq Sum
(f+
Sum
(f*
(binomial (+ N 3) (- N (* X 6)))
(berno (- N (* X 6)) ) ) ) ) )
Sum ) )
(de berno (N)
(cache '(NIL) N
(cond
((=0 N) (1 . 1))
((= 1 N) (-1 . 2))
((bit? 1 N) (0 . 1))
(T
(case (% N 6)
(0
(f/
(f-
(frac (+ N 3) 3)
(A N (/ N 6)) )
(binomial (+ N 3) N) ) )
(2
(f/
(f-
(frac (+ N 3) 3)
(A N (/ (- N 2) 6)) )
(binomial (+ N 3) N) ) )
(4
(f/
(f-
(f* (-1 . 1) (frac (+ N 3) 6))
(A N (/ (- N 4) 6)) )
(binomial (+ N 3) N) ) ) ) ) ) ) )
(de berno-brute (N)
(cache '(NIL) N
(let Sum (0 . 1)
(cond
((=0 N) (1 . 1))
((= 1 N) (-1 . 2))
((bit? 1 N) (0 . 1))
(T
(for (X 0 (> N X) (inc X))
(setq Sum
(f+
Sum
(f* (binomial (inc N) X) (berno-brute X)) ) ) )
(f/ (f* (-1 . 1) Sum) (binomial (inc N) N)) ) ) ) ) )
(for (N 0 (> 62 N) (inc N))
(if (or (= N 1) (not (bit? 1 N)))
(tab (2 4 -60) N " => " (sym (berno N))) ) )
(for (N 0 (> 400 N) (inc N))
(test (berno N) (berno-brute N)) )
(bye)</lang>
PL/I
<lang PL/I>Bern: procedure options (main); /* 4 July 2014 */
declare i fixed binary; declare B complex fixed (31);
Bernoulli: procedure (n) returns (complex fixed (31));
declare n fixed binary; declare anum(0:n) fixed (31), aden(0:n) fixed (31); declare (j, m) fixed; declare F fixed (31);
do m = 0 to n;
anum(m) = 1;
aden(m) = m+1;
do j = m to 1 by -1;
anum(j-1) = j*( aden(j)*anum(j-1) - aden(j-1)*anum(j) );
aden(j-1) = ( aden(j-1) * aden(j) );
F = gcd(abs(anum(j-1)), abs(aden(j-1)) );
if F ^= 1 then
do;
anum(j-1) = anum(j-1) / F;
aden(j-1) = aden(j-1) / F;
end;
end;
end;
return ( complex(anum(0), aden(0)) );
end Bernoulli;
do i = 0, 1, 2 to 36 by 2; /* 36 is upper limit imposed by hardware. */
B = Bernoulli(i);
put skip edit ('B(' , trim(i) , ')=' , real(B) , '/' , trim(imag(B)) )
(3 A, column(10), F(32), 2 A);
end;
end Bern;</lang> The above uses GCD (see Rosetta Code) extended for 31-digit working. Results obtained by this program:
B(0)= 1/1 B(1)= 1/2 B(2)= 1/6 B(4)= -1/30 B(6)= 1/42 B(8)= -1/30 B(10)= 5/66 B(12)= -691/2730 B(14)= 7/6 B(16)= -3617/510 B(18)= 43867/798 B(20)= -174611/330 B(22)= 854513/138 B(24)= -236364091/2730 B(26)= 8553103/6 B(28)= -23749461029/870 B(30)= 8615841276005/14322 B(32)= -7709321041217/510 B(34)= 2577687858367/6 B(36)= -26315271553053477373/1919190
Python
Python: Using task algorithm
<lang python>from fractions import Fraction as Fr
def bernoulli(n):
A = [0] * (n+1)
for m in range(n+1):
A[m] = Fr(1, m+1)
for j in range(m, 0, -1):
A[j-1] = j*(A[j-1] - A[j])
return A[0] # (which is Bn)
bn = [(i, bernoulli(i)) for i in range(61)] bn = [(i, b) for i,b in bn if b] width = max(len(str(b.numerator)) for i,b in bn) for i,b in bn:
print('B(%2i) = %*i/%i' % (i, width, b.numerator, b.denominator))</lang>
- Output:
B( 0) = 1/1 B( 1) = 1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
Python: Optimised task algorithm
Using the optimization mentioned in the Perl entry to reduce intermediate calculations we create and use the generator bernoulli2(): <lang python>def bernoulli2():
A, m = [], 0
while True:
A.append(Fr(1, m+1))
for j in range(m, 0, -1):
A[j-1] = j*(A[j-1] - A[j])
yield A[0] # (which is Bm)
m += 1
bn2 = [ix for ix in zip(range(61), bernoulli2())] bn2 = [(i, b) for i,b in bn2 if b] width = max(len(str(b.numerator)) for i,b in bn2) for i,b in bn2:
print('B(%2i) = %*i/%i' % (i, width, b.numerator, b.denominator))</lang>
Output is exactly the same as before.
Racket
This implements, firstly, the algorithm specified with the task... then the better performing bernoulli.3, which uses the "double sum formula" listed under REXX. The number generators all (there is also a bernoulli.2) use the same emmitter... it's just a matter of how long to wait for the emission.
<lang>#lang racket
- As described in task...
(define (bernoulli.1 n)
(define A (make-vector (add1 n)))
(for ((m (in-range 0 (add1 n))))
(vector-set! A m (/ (add1 m)))
(for ((j (in-range m (sub1 1) -1)))
(define new-A_j-1 (* j (- (vector-ref A (sub1 j)) (vector-ref A j))))
(vector-set! A (sub1 j) new-A_j-1)))
(vector-ref A 0))
(define (non-zero-bernoulli-indices s)
(sequence-filter (λ (n) (or (even? n) (= n 1))) s))
(define (bernoulli_0..n B N)
(for/list ((n (non-zero-bernoulli-indices (in-range (add1 N))))) (B n)))
- From REXX description / http://mathworld.wolfram.com/BernoulliNumber.html #33
- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- bernoulli.2 is for illustrative purposes, binomial is very costly if there is no memoisation
- (which math/number-theory doesn't do)
(require (only-in math/number-theory binomial)) (define (bernoulli.2 n)
(for/sum ((k (in-range 0 (add1 n))))
(* (/ (add1 k))
(for/sum ((r (in-range 0 (add1 k))))
(* (expt -1 r) (binomial k r) (expt r n))))))
- Three things to do
- 1. (expt -1 r)
- is 1 for even r, -1 for odd r... split the sum between those two.
- 2. splitting the sum might has arithmetic advantages, too. We're using rationals, so the smaller
- summations should require less normalisation of intermediate, fractional results
- 3. a memoised binomial... although the one from math/number-theory is fast, it is (and its
- factorials are) computed every time which is redundant
(define kCr-memo (make-hasheq)) (define !-memo (make-vector 1000 #f)) (vector-set! !-memo 0 1) ;; seed the memo (define (! k)
(cond [(vector-ref !-memo k) => values]
[else (define k! (* k (! (- k 1)))) (vector-set! !-memo k k!) k!]))
(define (kCr k r)
; If we want (kCr ... r>1000000) we'll have to reconsider this. However, until then... (define hash-key (+ (* 1000000 k) r)) (hash-ref! kCr-memo hash-key (λ () (/ (! k) (! r) (! (- k r))))))
(define (bernoulli.3 n)
(for/sum ((k (in-range 0 (add1 n))))
(define k+1 (add1 k))
(* (/ k+1)
(- (for/sum ((r (in-range 0 k+1 2))) (* (kCr k r) (expt r n)))
(for/sum ((r (in-range 1 k+1 2))) (* (kCr k r) (expt r n)))))))
(define (display/align-fractions caption/idx-fmt Bs)
;; widths are one more than the order of magnitude
(define oom+1 (compose add1 order-of-magnitude))
(define-values (I-width N-width D-width)
(for/fold ((I 0) (N 0) (D 0))
((b Bs) (n (non-zero-bernoulli-indices (in-naturals))))
(define +b (abs b))
(values (max I (oom+1 (max n 1)))
(max N (+ (oom+1 (numerator +b)) (if (negative? b) 1 0)))
(max D (oom+1 (denominator +b))))))
(define (~a/w/a n w a) (~a n #:width w #:align a))
(for ((n (non-zero-bernoulli-indices (in-naturals))) (b Bs))
(printf "~a ~a/~a~%"
(format caption/idx-fmt (~a/w/a n I-width 'right))
(~a/w/a (numerator b) N-width 'right)
(~a/w/a (denominator b) D-width 'left))))
(module+ main
(display/align-fractions "B(~a) =" (bernoulli_0..n bernoulli.3 60)))
(module+ test
(require rackunit)
; correctness and timing tests
(check-match (time (bernoulli_0..n bernoulli.1 60))
(list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...))
(check-match (time (bernoulli_0..n bernoulli.2 60))
(list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...))
(check-match (time (bernoulli_0..n bernoulli.3 60))
(list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...))
; timing only ...
(void (time (bernoulli_0..n bernoulli.3 100))))</lang>
- Output:
B( 0) = 1/1 B( 1) = -1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
REXX
The double sum formula used is number (33) from the entry Bernoulli number on Wolfram MathWorldTM.
- where is a binomial coefficient.
- where is a binomial coefficient.
<lang rexx>/*REXX program calculates a number of Bernoulli numbers (as fractions). */ parse arg N .; if N== then N=60 /*get N. If ¬ given, use default*/ !.=0; w=max(length(N),4); Nw=N+N%5 /*used for aligning the output. */ say 'B(n)' center('Bernoulli number expressed as a fraction', max(78-w,Nw)) say copies('─',w) copies('─', max(78-w,Nw+2*w))
do #=0 to N /*process numbers from 0 ──► N. */
b=bern(#); if b==0 then iterate /*calculate Bernoulli#, skip if 0*/
indent=max(0, nW-pos('/', b)) /*calculate alignment indentation*/
say right(#,w) left(,indent) b /*display the indented Bernoulli#*/
end /*#*/ /* [↑] align the Bernoulli number*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────BERN subroutine─────────────────────*/ bern: parse arg x /*obtain the subroutine argument.*/ if x==0 then return '1/1' /*handle the special case of zero*/ if x==1 then return '-1/2' /* " " " " " one.*/ if x//2 then return 0 /* " " " " " odds*/
/* [↓] process all #s up to X, */
do j=2 to x by 2; jp=j+1; d=j+j /* & set some shortcut vars.*/
if d>digits() then numeric digits d /*increase precision if needed. */
sn=1-j /*set the numerator. */
sd=2 /* " " denominator. */
do k=2 to j-1 by 2 /*calculate a SN/SD sequence. */
parse var @.k bn '/' ad /*get a previously calculated fra*/
an=comb(jp,k)*bn /*use COMBination for next term. */
lcm=lcm.(sd,ad) /*use Least Common Denominator. */
sn=lcm%sd*sn; sd=lcm /*calculate current numerator. */
an=lcm%ad*an; ad=lcm /* " next " */
sn=sn+an /* " current " */
end /*k*/ /* [↑] calculate SN/SD sequence.*/
sn=-sn /*adjust the sign for numerator. */
sd=sd*jp /*calculate the denomitator. */
if sn\==1 then do /*reduce the fraction if possible*/
_=gcd.(sn,sd) /*get Greatest Common Denominator*/
sn=sn%_; sd=sd%_ /*reduce numerator & denominator.*/
end /* [↑] done with the reduction.*/
@.j=sn'/'sd /*save the result for next round.*/
end /*j*/ /* [↑] done with calculating B#.*/
return sn'/'sd /*──────────────────────────────────COMB subroutine─────────────────────*/ comb: procedure expose !.; parse arg x,y; if x==y then return 1 if !.!c.x.y\==0 then return !.!c.x.y /*combination computed before ? */
if x-y<y then y=x-y; z=perm(x,y); do j=2 to y; z=z%j; end
!.!c.x.y=z; return z /*assign memoization; then return*/ /*──────────────────────────────────GCD. subroutine (simplified)────────*/ gcd.: procedure; parse arg x,y; x=abs(x)
do until _==0; _=x//y; x=y; y=_; end; return x
/*──────────────────────────────────LCM. subroutine (simplified)────────*/ lcm.: procedure; parse arg x,y; x=abs(x); return x*y/gcd.(x,y) /*──────────────────────────────────PERM subroutine─────────────────────*/ perm: procedure expose !.; parse arg x,y; z=1 if !.!p.x.y\==0 then return !.!p.x.y /*permutation computed before ? */
do j=x-y+1 to x; z=z*j; end; !.!p.x.y=z; return z</lang>
output when using the default input:
B(n) Bernoulli number expressed as a fraction ──── ──────────────────────────────────────────────────────────────────────────────── 0 1/1 1 -1/2 2 1/6 4 -1/30 6 1/42 8 -1/30 10 5/66 12 -691/2730 14 7/6 16 -3617/510 18 43867/798 20 -174611/330 22 854513/138 24 -236364091/2730 26 8553103/6 28 -23749461029/870 30 8615841276005/14322 32 -7709321041217/510 34 2577687858367/6 36 -26315271553053477373/1919190 38 2929993913841559/6 40 -261082718496449122051/13530 42 1520097643918070802691/1806 44 -27833269579301024235023/690 46 596451111593912163277961/282 48 -5609403368997817686249127547/46410 50 495057205241079648212477525/66 52 -801165718135489957347924991853/1590 54 29149963634884862421418123812691/798 56 -2479392929313226753685415739663229/870 58 84483613348880041862046775994036021/354 60 -1215233140483755572040304994079820246041491/56786730
Ruby
<lang ruby>bernoulli = Enumerator.new do |y|
ar, m = [], 0
loop do
ar << Rational(1, m+1)
m.downto(1){|j| ar[j-1] = j*(ar[j-1] - ar[j]) }
y << ar.first # yield
m += 1
end
end
b_nums = bernoulli.take(61) width = b_nums.map{|b| b.numerator.to_s.size}.max b_nums.each_with_index {|b,i| puts "B(%2i) = %*i/%i" % [i, width, b.numerator, b.denominator] unless b.zero? }
</lang>
- Output:
B( 0) = 1/1 B( 1) = 1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
Tcl
<lang tcl>proc bernoulli {n} {
for {set m 0} {$m <= $n} {incr m} {
lappend A [list 1 [expr {$m + 1}]] for {set j $m} {[set i $j] >= 1} {} { lassign [lindex $A [incr j -1]] a1 b1 lassign [lindex $A $i] a2 b2 set x [set p [expr {$i * ($a1*$b2 - $a2*$b1)}]] set y [set q [expr {$b1 * $b2}]] while {$q} {set q [expr {$p % [set p $q]}]} lset A $j [list [expr {$x/$p}] [expr {$y/$p}]] }
} return [lindex $A 0]
}
set len 0 for {set n 0} {$n <= 60} {incr n} {
set b [bernoulli $n]
if {[lindex $b 0]} {
lappend result $n {*}$b set len [expr {max($len, [string length [lindex $b 0]])}]
}
} foreach {n num denom} $result {
puts [format {B_%-2d = %*lld/%lld} $n $len $num $denom]
}</lang>
- Output:
B_0 = 1/1 B_1 = 1/2 B_2 = 1/6 B_4 = -1/30 B_6 = 1/42 B_8 = -1/30 B_10 = 5/66 B_12 = -691/2730 B_14 = 7/6 B_16 = -3617/510 B_18 = 43867/798 B_20 = -174611/330 B_22 = 854513/138 B_24 = -236364091/2730 B_26 = 8553103/6 B_28 = -23749461029/870 B_30 = 8615841276005/14322 B_32 = -7709321041217/510 B_34 = 2577687858367/6 B_36 = -26315271553053477373/1919190 B_38 = 2929993913841559/6 B_40 = -261082718496449122051/13530 B_42 = 1520097643918070802691/1806 B_44 = -27833269579301024235023/690 B_46 = 596451111593912163277961/282 B_48 = -5609403368997817686249127547/46410 B_50 = 495057205241079648212477525/66 B_52 = -801165718135489957347924991853/1590 B_54 = 29149963634884862421418123812691/798 B_56 = -2479392929313226753685415739663229/870 B_58 = 84483613348880041862046775994036021/354 B_60 = -1215233140483755572040304994079820246041491/56786730