Bernoulli numbers: Difference between revisions
(→{{header|Perl 6}}: better explanation, line up slash even if we go to B with 3 or 4 digits) |
(looks like a real task to me now) |
||
| Line 1: | Line 1: | ||
{{ |
{{task}} |
||
'''Bernoulli numbers''' |
'''Bernoulli numbers''' |
||
Revision as of 05:02, 12 March 2014
You are encouraged to solve this task according to the task description, using any language you may know.
Bernoulli numbers
Bernoulli numbers are used in some series expansions of serval functions (trigonometric, hyperbolic, gamma, etc.), and are extremely important in number theory and analysis.
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), and are expressed as Bn.
task requirements
- show the Bernoulli numbers B0 through B60.
- express the numbers as fractions (most are improper fractions).
- fractions should be reduced.
- suppress all 0 (zero) values.
- index each number in some way so that it can be discerned which number is being displayed.
- align the solidi [/] if used (extra credit).
(Other than B1, all odd Bernoulli have a value of 0 (zero).
- 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).
- Wiki entry Bernoulli number.
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!(i => tuple(i, i.bernoulli))
.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.
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.
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
Python
<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
REXX
The double sum formula used is #33 from the entry Bernoulli number on The Eric Weisstein's World of Mathematics (TM).
- 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*/ 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; parse arg x,y; if x==y then return 1
if x-y<y then y=x-y; z=perm(x,y); do j=2 to y; z=z/j;end; return z
/*──────────────────────────────────GCD subroutine──────────────────────*/ gcd: procedure; $=; do i=1 for arg(); $=$ arg(i); end; parse var $ x z .
if x=0 then x=z; x=abs(x); do j=2 to words($); y=abs(word($,j))
if y=0 then iterate; do until _==0; _=x//y; x=y; y=_; end; end
return x
/*──────────────────────────────────LCM subroutine──────────────────────*/ lcm: procedure; $=; do j=1 for arg(); $=$ arg(j); end; x=abs(word($,1))
do k=2 to words($); !=abs(word($,k)); if !=0 then return 0
x=x * ! / gcd(x,!); end
return x
/*──────────────────────────────────PERM subroutine─────────────────────*/ perm: procedure; parse arg x,y; z=1; do j=x-y+1 to x; z=z*j; end; 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