Bernoulli numbers: Difference between revisions
m
→{{header|Fōrmulæ}}
m (Added the Sidef language) |
|||
| (181 intermediate revisions by 57 users not shown) | |||
Line 1:
{{task|Mathematics}}
[[wp:Bernoulli number|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 n<sup>th</sup> Bernoulli number is expressed as '''B'''<sub>n</sub>.
<br>
;Task
:* show the Bernoulli numbers '''B'''<sub>0</sub> through '''B'''<sub>60</sub>.
:* suppress the output of values which are equal to zero. (Other than '''B'''<sub>1</sub> , all ''odd'' Bernoulli numbers have a value of zero.)
:* express the Bernoulli numbers as fractions (most are improper fractions).
:* the fractions should be reduced.
:* index each number in some way so that it can be discerned which Bernoulli number is being displayed.
:* align the solidi (<big><b>/</b></big>) if used (extra credit).
;An algorithm
Line 19 ⟶ 27:
;See also
* Sequence [
* Sequence [
* Entry [http://mathworld.wolfram.com/BernoulliNumber.html Bernoulli number] on The Eric Weisstein's World of Mathematics (TM).
* Luschny's [http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html The Bernoulli Manifesto] for a discussion on <
<br><br>
=={{header|Ada}}==
Using a GMP thick binding available at http://www.codeforge.com/article/422541
<syntaxhighlight lang="ada">WITH GMP.Rationals, GMP.Integers, Ada.Text_IO, Ada.Strings.Fixed, Ada.Strings;
USE GMP.Rationals, GMP.Integers, Ada.Text_IO, Ada.Strings.Fixed, Ada.Strings;
PROCEDURE Main IS
FUNCTION Bernoulli_Number (N : Natural) RETURN Unbounded_Fraction IS
FUNCTION "/" (Left, Right : Natural) RETURN Unbounded_Fraction IS
(To_Unbounded_Integer (Left) / To_Unbounded_Integer (Right));
A : ARRAY (0 .. N) OF Unbounded_Fraction;
BEGIN
FOR M IN 0 .. N LOOP
A (M) := 1 / (M + 1);
FOR J IN REVERSE 1 .. M LOOP
A (J - 1) := (J / 1 ) * (A (J - 1) - A (J));
END LOOP;
END LOOP;
RETURN A (0);
END Bernoulli_Number;
BEGIN
FOR I IN 0 .. 60 LOOP
IF I MOD 2 = 0 OR I = 1 THEN
DECLARE
B : Unbounded_Fraction := Bernoulli_Number (I);
S : String := Image (GMP.Rationals.Numerator (B));
BEGIN
Put_Line ("B (" & (IF I < 10 THEN " " ELSE "") & Trim (I'Img, Left)
& ")=" & (44 - S'Length) * " " & Image (B));
END;
END IF;
END LOOP;
END Main;</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Uses the LONG LONG INT mode of Algol 68G which allows large precision integers.
<syntaxhighlight lang="algol68">BEGIN
# Show Bernoulli numbers B0 to B60 as rational numbers #
# Uses code from the Arithmetic/Rational task modified to use #
# LONG LONG INT to allow for the large number of digits requried #
PR precision 100 PR # sets the precision of LONG LONG INT #
# Code from the Arithmetic/Rational task #
# ============================================================== #
MODE FRAC = STRUCT( LONG LONG INT num #erator#, den #ominator#);
PROC gcd = (LONG LONG INT a, b) LONG LONG INT: # greatest common divisor #
(a = 0 | b |: b = 0 | a |: ABS a > ABS b | gcd(b, a MOD b) | gcd(a, b MOD a));
PROC lcm = (LONG LONG INT a, b)LONG LONG INT: # least common multiple #
a OVER gcd(a, b) * b;
PRIO // = 9; # higher then the ** operator #
OP // = (LONG LONG INT num, den)FRAC: ( # initialise and normalise #
LONG LONG INT common = gcd(num, den);
IF den < 0 THEN
( -num OVER common, -den OVER common)
ELSE
( num OVER common, den OVER common)
FI
);
OP + = (FRAC a, b)FRAC: (
LONG LONG INT common = lcm(den OF a, den OF b);
FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
num OF result//den OF result
);
OP - = (FRAC a, b)FRAC: a + -b,
* = (FRAC a, b)FRAC: (
LONG LONG INT num = num OF a * num OF b,
den = den OF a * den OF b;
LONG LONG INT common = gcd(num, den);
(num OVER common) // (den OVER common)
);
OP - = (FRAC frac)FRAC: (-num OF frac, den OF frac);
# ============================================================== #
# end code from the Arithmetic/Rational task #
# Additional FRACrelated operators #
OP * = ( INT a, FRAC b )FRAC: ( num OF b * a ) // den OF b;
OP // = ( INT a, INT b )FRAC: LONG LONG INT( a ) // LONG LONG INT( b );
# returns the nth Bernoulli number, n must be >= 0 #
# Uses the algorithm suggested by the task, so B(1) is +1/2 #
PROC bernoulli = ( INT n )FRAC:
IF n < 0
THEN # n is out of range # 0 // 1
ELSE # n is valid #
[ 0 : n ]FRAC a;
FOR i FROM LWB a TO UPB a DO a[ i ] := 0 // 1 OD;
FOR m FROM 0 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 ] )
OD
OD;
a[ 0 ]
FI # bernoulli # ;
FOR n FROM 0 TO 60 DO
FRAC bn := bernoulli( n );
IF num OF bn /= 0 THEN
# have a non-0 Bn #
print( ( "B(", whole( n, -2 ), ") ", whole( num OF bn, -50 ), " / ", whole( den OF bn, 0 ), newline ) )
FI
OD
END
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|AppleScript}}==
To be able to handle numbers up to B(60) and beyond, this script represents the numbers with lists whose items are the digit values — which of course requires custom math routines.
<syntaxhighlight lang="applescript">on bernoullis(n) -- Return a list of "numerator / denominator" texts representing Bernoulli numbers B(0) to B(n).
set listMathScript to getListMathScript(10) -- Script object providing custom list math routines.
set output to {}
-- Akiyama–Tanigawa algorithm for the "second Bernoulli numbers".
-- List 'a' will contain {numerator, denominator} lists representing fractions.
-- The numerators and denominators will in turn be lists containing integers representing their (decimal) digits.
set a to {}
repeat with m from 0 to n
-- Append the structure for 1 / (m + 1) to the end of a.
set {numerator2, denominator2} to {{1}, listMathScript's intToList(m + 1)}
set a's end to result
repeat with j from m to 1 by -1
-- Retrieve the preceding numerator and denominator.
set {numerator1, denominator1} to a's item j
tell listMathScript
-- Get the two fractions' lowest common denominator and adjust the numerators accordingly.
set lcd to its lcm(denominator1, denominator2)
set numerator1 to its multiply(numerator1, its |div|(lcd, denominator1))
set numerator2 to its multiply(numerator2, its |div|(lcd, denominator2))
-- Subtract numerator2 from numerator1 and multiply the result by j.
-- Assign the results to numerator2 and denominator2 for the next iteration.
set numerator2 to its multiply(its subtract(numerator1, numerator2), its intToList(j))
set denominator2 to lcd
end tell
-- Also store them in a's slot j. No need to reduce them here.
set a's item j to {numerator2, denominator2}
end repeat
-- The fraction just stored in a's first slot is Bernoulli(m). Reduce it and append a text representation to the output.
tell listMathScript
set gcd to its hcf(numerator2, denominator2)
set numerator2 to its |div|(numerator2, gcd)
set denominator2 to its |div|(denominator2, gcd)
set end of output to its listToText(numerator2) & (" / " & its listToText(denominator2))
end tell
end repeat
return output
end bernoullis
on getListMathScript(base)
script
on multiply(lst1, lst2) -- Multiply lst1 by lst2.
set lst1Length to (count lst1)
set lst2Length to (count lst2)
set productLength to lst1Length + lst2Length - 1
set product to {}
repeat productLength times
set product's end to 0
end repeat
-- Long multiplication algorithm, updating product digits on the fly instead of summing rows at the end.
repeat with lst2Index from -1 to -lst2Length by -1
set lst2Digit to lst2's item lst2Index
if (lst2Digit is not 0) then
set carry to 0
set productIndex to lst2Index
repeat with lst1Index from lst1's length to 1 by -1
tell lst2Digit * (lst1's item lst1Index) + carry + (product's item productIndex)
set product's item productIndex to (it mod base)
set carry to (it div base)
end tell
set productIndex to productIndex - 1
end repeat
if (carry = 0) then
else if (productIndex < -productLength) then
set product's beginning to carry
else
set product's item productIndex to (product's item productIndex) + carry
end if
end if
end repeat
return product
end multiply
on subtract(lst1, lst2) -- Subtract lst2 from lst1.
set lst1Length to (count lst1)
set lst2Length to (count lst2)
-- Pad copies to equal lengths.
copy lst1 to lst1
repeat (lst2Length - lst1Length) times
set lst1's beginning to 0
end repeat
copy lst2 to lst2
repeat (lst1Length - lst2Length) times
set lst2's beginning to 0
end repeat
-- Is lst2's numeric value greater than lst1's?
set paddedLength to (count lst1)
repeat with i from 1 to paddedLength
set lst1Digit to lst1's item i
set lst2Digit to lst2's item i
set lst2Greater to (lst2Digit > lst1Digit)
if ((lst2Greater) or (lst1Digit > lst2Digit)) then exit repeat
end repeat
-- If so, set up to subtract lst1 from lst2 instead. We'll invert the result's sign at the end.
if (lst2Greater) then tell lst2
set lst2 to lst1
set lst1 to it
end tell
-- The subtraction at last!
set difference to {}
set borrow to 0
repeat with i from paddedLength to 1 by -1
tell (lst1's item i) + base - borrow - (lst2's item i)
set difference's beginning to (it mod base)
set borrow to 1 - (it div base)
end tell
end repeat
if (lst2Greater) then invert(difference)
return difference
end subtract
on |div|(lst1, lst2) -- List lst1 div lst2.
return divide(lst1, lst2)'s quotient
end |div|
on |mod|(lst1, lst2) -- List lst1 mod lst2.
return divide(lst1, lst2)'s remainder
end |mod|
on divide(lst1, lst2) -- Divide lst1 by lst2. Return a record containing separate lists for the quotient and remainder.
set dividend to trim(lst1)
set divisor to trim(lst2)
set dividendLength to (count dividend)
set divisorLength to (count divisor)
if (divisorLength > dividendLength) then return {quotient:{0}, remainder:dividend}
-- Note the dividend's and divisor's signs, but use absolute values in the division.
set dividendNegative to (dividend's beginning < 0)
if (dividendNegative) then invert(dividend)
set divisorNegative to (divisor's beginning < 0)
if (divisorNegative) then invert(divisor)
-- Long-division algorithm, but quotient digits are subtraction counts.
set quotient to {}
if (divisorLength > 1) then
set remainder to dividend's items 1 thru (divisorLength - 1)
else
set remainder to {}
end if
repeat with nextSlot from divisorLength to dividendLength
set remainder's end to dividend's item nextSlot
repeat with subtractionCount from 0 to base -- Only ever reaches base - 1.
set subtractionResult to trim(subtract(remainder, divisor))
if (subtractionResult's beginning < 0) then exit repeat
set remainder to subtractionResult
end repeat
set end of quotient to subtractionCount
end repeat
-- The quotient's negative if the input signs are different. Positive otherwise.
if (dividendNegative ≠ divisorNegative) then invert(quotient)
-- The remainder has the same sign as the dividend.
if (dividendNegative) then invert(remainder)
return {quotient:quotient, remainder:remainder}
end divide
on lcm(lst1, lst2) -- Lowest common multiple of lst1 and lst2.
return multiply(lst2, |div|(lst1, hcf(lst1, lst2)))
end lcm
on hcf(lst1, lst2) -- Highest common factor of lst1 and lst2.
set lst1 to trim(lst1)
set lst2 to trim(lst2)
repeat until (lst2 = {0})
set x to lst1
set lst1 to lst2
set lst2 to trim(|mod|(x, lst2))
end repeat
if (lst1's beginning < 0) then invert(lst1)
return lst1
end hcf
on invert(lst) -- Invert the sign of all lst's "digits".
repeat with thisDigit in lst
set thisDigit's contents to -thisDigit
end repeat
end invert
on trim(lst) -- Return a copy of lst with no leading zeros.
repeat with i from 1 to (count lst)
if (lst's item i is not 0) then exit repeat
end repeat
return lst's items i thru end
end trim
on intToList(n) -- Return a list of numbers representing n's digits.
set lst to {n mod base}
set n to n div base
repeat until (n = 0)
set beginning of lst to n mod base as integer
set n to n div base
end repeat
return lst
end intToList
on listToText(lst) -- Return the number represented by the input list as text.
-- This lazily assumes 2 <= base <= 10. :)
set lst to trim(lst)
if (lst's beginning < 0) then
invert(lst)
set lst's beginning to "-"
end if
return join(lst, "")
end listToText
end script
return result
end getListMathScript
on join(lst, delim)
set astid to AppleScript's text item delimiters
set AppleScript's text item delimiters to delim
set txt to lst as text
set AppleScript's text item delimiters to astid
return txt
end join
on task()
set maxN to 60
set output to {""}
set padding to " = "
set bernoulliNumbers to bernoullis(maxN)
repeat with n from 0 to maxN
set bernie to bernoulliNumbers's item (n + 1)
if (bernie does not start with "0") then
set Bn to "B(" & n & ")"
set output's end to Bn & ¬
text ((count Bn) - 3) thru (50 - (offset of "/" in bernie)) of padding & ¬
bernie
end if
end repeat
return join(output, linefeed)
end task
task()</syntaxhighlight>
{{output}}
<syntaxhighlight lang="applescript">"
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"</syntaxhighlight>
=={{header|Bracmat}}==
<syntaxhighlight lang="bracmat"> ( BernoulliList
= B Bs answer indLn indexLen indexPadding
, n numberPadding p solPos solidusPos sp
. ( B
= m A a j b
. -1:?m
& :?A
& whl
' ( 1+!m:~>!arg:?m
& ((!m+1:?j)^-1:?a)
map
$ ( (
= .(-1+!j:?j)*(!arg+-1*!a):?a
)
. !A
)
: ?A
)
& !A:? @?b
& !b
)
& -1:?n
& :?Bs
& whl
' ( 1+!n:~>!arg:?n
& B$!n !Bs:?Bs
)
& @(!arg:? [?indexLen)
& 1+!indexLen:?indexLen
& !Bs:%@(?:? "/" [?solidusPos ?) ?
& 1+!solidusPos:?solidusPos:?p
& :?sp
& whl
' (!p+-1:~<0:?p&" " !sp:?sp)
& :?answer
& whl
' ( !Bs:%?B ?Bs
& ( !B:0
| (!B:/|str$(!B "/1"):?B)
& @(!B:? "/" [?solPos ?)
& @(!arg:? [?indLn)
& !sp
: ? [(-1*!indexLen+!indLn) ?indexPadding
: ? [(-1*!solidusPos+!solPos) ?numberPadding
& "B("
!arg
")="
!indexPadding
!numberPadding
(!B:>0&" "|)
!B
\n
!answer
: ?answer
)
& -1+!arg:?arg
)
& str$!answer
)
& BernoulliList$60;</syntaxhighlight>
<pre>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</pre>
=={{header|C}}==
{{libheader|GMP}}
<syntaxhighlight lang="c">
#include <stdlib.h>
#include <gmp.h>
Line 78 ⟶ 646:
return 0;
}
</syntaxhighlight>
{{out}}
<pre>
Line 114 ⟶ 682:
B(60) = -1215233140483755572040304994079820246041491 / 56786730
</pre>
=={{header|C sharp|C#}}==
Line 120 ⟶ 687:
{{libheader|Mpir.NET}}
Translation of the C implementation
<
using Mpir.NET;
using System;
Line 174 ⟶ 741:
}
}
</syntaxhighlight>
{{out}}
<pre>
Line 212 ⟶ 779:
=== Using Math.NET ===
{{libheader|MathNet.Numerics}}
<
using System;
using System.Console;
Line 256 ⟶ 824:
static void Main()
{
Enumerable.Range(0,
.Select(n => new {N = n, BernoulliNumber = CalculateBernoulli(n)})
.Where(b => !b.BernoulliNumber.Numerator.IsZero)
Line 265 ⟶ 833:
}
}
</syntaxhighlight>
{{out}}
<pre>
Line 299 ⟶ 867:
B(56) = -2479392929313226753685415739663229 / 870
B(58) = 84483613348880041862046775994036021 / 354
B(60) = -1215233140483755572040304994079820246041491 / 56786730
</pre>
=== Using System.Numerics ===
{{libheader|System.Numerics}}
Algo based on the example provided in the header of this RC page (the one from Wikipedia). <br/> Extra feature - one can override the default of 60 by supplying a suitable number on the command line. The column widths are not hard-coded, but will adapt to the widths of the items listed.
<syntaxhighlight lang="csharp">using System;
using System.Numerics;
using System.Collections.Generic;
namespace bern
{
class Program
{
struct BerNum { public int index; public BigInteger Numer, Denomin; };
static int w1 = 1, w2 = 1; // widths for formatting output
static int max = 60; // default maximum, can override on command line
// returns nth Bernoulli number
static BerNum CalcBernoulli(int n)
{
BerNum res;
BigInteger f;
BigInteger[] nu = new BigInteger[n + 1],
de = new BigInteger[n + 1];
for (int m = 0; m <= n; m++)
{
nu[m] = 1; de[m] = m + 1;
for (int j = m; j > 0; j--)
if ((f = BigInteger.GreatestCommonDivisor(
nu[j - 1] = j * (de[j] * nu[j - 1] - de[j - 1] * nu[j]),
de[j - 1] *= de[j])) != BigInteger.One)
{ nu[j - 1] /= f; de[j - 1] /= f; }
}
res.index = n; res.Numer = nu[0]; res.Denomin = de[0];
w1 = Math.Max(n.ToString().Length, w1); // ratchet up widths
w2 = Math.Max(res.Numer.ToString().Length, w2);
if (max > 50) Console.Write("."); // progress dots appear for larger values
return res;
}
static void Main(string[] args)
{
List<BerNum> BNumbList = new List<BerNum>();
// defaults to 60 when no (or invalid) command line parameter is present
if (args.Length > 0) {
int.TryParse(args[0], out max);
if (max < 1 || max > Int16.MaxValue) max = 60;
if (args[0] == "0") max = 0;
}
for (int i = 0; i <= max; i++) // fill list with values
{
BerNum BNumb = CalcBernoulli(i);
if (BNumb.Numer != BigInteger.Zero) BNumbList.Add(BNumb);
}
if (max > 50) Console.WriteLine();
string strFmt = "B({0, " + w1.ToString() + "}) = {1, " + w2.ToString() + "} / {2}";
// display formatted list
foreach (BerNum bn in BNumbList)
Console.WriteLine(strFmt , bn.index, bn.Numer, bn.Denomin);
if (System.Diagnostics.Debugger.IsAttached) Console.Read();
}
}
}
</syntaxhighlight>
{{out}}
Default (nothing entered on command line):
<pre>.............................................................
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</pre>
Output with "8" entered on command line:
<pre>B(0) = 1 / 1
B(1) = 1 / 2
B(2) = 1 / 6
B(4) = -1 / 30
B(6) = 1 / 42
B(8) = -1 / 30</pre>
Output with "126" entered on the command line:
<pre style="height:64ex;overflow:scroll">...............................................................................................................................
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
B(102) = 3204019410860907078243020782116241775491817197152717450679002501086861530836678158791 / 4326
B(104) = -319533631363830011287103352796174274671189606078272738327103470162849568365549721224053 / 1590
B(106) = 36373903172617414408151820151593427169231298640581690038930816378281879873386202346572901 / 642
B(108) = -3469342247847828789552088659323852541399766785760491146870005891371501266319724897592306597338057 / 209191710
B(110) = 7645992940484742892248134246724347500528752413412307906683593870759797606269585779977930217515 / 1518
B(112) = -2650879602155099713352597214685162014443151499192509896451788427680966756514875515366781203552600109 / 1671270
B(114) = 21737832319369163333310761086652991475721156679090831360806110114933605484234593650904188618562649 / 42
B(116) = -309553916571842976912513458033841416869004128064329844245504045721008957524571968271388199595754752259 / 1770
B(118) = 366963119969713111534947151585585006684606361080699204301059440676414485045806461889371776354517095799 / 6
B(120) = -51507486535079109061843996857849983274095170353262675213092869167199297474922985358811329367077682677803282070131 / 2328255930
B(122) = 49633666079262581912532637475990757438722790311060139770309311793150683214100431329033113678098037968564431 / 6
B(124) = -95876775334247128750774903107542444620578830013297336819553512729358593354435944413631943610268472689094609001 / 30
B(126) = 5556330281949274850616324408918951380525567307126747246796782304333594286400508981287241419934529638692081513802696639 / 4357878
</pre>
=={{header|C++}}==
{{Works with|C++11}}
{{libheader|boost}}
<syntaxhighlight lang="cpp">/**
* Configured with: --prefix=/Library/Developer/CommandLineTools/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
* Apple LLVM version 9.1.0 (clang-902.0.39.1)
* Target: x86_64-apple-darwin17.5.0
* Thread model: posix
*/
#include <boost/multiprecision/cpp_int.hpp> // 1024bit precision
#include <boost/rational.hpp> // Rationals
#include <iostream> // formatting with std::cout
#include <vector> // Container
typedef boost::rational<boost::multiprecision::int1024_t> rational; // reduce boilerplate
rational bernoulli(size_t n) {
auto out = std::vector<rational>();
for (size_t m = 0; m <= n; m++) {
out.emplace_back(1, (m + 1)); // automatically constructs object
for (size_t j = m; j >= 1; j--) {
out[j - 1] = rational(j) * (out[j - 1] - out[j]);
}
}
return out[0];
}
int main() {
for (size_t n = 0; n <= 60; n += n >= 2 ? 2 : 1) {
auto b = bernoulli(n);
std::cout << "B(" << std::right << std::setw(2) << n << ") = ";
std::cout << std::right << std::setw(44) << b.numerator();
std::cout << " / " << b.denominator() << std::endl;
}
return 0;
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Clojure}}==
<syntaxhighlight lang="clojure">
ns test-project-intellij.core
(:gen-class))
(defn a-t [n]
" Used Akiyama-Tanigawa algorithm with a single loop rather than double nested loop "
" Clojure does fractional arithmetic automatically so that part is easy "
(loop [m 0
j m
A (vec (map #(/ 1 %) (range 1 (+ n 2))))] ; Prefil A(m) with 1/(m+1), for m = 1 to n
(cond ; Three way conditional allows single loop
(>= j 1) (recur m (dec j) (assoc A (dec j) (* j (- (nth A (dec j)) (nth A j))))) ; A[j-1] ← j×(A[j-1] - A[j]) ;
(< m n) (recur (inc m) (inc m) A) ; increment m, reset j = m
:else (nth A 0))))
(defn format-ans [ans]
" Formats answer so that '/' is aligned for all answers "
(if (= ans 1)
(format "%50d / %8d" 1 1)
(format "%50d / %8d" (numerator ans) (denominator ans))))
;; Generate a set of results for [0 1 2 4 ... 60]
(doseq [q (flatten [0 1 (range 2 62 2)])
:let [ans (a-t q)]]
(println q ":" (format-ans ans)))
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Common Lisp}}==
An implementation of the simple algorithm.
Line 306 ⟶ 1,187:
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 [http://oeis.org/A027641 the numerator list].
<
(loop with a = (make-array (list (1+ n)))
for m from 0 to n do
Line 348 ⟶ 1,229:
n
(numerator r)
(denominator r)))))</
{{out}}
Line 385 ⟶ 1,266:
B(60): -1215233140483755572040304994079820246041491/56786730
</pre>
=={{header|Crystal}}==
{{Trans|Ruby}}
<syntaxhighlight lang="ruby">require "big"
class Bernoulli
include Iterator(Tuple(Int32, BigRational))
def initialize
@a = [] of BigRational
@m = 0
end
def next
@a << BigRational.new(1, @m+1)
@m.downto(1) { |j| @a[j-1] = j*(@a[j-1] - @a[j]) }
v = @m.odd? && @m != 1 ? BigRational.new(0, 1) : @a.first
return {@m, v}
ensure
@m += 1
end
end
b = Bernoulli.new
bn = b.first(61).to_a
max_width = bn.map { |_, v| v.numerator.to_s.size }.max
bn.reject { |i, v| v.zero? }.each do |i, v|
puts "B(%2i) = %*i/%i" % [i, max_width, v.numerator, v.denominator]
end
</syntaxhighlight>
{{Trans|Python}}
Version 1: compute each number separately.
<syntaxhighlight lang="ruby">require "big"
def bernoulli(n)
ar = [] of BigRational
(0..n).each do |m|
ar << BigRational.new(1, m+1)
m.downto(1) { |j| ar[j-1] = j * (ar[j-1] - ar[j]) }
end
ar[0] # (which is Bn)
end
b_nums = (0..61).map { |i| bernoulli(i) }
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? }
</syntaxhighlight>
{{Trans|Python}}
Version 2: create faster generator to compute array of numbers once.
<syntaxhighlight lang="ruby">require "big"
def bernoulli2(limit)
ar = [] of BigRational
(0..limit).each do |m|
ar << BigRational.new(1, m+1)
m.downto(1) { |j| ar[j-1] = j * (ar[j-1] - ar[j]) }
yield ar[0] # use Bn value in required block
end
end
b_nums = [] of BigRational
bernoulli2(61){ |b| b_nums << b }
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? }
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|D}}==
This uses the D module from the Arithmetic/Rational task.
{{trans|Python}}
<
auto bernoulli(in uint n) pure nothrow /*@safe*/ {
Line 406 ⟶ 1,391:
foreach (immutable b; berns)
writefln("B(%2d) = %*d/%d", b[0], width, b[1].tupleof);
}</
The output is exactly the same as the Python entry.
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{libheader| Velthuis.BigRationals}}
{{Trans|Go}}
Thanks Rudy Velthuis for the [https://github.com/rvelthuis/DelphiBigNumbers Velthuis.BigRationals] library.<br>
<syntaxhighlight lang="delphi">
program Bernoulli_numbers;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
Velthuis.BigRationals;
function b(n: Integer): BigRational;
begin
var a: TArray<BigRational>;
SetLength(a, n + 1);
for var m := 0 to High(a) do
begin
a[m] := BigRational.Create(1, m + 1);
for var j := m downto 1 do
a[j - 1] := (a[j -
end;
Result := a[0];
end;
begin
for var n := 0 to 60 do
begin
var bb := b(n);
if bb.Numerator.BitLength > 0 then
writeln(format('B(%2d) =%45s/%s', [n, bb.Numerator.ToString, bb.Denominator.ToString]));
end;
readln;
end.</syntaxhighlight>
=={{header|EchoLisp}}==
{{improve|EchoLisp|<br> Try to show '''B<sub>1</sub>''' within the output proper as -1/2.}}
Only 'small' rationals are supported in EchoLisp, i.e numerator and demominator < 2^31. So, we create a class of 'large' rationals, supported by the bigint library, and then apply the magic formula.
<
(lib 'bigint) ;; lerge numbers
(lib 'gloops) ;; classes
Line 501 ⟶ 1,462:
(define-method div (Rational Rational) (lambda (r q)
(normalize (Rational (* r.a q.b) (* r.b q.a)))))
</syntaxhighlight>
{{Output}}
<
;; Bernoulli numbers
;; http://rosettacode.org/wiki/Bernoulli_numbers
Line 551 ⟶ 1,512:
(B 1) → 1 / 2
</syntaxhighlight>
=={{header|Elixir}}==
<syntaxhighlight lang="elixir">defmodule Bernoulli do
defmodule Rational do
import Kernel, except: [div: 2]
defstruct numerator: 0, denominator: 1
def new(numerator, denominator\\1) do
sign = if numerator * denominator < 0, do: -1, else: 1
{numerator, denominator} = {abs(numerator), abs(denominator)}
gcd = gcd(numerator, denominator)
%Rational{numerator: sign * Kernel.div(numerator, gcd),
denominator: Kernel.div(denominator, gcd)}
end
def sub(a, b) do
new(a.numerator * b.denominator - b.numerator * a.denominator,
a.denominator * b.denominator)
end
def mul(a, b) when is_integer(a) do
new(a * b.numerator, b.denominator)
end
defp gcd(a,0), do: a
defp gcd(a,b), do: gcd(b, rem(a,b))
end
def numbers(n) do
Stream.transform(0..n, {}, fn m,acc ->
acc = Tuple.append(acc, Rational.new(1,m+1))
if m>0 do
new =
Enum.reduce(m..1, acc, fn j,ar ->
put_elem(ar, j-1, Rational.mul(j, Rational.sub(elem(ar,j-1), elem(ar,j))))
end)
{[elem(new,0)], new}
else
{[elem(acc,0)], acc}
end
end) |> Enum.to_list
end
def task(n \\ 61) do
b_nums = numbers(n)
width = Enum.map(b_nums, fn b -> b.numerator |> to_string |> String.length end)
|> Enum.max
format = 'B(~2w) = ~#{width}w / ~w~n'
Enum.with_index(b_nums)
|> Enum.each(fn {b,i} ->
if b.numerator != 0, do: :io.fwrite format, [i, b.numerator, b.denominator]
end)
end
end
Bernoulli.task</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|F sharp|F#}}==
{{libheader|MathNet.Numerics.FSharp}}
<syntaxhighlight lang="fsharp">
open MathNet.Numerics
open System
open System.Collections.Generic
let calculateBernoulli n =
let ℚ(x) = BigRational.FromInt x
let A = Array.init<BigRational> (n+1) (fun x -> ℚ(x+1))
for m in [1..n] do
A.[m] <- ℚ(1) / (ℚ(m) + ℚ(1))
for j in [m..(-1)..1] do
A.[j-1] <- ℚ(j) * (A.[j-1] - A.[j])
A.[0]
[<EntryPoint>]
let main argv =
for n in [0..60] do
let bernoulliNumber = calculateBernoulli n
match bernoulliNumber.Numerator.IsZero with
| false ->
let formatedString = String.Format("B({0, 2}) = {1, 44} / {2}", n, bernoulliNumber.Numerator, bernoulliNumber.Denominator)
printfn "%s" formatedString
| true ->
printf ""
0
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Factor}}==
One could use the "bernoulli" word from the math.extras vocabulary as follows:
<syntaxhighlight lang="text">IN: scratchpad
[
0 1 1 "%2d : %d / %d\n" printf
1 -1 2 "%2d : %d / %d\n" printf
30 iota [
1 + 2 * dup bernoulli [ numerator ] [ denominator ] bi
"%2d : %d / %d\n" printf
] each
] time
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
Running time: 0.00489444 seconds</syntaxhighlight>
Alternatively a method described by Brent and Harvey (2011) in "Fast computation of Bernoulli, Tangent and Secant numbers" https://arxiv.org/pdf/1108.0286.pdf is shown.
<syntaxhighlight lang="text">:: bernoulli-numbers ( n -- )
n 1 + 0 <array> :> tab
1 1 tab set-nth
2 n [a,b] [| k |
k 1 - dup
tab nth *
k tab set-nth
] each
2 n [a,b] [| k |
k n [a,b] [| j |
j tab nth
j k - 2 + *
j 1 - tab nth
j k - * +
j tab set-nth
] each
] each
1 :> s!
1 n [a,b] [| k |
k 2 * dup
2^ dup 1 - *
k tab nth
swap / *
s * k tab set-nth
s -1 * s!
] each
0 1 1 "%2d : %d / %d\n" printf
1 -1 2 "%2d : %d / %d\n" printf
1 n [a,b] [| k |
k 2 * k tab nth
[ numerator ] [ denominator ] bi
"%2d : %d / %d\n" printf
] each
;</syntaxhighlight>
It gives the same result as the native implementation, but is slightly faster.
<syntaxhighlight lang="text">[ 30 bernoulli-numbers ] time
...
Running time: 0.004331652 seconds</syntaxhighlight>
=={{header|Fermat}}==
<syntaxhighlight lang="fermat">Func Bern(m) = Sigma<k=0,m>[Sigma<v=0,k>[(-1)^v*Bin(k,v)*(v+1)^m/(k+1)]].;
for i=0, 60 do b:=Bern(i); if b<>0 then !!(i,b) fi od;</syntaxhighlight>
{{out}}<pre>
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</pre>
=={{header|FreeBASIC}}==
{{libheader|GMP}}
<syntaxhighlight lang="freebasic">' version 08-10-2016
' compile with: fbc -s console
' uses gmp
#Include Once "gmp.bi"
#Define max 60
Dim As Long n
Dim As ZString Ptr gmp_str :gmp_str = Allocate(1000) ' 1000 char
Dim Shared As Mpq_ptr tmp, big_j
tmp = Allocate(Len(__mpq_struct)) :Mpq_init(tmp)
big_j = Allocate(Len(__mpq_struct)) :Mpq_init(big_j)
Dim Shared As Mpq_ptr a(max), b(max)
For n = 0 To max
A(n) = Allocate(Len(__mpq_struct)) :Mpq_init(A(n))
B(n) = Allocate(Len(__mpq_struct)) :Mpq_init(B(n))
Next
Function Bernoulli(n As Integer) As Mpq_ptr
Dim As Long m, j
For m = 0 To n
Mpq_set_ui(A(m), 1, m + 1)
For j = m To 1 Step - 1
Mpq_sub(tmp, A(j - 1), A(j))
Mpq_set_ui(big_j, j, 1) 'big_j = j
Mpq_mul(A(j - 1), big_j, tmp)
Next
Next
Return A(0)
End Function
' ------=< MAIN >=------
For n = 0 To max
Mpq_set(B(n), Bernoulli(n))
Mpq_get_str(gmp_str, 10, B(n))
If *gmp_str <> "0" Then
If *gmp_str = "1" Then *gmp_str = "1/1"
Print Using "B(##) = "; n;
Print Space(45 - InStr(*gmp_str, "/")); *gmp_str
End If
Next
' empty keyboard buffer
While Inkey <> "" :Wend
Print :Print "hit any key to end program"
Sleep
End</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|Frink}}==
<syntaxhighlight lang="frink">BernoulliNumber[n] :=
{
a = new array
for m = 0 to n
{
a@m = 1/(m+1)
for j = m to 1 step -1
a@(j-1) = j * (a@(j-1) - a@j)
}
return a@0
}
result = new array
for n=0 to 60
{
b = BernoulliNumber[n]
if b != 0
{
[num,den] = numeratorDenominator[b]
result.push[[n, num, "/", den]]
}
}
println[formatTable[result, "right"]]</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|FunL}}==
FunL has pre-defined function <code>B</code> in module <code>integers</code>, which is defined as:
<
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) )</
{{out}}
Line 597 ⟶ 1,984:
B(60) = -1215233140483755572040304994079820246041491/56786730
</pre>
=={{header|Fōrmulæ}}==
{{FormulaeEntry|page=https://formulae.org/?script=examples/Bernoulli_numbers}}
'''Solution.''' The following function reduces to the n-th Bernoulli number. It is a replica of the Akiyama–Tanigawa algorithm.
[[File:Fōrmulæ - Bernoulli numbers 01.png]]
'''Test case.''' Showing the Bernoulli numbers B<sub>0</sub> to B<sub>60</sub>
[[File:Fōrmulæ - Bernoulli numbers 02.png]]
[[File:Fōrmulæ - Bernoulli numbers 03.png]]
=={{header|GAP}}==
<syntaxhighlight lang="gap">for a in Filtered(List([0 .. 60], n -> [n, Bernoulli(n)]), x -> x[2] <> 0) do
Print(a, "\n");
od;
[ 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 ]</syntaxhighlight>
=={{header|Go}}==
<
import (
)
func b(n int) *big.Rat {
}
}
}
func main() {
}
}
}</
{{out}}
<pre>
Line 667 ⟶ 2,099:
B(60) = -1215233140483755572040304994079820246041491/56786730
</pre>
=={{header|Haskell}}==
====Task algorithm====
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.
<syntaxhighlight lang="haskell">import Data.Ratio
import System.Environment
main = getArgs >>=
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 =
where
berno i = 1 % (i + 1)
ulli _ [_] = []
ulli i (x:y:xs) = (i % 1) * (x - y) : ulli (i + 1) (y : xs)
{{Out}}
<pre>B(0) = 1/1
B(1) = 1/2
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(30) = 8615841276005/14322
B(
B(34) = 2577687858367/6
B(36) = -26315271553053477373/1919190
B(
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</pre>
====Derivation from Faulhaber's triangle====
<syntaxhighlight lang="haskell">import Data.Bool (bool)
import Data.Ratio (Ratio, denominator, numerator, (%))
-------------------- BERNOULLI NUMBERS -------------------
bernouillis :: Integer -> [Rational]
bernouillis =
fmap head
. tail
. scanl faulhaber []
. enumFromTo 0
faulhaber :: [Ratio Integer] -> Integer -> [Ratio Integer]
faulhaber rs n =
(:) =<< (-) 1 . sum $
zipWith ((*) . (n %)) [2 ..] rs
--------------------------- TEST -------------------------
main :: IO ()
main = do
let xs = bernouillis 60
w = length $ show (numerator (last xs))
putStrLn $
fTable
"Bernouillis from Faulhaber triangle:\n"
(show . fst)
(showRatio w . snd)
id
(filter ((0 /=) . snd) $ zip [0 ..] xs)
------------------------ FORMATTING ----------------------
fTable ::
String ->
(a -> String) ->
(b -> String) ->
(a -> b) ->
[a] ->
String
fTable s xShow fxShow f xs =
let w = maximum (length . xShow <$> xs)
in unlines $
s :
fmap
( ((<>) . rjust w ' ' . xShow)
<*> ((" -> " <>) . fxShow . f)
)
xs
showRatio :: Int -> Rational -> String
showRatio w r =
let d = denominator r
in rjust w ' ' $ show (numerator r)
<> bool [] (" / " <> show d) (1 /= d)
rjust :: Int -> a -> [a] -> [a]
rjust n c = drop . length <*> (replicate n c <>)</syntaxhighlight>
{{Out}}
<pre>Bernouillis from Faulhaber triangle:
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</pre>
=={{header|Icon}} and {{header|Unicon}}==
The following works in both languages:
<
procedure main(args)
Line 751 ⟶ 2,275:
procedure align(r,n)
return repl(" ",n-find("/",r))||r
end</
Sample run:
Line 790 ⟶ 2,314:
->
</pre>
=={{header|J}}==
'''Implementation:'''
See [[j:Essays/Bernoulli_Numbers|Bernoulli Numbers Essay]] on the J wiki.
<syntaxhighlight lang="j">B=: {.&1 %. (i. ! ])@>:@i.@x:</syntaxhighlight>
'''Task:'''
<syntaxhighlight lang="j"> 'B' ,. rplc&'r/_-'"1": (#~ 0 ~: {:"1)(i. ,. B) 61
B 0 1
B 1
B 2 1/6
B 4 -1/30
Line 834 ⟶ 2,355:
B56 -2479392929313226753685415739663229/870
B58 84483613348880041862046775994036021/354
B60 -1215233140483755572040304994079820246041491/56786730</
=={{header|Java}}==
<syntaxhighlight lang="java">import org.apache.commons.math3.fraction.BigFraction;
public class BernoulliNumbers {
public static void main(String[] args) {
for (int n = 0; n <= 60; n++) {
BigFraction b = bernouilli(n);
if (!b.equals(BigFraction.ZERO))
System.out.printf("B(%-2d) = %-1s%n", n , b);
}
}
static BigFraction bernouilli(int n) {
BigFraction[] A = new BigFraction[n + 1];
for (int m = 0; m <= n; m++) {
A[m] = new BigFraction(1, (m + 1));
for (int j = m; j >= 1; j--)
A[j - 1] = (A[j - 1].subtract(A[j])).multiply(new BigFraction(j));
}
return A[0];
}
}</syntaxhighlight>
<pre>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</pre>
=={{header|jq}}==
{{works with|jq|1.4}}
Line 846 ⟶ 2,420:
'''BigInt Stubs''':
<
# def lessOrEqual(x; y): # x <= y
# def long_add(x;y): # x+y
Line 884 ⟶ 2,458:
def long_mod(x;y):
((x|tonumber) % (y|tonumber)) | tostring;</
'''Fractions''':<
# A fraction is represented by [numerator, denominator] in reduced form, with the sign on top
Line 940 ⟶ 2,514:
| if $a == $b then ["0", "1"]
else add($a; [ ($b[0]|negate), $b[1] ] )
end ; </
'''Bernoulli Numbers''':
<
def bernoulli(n):
reduce range(0; n+1) as $m
Line 952 ⟶ 2,526:
.[$j-1] = multiply( [($j|tostring), "1"]; minus( .[$j-1] ; .[$j]) ) ))
| .[0] # (which is Bn)
;</
'''The task''':
<
| if . % 2 == 0 or . == 1 then "\(.): \(bernoulli(.) )" else empty end</
{{out}}
The following output was obtained using the previously mentioned BigInt library.
<
0: ["1","1"]
1: ["1","2"]
Line 991 ⟶ 2,565:
56: ["-2479392929313226753685415739663229","870"]
58: ["84483613348880041862046775994036021","354"]
60: ["-1215233140483755572040304994079820246041491","56786730"]</
=={{header|Julia}}==
<syntaxhighlight lang="julia">function bernoulli(n)
A = Vector{Rational{BigInt}}(undef, n + 1)
for m = 0 : n
A[m + 1] = 1 // (m + 1)
for j = m : -1 : 1
A[j] = j * (A[j] - A[j + 1])
end
end
return A[1]
end
function display(n)
B = map(bernoulli, 0 : n)
pad = mapreduce(x -> ndigits(numerator(x)) + Int(x < 0), max, B)
argdigits = ndigits(n)
for i = 0 : n
if numerator(B[i + 1]) & 1 == 1
println(
"B(", lpad(i, argdigits), ") = ",
lpad(numerator(B[i + 1]), pad), " / ", denominator(B[i + 1])
)
end
end
end
display(60)
# Alternative: Following the comment in the Perl section it is much more efficient
# to compute the list of numbers instead of one number after the other.
function BernoulliList(len)
A = Vector{Rational{BigInt}}(undef, len + 1)
B = similar(A)
for n in 0 : len
A[n + 1] = 1 // (n + 1)
for j = n : -1 : 1
A[j] = j * (A[j] - A[j + 1])
end
B[n + 1] = A[1]
end
return B
end
for (n, b) in enumerate(BernoulliList(60))
isodd(numerator(b)) && println("B($(n-1)) = $b")
end </syntaxhighlight>
Produces virtually the same output as the Python version.
=={{header|Kotlin}}==
{{trans|Java}}
{{works with|Commons Math|3.3.5}}
<syntaxhighlight lang="scala">import org.apache.commons.math3.fraction.BigFraction
object Bernoulli {
operator fun invoke(n: Int) : BigFraction {
val A = Array(n + 1, init)
for (m in 0..n)
for (j in m downTo 1)
A[j - 1] = A[j - 1].subtract(A[j]).multiply(integers[j])
return A.first()
}
val max = 60
private val init = { m: Int -> BigFraction(1, m + 1) }
private val integers = Array(max + 1, { m: Int -> BigFraction(m) } )
}
fun main(args: Array<String>) {
for (n in 0..Bernoulli.max)
if (n % 2 == 0 || n == 1)
System.out.printf("B(%-2d) = %-1s%n", n, Bernoulli(n))
}</syntaxhighlight>
{{out}}
Produces virtually the same output as the Java version.
=={{header|Lua}}==
LuaJIT version with FFI and GMP library
{{trans|C}}
{{libheader|luagmp}}
{{works with|LuaJIT|2.0-2.1}}
<syntaxhighlight lang="lua">#!/usr/bin/env luajit
local gmp = require 'gmp' ('libgmp')
local ffi = require'ffi'
local mpz, mpq = gmp.types.z, gmp.types.q
local function mpq_for(buf, op, n)
for i=0,n-1 do
op(buf[i])
end
end
local function bernoulli(rop, n)
local a=ffi.new("mpq_t[?]", n+1)
mpq_for(a, gmp.q_init, n+1)
for m=0,n do
gmp.q_set_ui(a[m],1, m+1)
for j=m,1,-1 do
gmp.q_sub(a[j-1], a[j], a[j-1])
gmp.q_set_ui(rop, j, 1)
gmp.q_mul(a[j-1], a[j-1], rop)
end
end
gmp.q_set(rop,a[0])
mpq_for(a, gmp.q_clear, n+1)
end
do --MAIN
local rop=mpq()
local n,d=mpz(),mpz()
gmp.q_init(rop)
gmp.z_inits(n, d)
local to=arg[1] and tonumber(arg[1]) or 60
local from=arg[2] and tonumber(arg[2]) or 0
if from~=0 then to,from=from,to end
for i=from,to do
bernoulli(rop, i)
if gmp.q_cmp_ui(rop, 0, 1)~=0 then
gmp.q_get_num(n, rop)
gmp.q_get_den(d, rop)
gmp.printf("B(%-2g) = %44Zd / %Zd\n", i, n, d)
end
end
gmp.z_clears(n,d)
gmp.q_clear(rop)
end</syntaxhighlight>
{{out}}
<pre>> time ./bernoulli_gmp.lua
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
./bernoulli_gmp.lua 0,02s user 0,00s system 97% cpu 0,022 total</pre>
Time compare: Python 0.591 sec, C 0.023 sec, Lua 0.022-0.025
=={{header|Maple}}==
<
{{out}}
<pre>[[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]]</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
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.
<
Do[
a[[m]] = 1/m;
Line 1,010 ⟶ 2,745:
, {m, 1, n + 1}];
]
bernoulli[60]</
{{out}}
<pre>{0,1}
Line 1,048 ⟶ 2,783:
(Note from task's author: nobody is forced to use any specific algorithm, the one shown is just a suggestion.)
<
Select[%, #[[2]] != 0 &] // TableForm</
{{out}}
<pre>0 1
Line 1,084 ⟶ 2,819:
60 -(1215233140483755572040304994079820246041491/56786730)</pre>
=={{header|Maxima}}==
Using built-in function bern
<syntaxhighlight lang="maxima">
block(makelist([sconcat("B","(",i,")","="),bern(i)],i,0,60),
sublist(%%,lambda([x],x[2]#0)),
table_form(%%))
</syntaxhighlight>
{{out}}
<pre>
matrix(
["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]
)
</pre>
=={{header|Nim}}==
<syntaxhighlight lang="nim">import bignum
import strformat
const Lim = 60
#---------------------------------------------------------------------------------------------------
proc bernoulli(n: Natural): Rat =
## Compute a Bernoulli number using Akiyama–Tanigawa algorithm.
var a = newSeq[Rat](n + 1)
for m in 0..n:
a[m] = newRat(1, m + 1)
for j in countdown(m, 1):
a[j-1] = j * (a[j] - a[j-1])
result = a[0]
#———————————————————————————————————————————————————————————————————————————————————————————————————
type Info = tuple
n: int # Number index in Bernoulli sequence.
val: Rat # Bernoulli number.
var values: seq[Info] # List of values as Info tuples.
var maxLen = -1 # Maximum length.
# First step: compute the values and prepare for display.
for n in 0..Lim:
# Compute value.
if n != 1 and (n and 1) == 1: continue # Ignore odd "n" except 1.
let b = bernoulli(n)
# Check numerator length.
let len = ($b.num).len
if len > maxLen: maxLen = len
# Store information for next step.
values.add((n, b))
# Second step: display the values with '/' aligned.
for (n, b) in values:
let s = fmt"{($b.num).alignString(maxLen, '>')} / {b.denom}"
echo fmt"{n:2}: {s}"</syntaxhighlight>
{{out}}
<pre> 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</pre>
=={{header|PARI/GP}}==
<
{{out}}
<pre>0 1
Line 1,119 ⟶ 2,976:
58 84483613348880041862046775994036021/354
60 -1215233140483755572040304994079820246041491/56786730</pre>
=={{header|Pascal|FreePascal}}==
{{libheader|BigDecimalMath}}
Tested with fpc 3.0.4
<syntaxhighlight lang="pascal">
(* Taken from the 'Ada 99' project, https://marquisdegeek.com/code_ada99 *)
program BernoulliForAda99;
uses BigDecimalMath; {library for arbitary high precision BCD numbers}
type
Fraction = object
private
numerator, denominator: BigDecimal;
public
procedure assign(n, d: Int64);
procedure subtract(rhs: Fraction);
procedure multiply(value: Int64);
procedure reduce();
procedure writeOutput();
end;
function gcd(a, b: BigDecimal):BigDecimal;
begin
if (b = 0) then begin
gcd := a;
end
else begin
gcd := gcd(b, a mod b);
end;
end;
procedure Fraction.writeOutput();
var sign : char;
begin
sign := ' ';
if (numerator<0) then sign := '-';
if (denominator<0) then sign := '-';
write(sign + BigDecimalToStr(abs(numerator)):45);
write(' / ');
write(BigDecimalToStr(abs(denominator)));
end;
procedure Fraction.assign(n, d: Int64);
begin
numerator := n;
denominator := d;
end;
procedure Fraction.subtract(rhs: Fraction);
begin
numerator := numerator * rhs.denominator;
numerator := numerator - (rhs.numerator * denominator);
denominator := denominator * rhs.denominator;
end;
procedure Fraction.multiply(value: Int64);
var
temp :BigDecimal;
begin
temp := value;
numerator := numerator * temp;
end;
procedure Fraction.reduce();
var gcdResult: BigDecimal;
begin
gcdResult := gcd(numerator, denominator);
begin
numerator := numerator div gcdResult; (* div is Int64 division *)
denominator := denominator div gcdResult; (* could also use round(d/r) *)
end;
end;
function calculateBernoulli(n: Int64) : Fraction;
var
m, j: Int64;
results: array of Fraction;
begin
setlength(results, 60) ; {largest value 60}
for m:= 0 to n do
begin
results[m].assign(1, m+1);
for j:= m downto 1 do
begin
results[j-1].subtract(results[j]);
results[j-1].multiply(j);
results[j-1].reduce();
end;
end;
calculateBernoulli := results[0];
end;
(* Main program starts here *)
var
b: Int64;
result: Fraction;
begin
writeln('Calculating Bernoulli numbers...');
writeln('B( 0) : 1 / 1');
for b:= 1 to 60 do
begin
if (b<3) or ((b mod 2) = 0) then begin
result := calculateBernoulli(b);
write('B(',b:2,')');
write(' : ');
result.writeOutput();
writeln;
end;
end;
end.
</syntaxhighlight>
{{out}}
<pre>
Calculating Bernoulli numbers...
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
</pre>
=={{header|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.
Line 1,125 ⟶ 3,145:
Instead of doing the same calculations over and over again, I retain the A array until the final Bernoulli number is produced.
<
use strict;
use warnings;
Line 1,149 ⟶ 3,169:
bernoulli_print();
</syntaxhighlight>
The output is exactly the same as the Python entry.
We can also use modules for faster results. E.g.
{{libheader|ntheory}}
<
for my $n (0 .. 60) {
my($num,$den) = bernfrac($n);
printf "B(%2d) = %44s/%s\n", $n, $num, $den if $num != 0;
}</
with identical output. Or:
<
for my $n (0 .. 60) {
my($num,$den) = split "/", bernfrac($n);
printf("B(%2d) = %44s/%s\n", $n, $num, $den||1) if $num != 0;
}</
with the difference being that Pari chooses <math>B_1</math> = -½.
=={{header|Phix}}==
{{libheader|Phix/mpfr}}
{{trans|C}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #000000;">builtins</span><span style="color: #0000FF;">/</span><span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">bernoulli</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpq</span> <span style="color: #000000;">rop</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">mpq_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">m</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">m</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">mpq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">])</span>
<span style="color: #7060A8;">mpq_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">rop</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">mpq_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_free</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #004080;">mpq</span> <span style="color: #000000;">rop</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_init</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">60</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">bernoulli</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpq_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">mpq_get_num</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rop</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpq_get_den</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rop</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">ns</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">ds</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"B(%2d) = %44s / %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ns</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ds</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_free</span><span style="color: #0000FF;">({</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">rop</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_free</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(
B(58) = 84483613348880041862046775994036021 / 354
B(60) = -1215233140483755572040304994079820246041491 / 56786730
</pre>
=={{header|PicoLisp}}==
Brute force and method by Srinivasa Ramanujan.
<
(de fact (N)
Line 1,377 ⟶ 3,334:
(test (berno N) (berno-brute N)) )
(bye)</
=={{header|PL/I}}==
<
declare i fixed binary;
declare B complex fixed (31);
Line 1,412 ⟶ 3,368:
(3 A, column(10), F(32), 2 A);
end;
end Bern;</
The above uses GCD (see Rosetta Code) extended for 31-digit working.
Results obtained by this program are limited to the entries shown below due to the restrictions imposed by storing numbers in fixed decimal (31 digits).
<pre>
B(0)= 1/1
Line 1,437 ⟶ 3,394:
B(36)= -26315271553053477373/1919190
</pre>
=={{header|Python}}==
===Python: Using task algorithm===
<
def bernoulli(n):
Line 1,454 ⟶ 3,410:
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))</
{{out}}
Line 1,492 ⟶ 3,448:
===Python: Optimised task algorithm===
Using the optimization mentioned in the Perl entry to reduce intermediate calculations we create and use the generator bernoulli2():
<
A, m = [], 0
while True:
Line 1,505 ⟶ 3,461:
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))</
Output is exactly the same as before.
=={{header|Quackery}}==
<syntaxhighlight lang="quackery"> $ "bigrat.qky" loadfile
[ 1+
' [ [] ] over of swap
times
[ i^ 1+ n->v 1/v
join swap i^ poke
i^ times
[ dup i 1+ peek do
dip over swap i peek do
v- i 1+ n->v v*
join swap i poke ] ]
1 split drop do ] is bernoulli ( n --> n/d )
61 times
[ i^ bernoulli
2dup v0= iff
2drop
else
[ i^ 10 < if sp
i^ echo sp
vulgar$
char / over find
44 swap - times sp
echo$ cr ] ]</syntaxhighlight>
{{out}}
<pre> 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
</pre>
=={{header|R}}==
{{incorrect|rsplus|This example is incorrect: It is not executable and if made executable (with 'library(gmp)') it returns completely different and wrong results -- not the ones shown here. The R code needs complete rewrite and the 'pracma' library will not be of any help.}}
<syntaxhighlight lang="rsplus">
library(pracma)
for (idx in c(1,2*0:30)) {
b <- bernoulli(idx)
d <- as.character(denominator(b))
n <- as.character(numerator(b))
cat("B(",idx,") = ",n,"/",d,"\n", sep = "")
}
</syntaxhighlight>
{{out}}
<pre>
B(1) = 1/2
B(0) = 1/1
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
</pre>
=={{header|Racket}}==
Line 1,515 ⟶ 3,582:
use the same emmitter... it's just a matter of how long to wait for the emission.
<syntaxhighlight lang="text">#lang racket
;; For: http://rosettacode.org/wiki/Bernoulli_numbers
Line 1,598 ⟶ 3,665:
(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))))</
{{out}}
Line 1,633 ⟶ 3,700:
B(58) = 84483613348880041862046775994036021/354
B(60) = -1215233140483755572040304994079820246041491/56786730</pre>
=={{header|Raku}}==
(formerly Perl 6)
===Simple===
First, a straighforward implementation of the naïve algorithm in the task description.
{{works with|Rakudo|2015.12}}
<syntaxhighlight lang="raku" line>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;</syntaxhighlight>
{{out}}
<pre>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</pre>
===With memoization===
Here is a much faster way, following the Perl solution that avoids recalculating previous values each time through the function. We do this in Raku 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 (<tt>...</tt>) 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.
{{works with|Rakudo|2015.12}}
<syntaxhighlight lang="raku" line>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;</syntaxhighlight>
{{out}}
<pre>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</pre>
===Functional===
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 <tt>=></tt> 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.
{{works with|Rakudo|2016.12}}
<syntaxhighlight lang="raku" line>sub infix:<bop>(\prev, \this) {
this.key => this.key * (this.value - prev.value)
}
sub next-bernoulli ( (:key($pm), :value(@pa)) ) {
$pm + 1 => [
map *.value,
[\bop] ($pm + 2 ... 1) Z=> FatRat.new(1, $pm + 2), |@pa
]
}
constant bernoulli =
grep *.value,
map { .key => .value[*-1] },
(0 => [FatRat.new(1,1)], &next-bernoulli ... *)
;
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;</syntaxhighlight>
Same output as memoization example
=={{header|REXX}}==
The double sum formula used is number '''(33)''' from the entry [http://mathworld.wolfram.com/BernoulliNumber.html Bernoulli number] on Wolfram MathWorld<sup>TM</sup>.
<br><br>
::::::: <big><big> <math> B_n = \sum_{k=0}^n \frac{1}{k+1} \sum_{r=0}^k (-1)^r \binom kr r^n </math> </big></big>
<br>
<
parse arg N .; if N=='' | N=="," then N=
numeric digits max(9, n*2) /*increase the decimal digits if needed*/
w= max(length(N), 4); Nw= N + w + N % 4 /*used for aligning (output) fractions.*/
say
say copies('─',w) copies("─", max(78-w,Nw+2*w)) /*display 2nd line of title, separators*/
!.= .; do #=0 to N /*process the numbers from 0 ──► N. */
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
bern: parse arg x; if x==0 then return '1/1'
if x==1 then return '-1/2'
do j=2 to x by 2; jp=
sd= 2 /* " " denominator. */
sn=
sn= -sn
sd= sd * jp
@.j= sn'/'sd
end /*j*/
return sn'/'sd
/*──────────────────────────────────────────────────────────────────────────────────────*/
comb: procedure expose !.; parse arg x,y; if x==y then return 1
if !.C.x.y\==. then return !.C.x.y /*combination computed before?*/
if x-y < y then y= x-y /*x-y < y? Then use a new Y.*/
z= perm(x, y); do j=2 for y-1; z= z % j
end /*j*/
!.C.x.y= z; return z /*assign memoization & return.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
GCD: procedure; parse arg x,y; x= abs(x)
do until y==0; parse value x//y y with y x; end; return x
/*──────────────────────────────────────────────────────────────────────────────────────*/
LCM: procedure; parse arg x,y /*X=ABS(X); Y=ABS(Y) not needed for Bernoulli #s.*/
/*IF Y==0 THEN RETURN 0 " " " " " */
$= x * y /*calculate part of the LCM here. */
do until y==0; parse value x//y y with y x
end /*until*/ /* [↑] this is a short & fast GCD*/
return $ % x /*divide the pre─calculated value.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
perm: procedure expose !.; parse arg x,y; if !.P.x.y\==. then return !.P.x.y
z= 1; do j=x-y+1 to x; z= z*j; end; !.P.x.y= z; return z</syntaxhighlight>
{{out|output|text= when using the default input:}}
<pre>
B(n) Bernoulli numbers expressed as vulgar fractions
──── ───────────────────────────────────────────────────────────────────────────────────────
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
</pre>
Output notes: This version of REXX can compute and display all values up to '''B<sub>110</sub>''' in sub─second.
<br><br>
=={{header|RPL}}==
Fractions such as a/b are here handled through the complex number data structure <code>(a,b)</code>.
2 local words support the algorithm suggested by the task: <code>FracSub</code> substracts 2 fractions and <code>FracSimp</code> make a fraction irreducible
Unfortunately, floating-point precision prevents from going beyond B22.
{{works with|Halcyon Calc|4.2.7}}
{| class="wikitable"
! RPL code
! Comment
|-
|
≪ DUP2 IM * ROT IM ROT RE * -
≫ '<span style="color:blue">FracSub</span>' STO
≪ DUP C→R ABS SWAP ABS DUP2 < ≪ SWAP ≫ IFT
'''WHILE''' DUP '''REPEAT''' SWAP OVER MOD '''END'''
DROP /
≫ '<span style="color:blue">FracSimp</span>' STO
≪
{ } 1 ROT 1 + '''FOR''' m
1 m R→C +
'''IF''' m 2 ≥ '''THEN''' m 2 '''FOR''' j
DUP j 1 - GET OVER j GET <span style="color:blue">FracSub</span>
C→R SWAP j 1 - * SWAP R→C
<span style="color:blue">FracSimp</span> j 1 - SWAP PUT
-1 '''STEP END'''
'''NEXT''' 1 GET DUP RE SWAP 0 '''IFTE'''
≫ '<span style="color:blue">BRNOU</span>' STO
|
''( (a,b) (c,d) -- (e,f) )'' with e/f = a/b - c/d
''( (a,b) -- (c,d) )'' with c/d simplified fraction of a/b
get GCD of a and b
divide (a,b) by GCD
''( n -- (a,b) )'' with a/b = B(n)
For m from 1 to n+1 do
A[m] ← 1/m
for j from m to 2 do
(A[j-1] - A[j])
. * (j-1)
A[j-1] ←
return A[1] as a fraction or zero
|}
5 <span style="color:blue">BRNOU</span>
22 <span style="color:blue">BRNOU</span>
{{out}}
<pre>
2: 0
1: (854513,138)
</pre>
===HP-49+ version===
Latest RPL implementations can natively handle long fractions and generate Bernoulli numbers.
{{works with|HP|49}}
≪ { }
0 ROT '''FOR''' n
'''IF''' n 2 > LASTARG MOD AND NOT '''THEN''' n IBERNOULLI + '''END'''
'''NEXT'''
≫ '<span style="color:blue">TASK</span>' STO
60 <span style="color:blue">TASK</span>
{{out}}
<pre>
1: {1 -1/2 1/6 -1/30 1/42 -1/30 5/66 -691/2730 7/6 -3617/510 43867/798 -174611/330 854513/138 -236364091/2730 8553103/6 -23749461029/870 8615841276005/14322 -7709321041217/510 2577687858367/6 -26315271553053477373/1919190 2929993913841559/6 -261082718496449122051/13530 1520097643918070802691/1806 -27833269579301024235023/690 596451111593912163277961/282 -5609403368997817686249127547/46410 495057205241079648212477525/66 -801165718135489957347924991853/1590 9149963634884862421418123812691/798 -2479392929313226753685415739663229/870 84483613348880041862046775994036021/354 -1215233140483755572040304994079820246041491/56786730}
</pre>
Runs in 3 minutes 40 on a HP-50g, against 1 hour and 30 minutes if calculating Bernoulli numbers with the above function.
=={{header|Ruby}}==
{{trans|Python}}
<
ar
ar << Rational(1, m+1)
m.downto(1){|j| ar[j-1] = j*(ar[j-1] - ar[j]) }
y << ar.first # yield
end
end
Line 1,747 ⟶ 4,054:
b_nums.each_with_index {|b,i| puts "B(%2i) = %*i/%i" % [i, width, b.numerator, b.denominator] unless b.zero? }
</syntaxhighlight>
{{out}}
<pre>
Line 1,784 ⟶ 4,091:
</pre>
=={{header|Rust}}==
<syntaxhighlight lang="rust">// 2.5 implementations presented here: naive, optimized, and an iterator using
// the optimized function. The speeds vary significantly: relative
// speeds of optimized:iterator:naive implementations is 625:25:1.
#![feature(test)]
extern crate num;
extern crate test;
use num::bigint::{BigInt, ToBigInt};
use num::rational::{BigRational};
use std::cmp::max;
use std::env;
use std::ops::{Mul, Sub};
use std::process;
struct Bn {
value: BigRational,
index: i32
}
struct Context {
bigone_const: BigInt,
a: Vec<BigRational>,
index: i32 // Counter for iterator implementation
}
impl Context {
pub fn new() -> Context {
let bigone = 1.to_bigint().unwrap();
let a_vec: Vec<BigRational> = vec![];
Context {
bigone_const: bigone,
a: a_vec,
index: -1
}
}
}
impl Iterator for Context {
type Item = Bn;
fn next(&mut self) -> Option<Bn> {
self.index += 1;
Some(Bn { value: bernoulli(self.index as usize, self), index: self.index })
}
}
fn help() {
println!("Usage: bernoulli_numbers <up_to>");
}
fn main() {
let args: Vec<String> = env::args().collect();
let mut up_to: usize = 60;
match args.len() {
1 => {},
2 => {
up_to = args[1].parse::<usize>().unwrap();
},
_ => {
help();
process::exit(0);
}
}
let context = Context::new();
// Collect the solutions by using the Context iterator
// (this is not as fast as calling the optimized function directly).
let res = context.take(up_to + 1).collect::<Vec<_>>();
let width = res.iter().fold(0, |a, r| max(a, r.value.numer().to_string().len()));
for r in res.iter().filter(|r| *r.value.numer() != ToBigInt::to_bigint(&0).unwrap()) {
println!("B({:>2}) = {:>2$} / {denom}", r.index, r.value.numer(), width,
denom = r.value.denom());
}
}
// Implementation with no reused calculations.
fn _bernoulli_naive(n: usize, c: &mut Context) -> BigRational {
for m in 0..n + 1 {
c.a.push(BigRational::new(c.bigone_const.clone(), (m + 1).to_bigint().unwrap()));
for j in (1..m + 1).rev() {
c.a[j - 1] = (c.a[j - 1].clone().sub(c.a[j].clone())).mul(
BigRational::new(j.to_bigint().unwrap(), c.bigone_const.clone())
);
}
}
c.a[0].reduced()
}
// Implementation with reused calculations (does not require sequential calls).
fn bernoulli(n: usize, c: &mut Context) -> BigRational {
for i in 0..n + 1 {
if i >= c.a.len() {
c.a.push(BigRational::new(c.bigone_const.clone(), (i + 1).to_bigint().unwrap()));
for j in (1..i + 1).rev() {
c.a[j - 1] = (c.a[j - 1].clone().sub(c.a[j].clone())).mul(
BigRational::new(j.to_bigint().unwrap(), c.bigone_const.clone())
);
}
}
}
c.a[0].reduced()
}
#[cfg(test)]
mod tests {
use super::{Bn, Context, bernoulli, _bernoulli_naive};
use num::rational::{BigRational};
use std::str::FromStr;
use test::Bencher;
// [tests elided]
#[bench]
fn bench_bernoulli_naive(b: &mut Bencher) {
let mut context = Context::new();
b.iter(|| {
let mut res: Vec<Bn> = vec![];
for n in 0..30 + 1 {
let b = _bernoulli_naive(n, &mut context);
res.push(Bn { value:b.clone(), index: n as i32});
}
});
}
#[bench]
fn bench_bernoulli(b: &mut Bencher) {
let mut context = Context::new();
b.iter(|| {
let mut res: Vec<Bn> = vec![];
for n in 0..30 + 1 {
let b = bernoulli(n, &mut context);
res.push(Bn { value:b.clone(), index: n as i32});
}
});
}
#[bench]
fn bench_bernoulli_iter(b: &mut Bencher) {
b.iter(|| {
let context = Context::new();
let _res = context.take(30 + 1).collect::<Vec<_>>();
});
}
}
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Scala}}==
'''With Custom Rational Number Class'''<br/>
(code will run in Scala REPL with a cut-and-paste without need for a third-party library)
<syntaxhighlight lang="scala">/** Roll our own pared-down BigFraction class just for these Bernoulli Numbers */
case class BFraction( numerator:BigInt, denominator:BigInt ) {
require( denominator != BigInt(0), "Denominator cannot be zero" )
val gcd = numerator.gcd(denominator)
val num = numerator / gcd
val den = denominator / gcd
def unary_- = BFraction(-num, den)
def -( that:BFraction ) = that match {
case f if f.num == BigInt(0) => this
case f if f.den == this.den => BFraction(this.num - f.num, this.den)
case f => BFraction(((this.num * f.den) - (f.num * this.den)), this.den * f.den )
}
def *( that:Int ) = BFraction( num * that, den )
override def toString = num + " / " + den
}
def bernoulliB( n:Int ) : BFraction = {
val aa : Array[BFraction] = Array.ofDim(n+1)
for( m <- 0 to n ) {
aa(m) = BFraction(1,(m+1))
for( n <- m to 1 by -1 ) {
aa(n-1) = (aa(n-1) - aa(n)) * n
}
}
aa(0)
}
assert( {val b12 = bernoulliB(12); b12.num == -691 && b12.den == 2730 } )
val r = for( n <- 0 to 60; b = bernoulliB(n) if b.num != 0 ) yield (n, b)
val numeratorSize = r.map(_._2.num.toString.length).max
// Print the results
r foreach{ case (i,b) => {
val label = f"b($i)"
val num = (" " * (numeratorSize - b.num.toString.length)) + b.num
println( f"$label%-6s $num / ${b.den}" )
}}
</syntaxhighlight>
{{out}}
<pre>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
</pre>
=={{header|Scheme}}==
{{works with|Chez Scheme}}
<syntaxhighlight lang="scheme">; Return the n'th Bernoulli number.
(define bernoulli
(lambda (n)
(let ((a (make-vector (1+ n))))
(do ((m 0 (1+ m)))
((> m n))
(vector-set! a m (/ 1 (1+ m)))
(do ((j m (1- j)))
((< j 1))
(vector-set! a (1- j) (* j (- (vector-ref a (1- j)) (vector-ref a j))))))
(vector-ref a 0))))
; Convert a rational to a string. If an integer, ends with "/1".
(define rational->string
(lambda (rational)
(format "~a/~a" (numerator rational) (denominator rational))))
; Returns the string length of the numerator of a rational.
(define rational-numerator-length
(lambda (rational)
(string-length (format "~a" (numerator rational)))))
; Formats a rational with left-padding such that total length to the slash is as given.
(define rational-padded
(lambda (rational total-length-to-slash)
(let* ((length-padding (- total-length-to-slash (rational-numerator-length rational)))
(padding-string (make-string length-padding #\ )))
(string-append padding-string (rational->string rational)))))
; Return the Bernoulli numbers 0 through n in a list.
(define make-bernoulli-list
(lambda (n)
(if (= n 0)
(list (bernoulli n))
(append (make-bernoulli-list (1- n)) (list (bernoulli n))))))
; Print the non-zero Bernoulli numbers 0 through 60 aligning the slashes.
(let* ((bernoullis-list (make-bernoulli-list 60))
(numerator-lengths (map rational-numerator-length bernoullis-list))
(max-numerator-length (apply max numerator-lengths)))
(let print-bernoulli ((index 0) (numbers bernoullis-list))
(cond
((null? numbers))
((= 0 (car numbers))
(print-bernoulli (1+ index) (cdr numbers)))
(else
(printf "B(~2@a) = ~a~%" index (rational-padded (car numbers) max-numerator-length))
(print-bernoulli (1+ index) (cdr numbers))))))</syntaxhighlight>
{{out}}
<pre>$ scheme --script bernoulli.scm
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</pre>
=={{header|Seed7}}==
The program below uses [http://seed7.sourceforge.net/manual/types.htm#bigRational bigRational]
numbers. The Bernoulli numbers are written as fraction and as decimal number, with possible repeating decimals.
The conversion of a bigRational number to [http://seed7.sourceforge.net/manual/types.htm#string string] is done
with the function [http://seed7.sourceforge.net/libraries/bigrat.htm#str(in_bigRational) str]. This
function automatically writes repeating decimals in parentheses, when necessary.
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i";
include "bigrat.s7i";
const func bigRational: bernoulli (in integer: n) is func
result
var bigRational: bernoulli is bigRational.value;
local
var integer: m is 0;
var integer: j is 0;
var array bigRational: a is 0 times bigRational.value;
begin
a := [0 .. n] times bigRational.value;
for m range 0 to n do
a[m] := 1_ / bigInteger(succ(m));
for j range m downto 1 do
a[pred(j)] := bigRational(j) * (a[j] - a[pred(j)]);
end for;
end for;
bernoulli := a[0];
end func;
const proc: main is func
local
var bigRational: bernoulli is bigRational.value;
var integer: i is 0;
begin
for i range 0 to 60 do
bernoulli := bernoulli(i);
if bernoulli <> bigRational.value then
writeln("B(" <& i lpad 2 <& ") = " <& bernoulli.numerator lpad 44 <&
" / " <& bernoulli.denominator rpad 8 <& " " <& bernoulli);
end if;
end for;
end func;</syntaxhighlight>
{{out}}
<pre>
B( 0) = 1 / 1 1.0
B( 1) = -1 / 2 -0.5
B( 2) = 1 / 6 0.1(6)
B( 4) = -1 / 30 -0.0(3)
B( 6) = 1 / 42 0.0(238095)
B( 8) = -1 / 30 -0.0(3)
B(10) = 5 / 66 0.0(75)
B(12) = -691 / 2730 -0.2(531135)
B(14) = 7 / 6 1.1(6)
B(16) = -3617 / 510 -7.0(9215686274509803)
B(18) = 43867 / 798 54.9(711779448621553884)
B(20) = -174611 / 330 -529.1(24)
B(22) = 854513 / 138 6192.1(2318840579710144927536)
B(24) = -236364091 / 2730 -86580.2(531135)
B(26) = 8553103 / 6 1425517.1(6)
B(28) = -23749461029 / 870 -27298231.0(6781609195402298850574712643)
B(30) = 8615841276005 / 14322 601580873.9(006423683843038681748359167714)
B(32) = -7709321041217 / 510 -15116315767.0(9215686274509803)
B(34) = 2577687858367 / 6 429614643061.1(6)
B(36) = -26315271553053477373 / 1919190 -13711655205088.3(327721590879485616)
B(38) = 2929993913841559 / 6 488332318973593.1(6)
B(40) = -261082718496449122051 / 13530 -19296579341940068.1(4863266814)
B(42) = 1520097643918070802691 / 1806 841693047573682615.0(005537098560354374307862679955703211517165)
B(44) = -27833269579301024235023 / 690 -40338071854059455413.0(7681159420289855072463)
B(46) = 596451111593912163277961 / 282 2115074863808199160560.1(4539007092198581560283687943262411347517730496)
B(48) = -5609403368997817686249127547 / 46410 -120866265222965259346027.3(119370825253178194354664942900237017884076707606)
B(50) = 495057205241079648212477525 / 66 7500866746076964366855720.0(75)
B(52) = -801165718135489957347924991853 / 1590 -503877810148106891413789303.0(5220125786163)
B(54) = 29149963634884862421418123812691 / 798 36528776484818123335110430842.9(711779448621553884)
B(56) = -2479392929313226753685415739663229 / 870 -2849876930245088222626914643291.0(6781609195402298850574712643)
B(58) = 84483613348880041862046775994036021 / 354 238654274996836276446459819192192.1(4971751412429378531073446327683615819209039548022598870056)
B(60) = -1215233140483755572040304994079820246041491 / 56786730 -21399949257225333665810744765191097.3(926741511617238745742183076926598872659158222352299560126106)
</pre>
=={{header|Sidef}}==
Built-in:
<syntaxhighlight lang="ruby">say bernoulli(42).as_frac #=> 1520097643918070802691/1806</syntaxhighlight>
Recursive solution (with auto-memoization):
<syntaxhighlight lang="ruby">func bernoulli_number(n) is cached {
n.is_one && return
n.is_odd && return 0
1 - sum(^n,
binomial(n,k) * __FUNC__(k) / (n - k + 1)
})
}
for n in (0..60) {
var Bn = bernoulli_number(n) || next
printf("B(%2d) = %44s / %s\n", n, Bn.nude)
}</syntaxhighlight>
Using Ramanujan's congruences (pretty fast):
<syntaxhighlight lang="ruby">func ramanujan_bernoulli_number(n) is cached {
return 1/2 if n.is_one
return 0 if n.is_odd
((n%6 == 4 ? -1/2 : 1) * (n+3)/3 - sum(1 .. (n - n%6)/6, {|k|
binomial(n+3, n - 6*k) * __FUNC__(n - 6*k)
})) / binomial(n+3, n)
}</syntaxhighlight>
Using Euler's product formula for the Riemann zeta function and the Von Staudt–Clausen theorem (very fast):
<syntaxhighlight lang="ruby">func bernoulli_number_from_zeta(n) {
n.is_zero && return 1
n.is_one && return 1/2
n.is_odd && return 0
var log2B = (log(4*Num.tau*n)/2 + n*log(n) - n*log(Num.tau) - n)/log(2)
local Num!PREC = *(int(n + log2B) + (n <= 90 ? 18 : 0))
var K = 2*(n! / Num.tau**n)
var d = n.divisors.grep {|k| is_prime(k+1) }.prod {|k| k+1 }
var z = ceil((K*d).root(n-1)).primes.prod {|p| 1 - p.float**(-n) }
(-1)**(n/2 + 1) * int(ceil(d*K / z)) / d
}</syntaxhighlight>
The Akiyama–Tanigawa algorithm:
<syntaxhighlight lang="ruby">func bernoulli_print {
var a = []
for m in (0..60) {
a << 1/(m+1)
for j in (1..m -> flip) {
(a[j-1] -= a[j]) *= j
}
a[0] || next
printf("B(%2d) = %44s / %s\n", m, a[0].nude)
}
}
bernoulli_print()</syntaxhighlight>
{{out}}
<pre>
B( 0) = 1
B( 1) = 1 /
B( 2) = 1 /
B( 4) =
B( 6) = 1 /
B( 8) =
B(10) = 5 /
B(12) =
B(14) = 7 /
B(16) =
B(18) = 43867 /
B(20) =
B(22) = 854513 /
B(24) =
B(26) = 8553103 /
B(28) =
B(30) = 8615841276005 /
B(32) =
B(34) = 2577687858367 /
B(36) =
B(38) = 2929993913841559 /
B(40) =
B(42) = 1520097643918070802691 /
B(44) =
B(46) = 596451111593912163277961 /
B(48) =
B(50) = 495057205241079648212477525 /
B(52) =
B(54) = 29149963634884862421418123812691 /
B(56) =
B(58) = 84483613348880041862046775994036021 /
B(60) =
</pre>
=={{header|SPAD}}==
{{works with|FriCAS, OpenAxiom, Axiom}}
<syntaxhighlight lang="spad">
for n in 0..60 | (b:=bernoulli(n)$INTHEORY; b~=0) repeat print [n,b]
</syntaxhighlight>
Package:[http://fricas.github.io/api/IntegerNumberTheoryFunctions.html?highlight=bernoulli IntegerNumberTheoryFunctions]
{{out}}
<pre>
===============
Format: [n,B_n]
===============
[0,1]
1
[1,- -]
2
1
[2,-]
6
1
[4,- --]
30
1
[6,--]
42
1
[8,- --]
30
5
[10,--]
66
691
[12,- ----]
2730
7
[14,-]
6
3617
[16,- ----]
510
43867
[18,-----]
798
174611
[20,- ------]
330
854513
[22,------]
138
236364091
[24,- ---------]
2730
8553103
[26,-------]
6
23749461029
[28,- -----------]
870
8615841276005
[30,-------------]
14322
7709321041217
[32,- -------------]
510
2577687858367
[34,-------------]
6
26315271553053477373
[36,- --------------------]
1919190
2929993913841559
[38,----------------]
6
261082718496449122051
[40,- ---------------------]
13530
1520097643918070802691
[42,----------------------]
1806
27833269579301024235023
[44,- -----------------------]
690
596451111593912163277961
[46,------------------------]
282
5609403368997817686249127547
[48,- ----------------------------]
46410
495057205241079648212477525
[50,---------------------------]
66
801165718135489957347924991853
[52,- ------------------------------]
1590
29149963634884862421418123812691
[54,--------------------------------]
798
2479392929313226753685415739663229
[56,- ----------------------------------]
870
84483613348880041862046775994036021
[58,-----------------------------------]
354
1215233140483755572040304994079820246041491
[60,- -------------------------------------------]
56786730
Type: Void
</pre>
=={{header|Swift}}==
{{libheader|AttaSwift BigInt}}
Uses the Frac type defined in the [http://rosettacode.org/wiki/Arithmetic/Rational#Swift Rational] task.
<syntaxhighlight lang="swift">import BigInt
public func bernoulli<T: BinaryInteger & SignedNumeric>(n: Int) -> Frac<T> {
guard n != 0 else {
return 1
}
var arr = [Frac<T>]()
for m in 0...n {
arr.append(Frac(numerator: 1, denominator: T(m) + 1))
for j in stride(from: m, through: 1, by: -1) {
arr[j-1] = (arr[j-1] - arr[j]) * Frac(numerator: T(j), denominator: 1)
}
}
return arr[0]
}
for n in 0...60 {
let b = bernoulli(n: n) as Frac<BigInt>
guard b != 0 else {
continue
}
print("B(\(n)) = \(b)")
}</syntaxhighlight>
{{out}}
<pre>B(0) = Frac(1 / 1)
B(1) = Frac(1 / 2)
B(2) = Frac(1 / 6)
B(4) = Frac(-1 / 30)
B(6) = Frac(1 / 42)
B(8) = Frac(-1 / 30)
B(10) = Frac(5 / 66)
B(12) = Frac(-691 / 2730)
B(14) = Frac(7 / 6)
B(16) = Frac(-3617 / 510)
B(18) = Frac(43867 / 798)
B(20) = Frac(-174611 / 330)
B(22) = Frac(854513 / 138)
B(24) = Frac(-236364091 / 2730)
B(26) = Frac(8553103 / 6)
B(28) = Frac(-23749461029 / 870)
B(30) = Frac(8615841276005 / 14322)
B(32) = Frac(-7709321041217 / 510)
B(34) = Frac(2577687858367 / 6)
B(36) = Frac(-26315271553053477373 / 1919190)
B(38) = Frac(2929993913841559 / 6)
B(40) = Frac(-261082718496449122051 / 13530)
B(42) = Frac(1520097643918070802691 / 1806)
B(44) = Frac(-27833269579301024235023 / 690)
B(46) = Frac(596451111593912163277961 / 282)
B(48) = Frac(-5609403368997817686249127547 / 46410)
B(50) = Frac(495057205241079648212477525 / 66)
B(52) = Frac(-801165718135489957347924991853 / 1590)
B(54) = Frac(29149963634884862421418123812691 / 798)
B(56) = Frac(-2479392929313226753685415739663229 / 870)
B(58) = Frac(84483613348880041862046775994036021 / 354)
B(60) = Frac(-1215233140483755572040304994079820246041491 / 56786730)</pre>
=={{header|Tcl}}==
<
for {set m 0} {$m <= $n} {incr m} {
lappend A [list 1 [expr {$m + 1}]]
Line 1,862 ⟶ 4,834:
foreach {n num denom} $result {
puts [format {B_%-2d = %*lld/%lld} $n $len $num $denom]
}</
{{out}}
<pre>
Line 1,897 ⟶ 4,869:
B_58 = 84483613348880041862046775994036021/354
B_60 = -1215233140483755572040304994079820246041491/56786730
</pre>
=={{header|Visual Basic .NET}}==
{{works with|Visual Basic .NET|2013}}
{{libheader|System.Numerics}}
<syntaxhighlight lang="vbnet">' Bernoulli numbers - vb.net - 06/03/2017
Imports System.Numerics 'BigInteger
Module Bernoulli_numbers
Function gcd_BigInt(ByVal x As BigInteger, ByVal y As BigInteger) As BigInteger
Dim y2 As BigInteger
x = BigInteger.Abs(x)
Do
y2 = BigInteger.Remainder(x, y)
x = y
y = y2
Loop Until y = 0
Return x
End Function 'gcd_BigInt
Sub bernoul_BigInt(n As Integer, ByRef bnum As BigInteger, ByRef bden As BigInteger)
Dim j, m As Integer
Dim f As BigInteger
Dim anum(), aden() As BigInteger
ReDim anum(n + 1), aden(n + 1)
For m = 0 To n
anum(m + 1) = 1
aden(m + 1) = m + 1
For j = m To 1 Step -1
anum(j) = j * (aden(j + 1) * anum(j) - aden(j) * anum(j + 1))
aden(j) = aden(j) * aden(j + 1)
f = gcd_BigInt(BigInteger.Abs(anum(j)), BigInteger.Abs(aden(j)))
If f <> 1 Then
anum(j) = anum(j) / f
aden(j) = aden(j) / f
End If
Next
Next
bnum = anum(1) : bden = aden(1)
End Sub 'bernoul_BigInt
Sub bernoulli_BigInt()
Dim i As Integer
Dim bnum, bden As BigInteger
bnum = 0 : bden = 0
For i = 0 To 60
bernoul_BigInt(i, bnum, bden)
If bnum <> 0 Then
Console.WriteLine("B(" & i & ")=" & bnum.ToString("D") & "/" & bden.ToString("D"))
End If
Next i
End Sub 'bernoulli_BigInt
End Module 'Bernoulli_numbers</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Wren}}==
{{libheader|Wren-fmt}}
{{libheader|Wren-big}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
import "./big" for BigRat
var bernoulli = Fn.new { |n|
if (n < 0) Fiber.abort("Argument must be non-negative")
var a = List.filled(n+1, null)
for (m in 0..n) {
a[m] = BigRat.new(1, m+1)
var j = m
while (j >= 1) {
a[j-1] = (a[j-1] - a[j]) * BigRat.new(j, 1)
j = j - 1
}
}
return (n != 1) ? a[0] : -a[0] // 'first' Bernoulli number
}
for (n in 0..60) {
var b = bernoulli.call(n)
if (b != BigRat.zero) Fmt.print("B($2d) = $44i / $i", n, b.num, b.den)
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|zkl}}==
{{trans|EchoLisp}}
Uses lib GMP (GNU MP Bignum Library).
<syntaxhighlight lang="zkl">class Rational{ // Weenie Rational class, can handle BigInts
fcn init(_a,_b){ var a=_a, b=_b; normalize(); }
fcn toString{ "%50d / %d".fmt(a,b) }
fcn normalize{ // divide a and b by gcd
g:= a.gcd(b);
a/=g; b/=g;
if(b<0){ a=-a; b=-b; } // denominator > 0
self
}
fcn __opAdd(n){
if(Rational.isChildOf(n)) self(a*n.b + b*n.a, b*n.b); // Rat + Rat
else self(b*n + a, b); // Rat + Int
}
fcn __opSub(n){ self(a*n.b - b*n.a, b*n.b) } // Rat - Rat
fcn __opMul(n){
if(Rational.isChildOf(n)) self(a*n.a, b*n.b); // Rat * Rat
else self(a*n, b); // Rat * Int
}
fcn __opDiv(n){ self(a*n.b,b*n.a) } // Rat / Rat
}</syntaxhighlight>
<syntaxhighlight lang="zkl">var [const] BN=Import.lib("zklBigNum"); // libGMP (GNU MP Bignum Library)
fcn B(N){ // calculate Bernoulli(n)
var A=List.createLong(100,0); // aka static aka not thread safe
foreach m in (N+1){
A[m]=Rational(BN(1),BN(m+1));
foreach j in ([m..1, -1]){ A[j-1]= (A[j-1] - A[j])*j; }
}
A[0]
}</syntaxhighlight>
<syntaxhighlight lang="zkl">foreach b in ([0..1].chain([2..60,2])){ println("B(%2d)%s".fmt(b,B(b))) }</syntaxhighlight>
{{out}}
<pre>
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
</pre>
| |||