# Bernoulli Numbers

The Bernoulli numbers (specifically, the second Bernoulli numbers) are defined by the following recursive definition:

Where denotes a combination.

Given a nonnegative integer m as input, output the decimal representation OR a reduced fraction for the mth second Bernoulli number. If you output a decimal representation, you must have at least 6 decimal places (digits after the decimal point) of precision, and it must be accurate when rounded to 6 decimal places. For example, for m = 2, 0.166666523 is acceptable because it rounds to 0.166667. 0.166666389 is not acceptable, because it rounds to 0.166666. Trailing zeroes may be omitted. Scientific notation may be used for decimal representations.

Here is the input and expected output for m up to and including 60, in scientific notation rounded to 6 decimal places, and as reduced fractions:

0 -> 1.000000e+00 (1/1)
1 -> 5.000000e-01 (1/2)
2 -> 1.666667e-01 (1/6)
3 -> 0.000000e+00 (0/1)
4 -> -3.333333e-02 (-1/30)
5 -> 0.000000e+00 (0/1)
6 -> 2.380952e-02 (1/42)
7 -> 0.000000e+00 (0/1)
8 -> -3.333333e-02 (-1/30)
9 -> 0.000000e+00 (0/1)
10 -> 7.575758e-02 (5/66)
11 -> 0.000000e+00 (0/1)
12 -> -2.531136e-01 (-691/2730)
13 -> 0.000000e+00 (0/1)
14 -> 1.166667e+00 (7/6)
15 -> 0.000000e+00 (0/1)
16 -> -7.092157e+00 (-3617/510)
17 -> 0.000000e+00 (0/1)
18 -> 5.497118e+01 (43867/798)
19 -> 0.000000e+00 (0/1)
20 -> -5.291242e+02 (-174611/330)
21 -> 0.000000e+00 (0/1)
22 -> 6.192123e+03 (854513/138)
23 -> 0.000000e+00 (0/1)
24 -> -8.658025e+04 (-236364091/2730)
25 -> 0.000000e+00 (0/1)
26 -> 1.425517e+06 (8553103/6)
27 -> 0.000000e+00 (0/1)
28 -> -2.729823e+07 (-23749461029/870)
29 -> 0.000000e+00 (0/1)
30 -> 6.015809e+08 (8615841276005/14322)
31 -> 0.000000e+00 (0/1)
32 -> -1.511632e+10 (-7709321041217/510)
33 -> 0.000000e+00 (0/1)
34 -> 4.296146e+11 (2577687858367/6)
35 -> 0.000000e+00 (0/1)
36 -> -1.371166e+13 (-26315271553053477373/1919190)
37 -> 0.000000e+00 (0/1)
38 -> 4.883323e+14 (2929993913841559/6)
39 -> 0.000000e+00 (0/1)
40 -> -1.929658e+16 (-261082718496449122051/13530)
41 -> 0.000000e+00 (0/1)
42 -> 8.416930e+17 (1520097643918070802691/1806)
43 -> 0.000000e+00 (0/1)
44 -> -4.033807e+19 (-27833269579301024235023/690)
45 -> 0.000000e+00 (0/1)
46 -> 2.115075e+21 (596451111593912163277961/282)
47 -> 0.000000e+00 (0/1)
48 -> -1.208663e+23 (-5609403368997817686249127547/46410)
49 -> 0.000000e+00 (0/1)
50 -> 7.500867e+24 (495057205241079648212477525/66)
51 -> 0.000000e+00 (0/1)
52 -> -5.038778e+26 (-801165718135489957347924991853/1590)
53 -> 0.000000e+00 (0/1)
54 -> 3.652878e+28 (29149963634884862421418123812691/798)
55 -> 0.000000e+00 (0/1)
56 -> -2.849877e+30 (-2479392929313226753685415739663229/870)
57 -> 0.000000e+00 (0/1)
58 -> 2.386543e+32 (84483613348880041862046775994036021/354)
59 -> 0.000000e+00 (0/1)
60 -> -2.139995e+34 (-1215233140483755572040304994079820246041491/56786730)


Reference implementation (in Python 3):

def factorial(n):
if n < 1:
return 1
else:
return n * factorial(n - 1)

def combination(m,k):
if k <= m:
return factorial(m)/(factorial(k) * factorial(m - k))
else:
return 0

def Bernoulli(m):
if m == 0:
return 1
else:
t = 0
for k in range(0, m):
t += combination(m, k) * Bernoulli(k) / (m - k + 1)
return 1 - t


### Rules

• This is , so shortest code in bytes wins
• You may not use any functions, either built-in or included in an external library, that compute either type of Bernoulli numbers or Bernoulli polynomials.
• Your answer must give correct output for all inputs up to and including 60.

The Stack Snippet at the bottom of this post generates the leaderboard from the answers a) as a list of shortest solution per language and b) as an overall leaderboard.

## Language Name, N bytes


where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:

## Ruby, <s>104</s> <s>101</s> 96 bytes


If there you want to include multiple numbers in your header (e.g. because your score is the sum of two files or you want to list interpreter flag penalties separately), make sure that the actual score is the last number in the header:

## Perl, 43 + 2 (-p flag) = 45 bytes


You can also make the language name a link which will then show up in the snippet:

## [><>](http://esolangs.org/wiki/Fish), 121 bytes


<style>body { text-align: left !important} #answer-list { padding: 10px; width: 290px; float: left; } #language-list { padding: 10px; width: 290px; float: left; } table thead { font-weight: bold; } table td { padding: 5px; }</style><script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"> <div id="language-list"> <h2>Shortest Solution by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr> </thead> <tbody id="languages"> </tbody> </table> </div> <div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr> </thead> <tbody id="answers"> </tbody> </table> </div> <table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr> </tbody> </table> <table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr> </tbody> </table><script>var QUESTION_ID = 65382; var ANSWER_FILTER = "!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe"; var COMMENT_FILTER = "!)Q2B_A2kjfAiU78X(md6BoYk"; var OVERRIDE_USER = 45941; var answers = [], answers_hash, answer_ids, answer_page = 1, more_answers = true, comment_page; function answersUrl(index) { return "https://api.stackexchange.com/2.2/questions/" + QUESTION_ID + "/answers?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + ANSWER_FILTER; } function commentUrl(index, answers) { return "https://api.stackexchange.com/2.2/answers/" + answers.join(';') + "/comments?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + COMMENT_FILTER; } function getAnswers() { jQuery.ajax({ url: answersUrl(answer_page++), method: "get", dataType: "jsonp", crossDomain: true, success: function (data) { answers.push.apply(answers, data.items); answers_hash = []; answer_ids = []; data.items.forEach(function(a) { a.comments = []; var id = +a.share_link.match(/\d+/); answer_ids.push(id); answers_hash[id] = a; }); if (!data.has_more) more_answers = false; comment_page = 1; getComments(); } }); } function getComments() { jQuery.ajax({ url: commentUrl(comment_page++, answer_ids), method: "get", dataType: "jsonp", crossDomain: true, success: function (data) { data.items.forEach(function(c) { if (c.owner.user_id === OVERRIDE_USER) answers_hash[c.post_id].comments.push(c); }); if (data.has_more) getComments(); else if (more_answers) getAnswers(); else process(); } }); } getAnswers(); var SCORE_REG = /<h\d>\s*([^\n,<]*(?:<(?:[^\n>]*>[^\n<]*<\/[^\n>]*>)[^\n,<]*)*),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/; var OVERRIDE_REG = /^Override\s*header:\s*/i; function getAuthorName(a) { return a.owner.display_name; } function process() { var valid = []; answers.forEach(function(a) { var body = a.body; a.comments.forEach(function(c) { if(OVERRIDE_REG.test(c.body)) body = '<h1>' + c.body.replace(OVERRIDE_REG, '') + '</h1>'; }); var match = body.match(SCORE_REG); if (match) valid.push({ user: getAuthorName(a), size: +match[2], language: match[1], link: a.share_link, }); else console.log(body); }); valid.sort(function (a, b) { var aB = a.size, bB = b.size; return aB - bB }); var languages = {}; var place = 1; var lastSize = null; var lastPlace = 1; valid.forEach(function (a) { if (a.size != lastSize) lastPlace = place; lastSize = a.size; ++place; var answer = jQuery("#answer-template").html(); answer = answer.replace("{{PLACE}}", lastPlace + ".") .replace("{{NAME}}", a.user) .replace("{{LANGUAGE}}", a.language) .replace("{{SIZE}}", a.size) .replace("{{LINK}}", a.link); answer = jQuery(answer); jQuery("#answers").append(answer); var lang = a.language; lang = jQuery('<a>'+lang+'</a>').text(); languages[lang] = languages[lang] || {lang: a.language, lang_raw: lang.toLowerCase(), user: a.user, size: a.size, link: a.link}; }); var langs = []; for (var lang in languages) if (languages.hasOwnProperty(lang)) langs.push(languages[lang]); langs.sort(function (a, b) { if (a.lang_raw > b.lang_raw) return 1; if (a.lang_raw < b.lang_raw) return -1; return 0; }); for (var i = 0; i < langs.length; ++i) { var language = jQuery("#language-template").html(); var lang = langs[i]; language = language.replace("{{LANGUAGE}}", lang.lang) .replace("{{NAME}}", lang.user) .replace("{{SIZE}}", lang.size) .replace("{{LINK}}", lang.link); language = jQuery(language); jQuery("#languages").append(language); } }</script>

• @MorganThrapp The reference implementation is only to make clearer the definition of a Bernoulli number, not to actually solve the problem.
– user45941
Nov 30, 2015 at 21:04
• Ah, gotcha. I thought it was a fully functional implementation. Nov 30, 2015 at 21:05
• @Mego No standard float (not even quad precision) can store B_60 to quad precision. Should we then use an extended precision format if we're outputting as decimals? Nov 30, 2015 at 22:39
• I don't like the precision requirement. Some languages don't have the tools to work with floats to sufficient accuracy for B_60, and I'd rather not deal with such issues when golfing a mathematical problem. It's frustrating to write up a solution and then find that it's invalid due to what seems like a technicality.
– xnor
Dec 1, 2015 at 3:28
• @xnor 6 digits of accuracy seems incredibly lax already. Dec 3, 2015 at 12:58

# Mathematica, 402823 22 bytes

Using the famous formula n*ζ(1−n)=−Bn, where ζ is the Riemann Zeta function.

If[#>0,-Zeta[1-#]#,1]&


The same length:

B@0=1;B@n_=-Zeta[1-n]n


Original solution, 40 bytes, using the generating function of Bernoulli numbers.

#!SeriesCoefficient[t/(1-E^-t),{t,0,#}]&

• +2 if I could... Dec 1, 2015 at 12:04

# Julia, 58 bytes

B(m)=m<1?1:1-sum(k->big(binomial(m,k))*B(k)/(m-k+1),0:m-1)


This creates a recursive function B that accepts an integer and returns a BigFloat (i.e. high precision floating point).

Ungolfed:

function B(m::Integer)
m == 0 && return 1
return 1 - sum(k -> big(binomial(m, k)) * B(k) / (m-k+1), 0:m-1)
end


## Minkolang 0.14, 97 bytes

I actually tried doing it recursively first, but my interpreter, as currently designed, actually can't do it. If you try to recurse from within a for loop, it starts a new recursion. So I went for the tabulating approach...which had precision problems. So I did the whole thing with fractions. With no built-in support for fractions. [sigh]

n1+[xxi$z0z2%1+F0c0=$1&$d4Mdm:1R:r$dz1Az0A]$:N. 11z[i0azi6M*i1azi-1+*d0c*1c2c*-1c3c*4$X]f
z1=d1+f


Try it here. Bonus: the array has all fractions for every previous Bernoulli number!

### Explanation (in a bit)

n1+                 Take number from input (N) and add 1
[                Open for loop that runs N+1 times (starts at zero)
xx              Dump the top two values of the stack
i$z Store the loop counter in the register (m) 0 Push 0 z2%1+ Push 1 if N is even, 2 if odd F Gosub; pops y,x then goes to codebox(x,y), to be returned to 0c Copy the first item on the stack , 1 if equal to 0, 0 otherwise$1&                             Jump 11 spaces if top of stack is not 0

(If top of stack is not 0, then...)
$d Duplicate whole stack 4M Pop b,a and push GCD(a,b) dm Duplicate and merge (a,b,c,c -> a,c,b,c) : Divide 1R Rotate 1 item to the right (0G works too) : Divide r Reverse stack (In both cases...)$d                Duplicate whole stack
z1A             Store denominator of B_m in array
z0A         Store numerator of B_m in array
]        Close for loop
$: Divide (float division) N. Output as number and stop. 11 Push two 1s (a, b) z[ Open a for loop that repeats m times i0a Retrieve numerator of B_k (p) zi Push m, k 6M Pop k,m and push mCk (binomial) (x) * p*x (c) i1a Retrieve denominator of B_k (q) zi-1+ m-k+1 (y) * q*y (d) d Duplicate top of stack 0c* a*d 1c2c* b*c - a*d-b*c 1c3c* b*d 4$X      Dump the bottom four items of stack
]f    Jump back to F

z          m
1=        0 (if m is 1) or 1 (otherwise)
d1+     Duplicate and add 1 (1 or 2)
f    Jump back to F


The third line is responsible for 1/2 if m is 1 and 0/1 if m is an odd number greater than 1. The second line calculates B_m with the summation formula given in the question, and does so entirely with numerators and denominators. Otherwise it'd be a lot shorter. The first half of the first line does some bookkeeping and chooses whether to execute the second or third line, and the second half divides the numerator and denominator by their GCD (if applicable) and stores those values. And outputs the answer at the end.

# Julia, 23 20 bytes

Saved 3 bytes thanks to Alex A.

It uses the same formula as my Mathematica solution and PARI/GP solution.

n->n>0?-zeta(1-n)n:1

• 20 bytes: n->n>0?-zeta(1-n)n:1 Dec 2, 2015 at 21:02
• @AlexA. I don't know why, but zeta(n) throws an error when n is a negative integer. I'm using Julia 0.2.1 on linux. Dec 3, 2015 at 3:12
• Oh my, your version of Julia is quite outdated. It works just fine for me on 0.4.1. Dec 3, 2015 at 3:15

# Python 2, 118 bytes

Saved 6 bytes due to xsot.
Saved 6 10 more due to Peter Taylor.

n=input()
a=[n%4-1,n<2]*n;exec"a=[(a[j-1]+a[j+1])*j/2for j in range(len(a)-2)];"*~-n
print+(n<1)or-n/(2.**n-4**n)*a[1]


Uses the following identity:

where An is the nth Alternating Number, which can be formally defined as the number of alternating permutations on a set of size n, halved (see also: A000111).

The algorithm used was originally given by Knuth and Buckholtz (1967):

Let T1,k = 1 for all k = 1..n

Subsequent values of T are given by the recurrence relation:

Tn+1,k = 1/2 [ (k - 1) Tn,k-1 + (k + 1) Tn,k+1 ]

An is then given by Tn,1

## Python 2, 152 bytes

from fractions import*
n=input()
a=[n%4-1,n<2]*n
for k in range(n-1):a=[(a[j-1]+a[j+1])*j/2for j in range(n-k)]
print+(n<1)or Fraction(n*a[1],4**n-2**n)


Prints the exact fractional representation, necessary for values greater than 200 or so.

• If you change range(2,n) to range(n-2) you can shorten n-k+1 to n+~k. Also, is there a reason you use >>1 instead of /2? Lastly, a trivial improvement, but you can save a few bytes by aliasing range.
– xsot
Dec 1, 2015 at 10:49
• Thanks for the suggestions. I originally had two expressions, when I joined them I overlooked changing >>1 with /2. Dec 1, 2015 at 11:26
• There's a one-char saving in the output line: print+(n<1)or-(-1.)**(n+n/2)*n/(4**n-2**n)*a[n%2^1%n]. And the calculation can be done for the same char count as a=[1]*(n+1);exec"a=[(a[j-1]+a[j+1])*j/2for j in range(len(a)-1)];"*(n-1) Dec 2, 2015 at 23:13
• @PeterTaylor the n+n/2 is clever; 1 doesn't need to be singled out, because all other odd values are zero anyway. The alternative calculation is actually 4 bytes shorter with bit inversions, but also considerably slower, for some reason. Dec 3, 2015 at 1:55
• I was working from the OEIS table, and thought that you'd found range and skipping one iteration to be a clever way to shorten the initialisation. The way you've now split out even and odd indices is very nice, and allows a further saving by pulling the sign into the definition of a: a=[(-1)**(n/2),n<2]*n. Then the return value is +(n<1)or-n/(2.**n-4**n)*a[1]. You've also got a stray semicolon at the end of line 2. Dec 3, 2015 at 8:44

# PARI/GP, 52 23 bytes

Using the famous formula n*ζ(1−n)=−Bn, where ζ is the Riemann Zeta function.

n->if(n,-n*zeta(1-n),1)


Original solution, 52 bytes, using the generating function of Bernoulli numbers.

n->n!*polcoeff(-x/sum(i=1,n+1,(-x)^i/i!)+O(x^n*x),n)

• Can only upvote once. It is a shame it's not exact, though. Dec 1, 2015 at 10:05
• According to the documentation the zeta function is calculated using Bernoulli numbers, in fact. Dec 4, 2015 at 12:59
• @primo, yes, I regard all of the answers which use built-in zeta as cheating. Dec 4, 2015 at 17:59
• Even easier, bernfrac and bernreal are 8 bytes each and they're already functions, so no need for the n->. But +1 for a good solution. Aug 26, 2016 at 14:08

# Python 3, 112 bytes

Edit: I cleaned up this answer. If you want to see all the other ways I thought of to answer this question in Python 2 and 3, look in the revisions.

If I don't use the lookup table (and I use memoization instead), I manage to get the recursive definition to 112 bytes! WOO! Note that b(m) returns a Fraction. As usual, the byte count and a link for testing.

from fractions import*
def b(m):
s=k=0;p=1
while k<m:a=m-k;s+=Fraction(p*b(k))/-~a;p=p*a//-~k;k+=1
return 1-s


And a function that does use a lookup table, and returns the whole table of fractions from b(0) to b(m), inclusive.

from fractions import*
def b(m,r=[]):
s=k=0;p=1
while k<m:
if k>=len(r):r=b(k,r)
a=m-k;s+=Fraction(p*r[k])/-~a;p=p*a//-~k;k+=1
return r+[1-s]

• I think you can omit the trailing zero on float literals, e.g. 1. instead of 1.0. Nov 30, 2015 at 22:06
• @AlexA. Done. Removed .0 from s entirely, because it will quickly become a float later. Nov 30, 2015 at 22:13
• Can you use p=v=1;exec('[...];p+=1'*k) instead of your innermost loop? Dec 20, 2015 at 20:00

## CJam, 69 49 34 33 bytes

{_),:):R_:*_@f/@{_(;.-R.*}*0=\d/}


Online demo

Thanks to Cabbie407, whose answer made me aware of the Akiyama–Tanigawa algorithm.

### Dissection

{           e# Function: takes n on the stack
_),:)     e# Stack: n [1 2 3 ... n+1]
:R        e# Store that array in R
_:*       e# Stack: n [1 2 3 ... n+1] (n+1)!
_@f/      e# Stack: n (n+1)! [(n+1)!/1 (n+1)!/2 (n+1)!/3 ... (n+1)!/(n+1)]
e#   representing [1/1 1/2 ... 1/(n+1)] but avoiding floating point
@{        e# Repeat n times:
_(;.-   e#   Take pairwise differences
R.*     e#   Pointwise multiply by 1-based indices
}*        e#   Note that the tail of the array accumulates junk, but we don't care
0=\d/     e# Take the first element and divide by (n+1)!
}

• Multiplying through by n! to prevent precision loss is clever, if not slightly ridiculous. I wonder if the algorithm couldn't be refactored slightly to avoid this. Dec 5, 2015 at 3:37
• I don't think refactoring could avoid the need to scale for the simple reason that since we know that every other Bernoulli number is 0 there's obviously a lot of subtraction of similar values going on, so there's a lot of places where catastrophic loss of significance can occur. Dec 5, 2015 at 9:22

# PARI/GP, 45 bytes

n->if(n,2*n/(2^n-4^n)*real(polylog(1-n,I)),1)


Using the same formula as my Python answer, with An generated via polylog.

Test Script

Run gp, at the prompt paste the following:

n->if(n,2*n/(2^n-4^n)*real(polylog(1-n,I)),1)
for(i=0, 60, print(i, ": ", %(i)))

• Thank you for providing a test script - it made testing this a lot easier!
– user45941
Dec 1, 2015 at 11:40
• @Mego for you and me both ;) Dec 1, 2015 at 11:59

## Mathematica, 5248 42 bytes

1-Sum[#~Binomial~k#0@k/(#-k+1),{k,0,#-1}]&


Unnamed function that uses the literal definition.

• Is the Sign@# necessary? Dec 1, 2015 at 2:25
• I tested it on my computer. After removing the Sign@#, it still returns the correct answer for 0 . Dec 1, 2015 at 11:56

# Python 2, 132 130 bytes

import math,fractions
f=math.factorial
B=lambda m:~-m*m%2or 1+sum(B(k)*f(m)/f(k)/f(m-k)/fractions.Fraction(k+~m)for k in range(m))


This is just a golfed version of the reference implementation.

This is a little slow in practice, but can be sped up significantly with memoization:

import math,fractions
f=math.factorial

def memoize(f):
memo = {}
def helper(x):
if x not in memo:
memo[x] = f(x)
return memo[x]
return helper

@memoize
def B(m):
return~-m*m%2or 1+sum(B(k)*f(m)/f(k)/f(m-k)/fractions.Fraction(k+~m)for k in range(m))

for m in range(61):
print(B(m))


You can try this version online on Ideone.

# gawk4, 79 bytes

77 bytes code + 2 bytes for -M flag

PREC^=2{for(n=$0;m++<=n;)for($(j=m)=1/m;j>1;)$j=(-$j+$--j)*j;printf"%.6f",$1}


Had some trouble with the "6-decimal-digits-rule", because this prints the whole number and then 6 digits, but there is no list here to compare the results.

A flaw is that this prints a minus sign in front of the 0.000000 al lot of times, but I don't think that's wrong.

## Usage example

echo 58 | awk -M 'PREC^=2{for(n=$0;m++<=n;)for($(j=m)=1/m;j>1;)$j=(-$j+$--j)*j;printf"%.6f",$1}'


## Output from 0 to 60

0 -> 1.000000
1 -> 0.500000
2 -> 0.166667
3 -> -0.000000
4 -> -0.033333
5 -> 0.000000
6 -> 0.023810
7 -> 0.000000
8 -> -0.033333
9 -> 0.000000
10 -> 0.075758
11 -> -0.000000
12 -> -0.253114
13 -> -0.000000
14 -> 1.166667
15 -> -0.000000
16 -> -7.092157
17 -> -0.000000
18 -> 54.971178
19 -> -0.000000
20 -> -529.124242
21 -> -0.000000
22 -> 6192.123188
23 -> 0.000000
24 -> -86580.253114
25 -> 0.000000
26 -> 1425517.166667
27 -> 0.000000
28 -> -27298231.067816
29 -> 0.000000
30 -> 601580873.900642
31 -> 0.000000
32 -> -15116315767.092157
33 -> 0.000000
34 -> 429614643061.166667
35 -> 0.000000
36 -> -13711655205088.332772
37 -> 0.000000
38 -> 488332318973593.166667
39 -> -0.000000
40 -> -19296579341940068.148633
41 -> -0.000000
42 -> 841693047573682615.000554
43 -> -0.000000
44 -> -40338071854059455413.076812
45 -> -0.000000
46 -> 2115074863808199160560.145390
47 -> -0.000000
48 -> -120866265222965259346027.311937
49 -> -0.000000
50 -> 7500866746076964366855720.075758
51 -> -0.000000
52 -> -503877810148106891413789303.052201
53 -> -0.000000
54 -> 36528776484818123335110430842.971178
55 -> -0.000000
56 -> -2849876930245088222626914643291.067816
57 -> -0.000000
58 -> 238654274996836276446459819192192.149718
59 -> -0.000000
60 -> -21399949257225333665810744765191097.392674

• Would printf"%e" work? Dec 4, 2015 at 0:45
• No, it wouldn't, because the 0.00000s are only very small and not really zero. Dec 5, 2015 at 14:10

# R, 93 bytes

function(m){if(m==0){1}else{v=c();for(k in 0:(m-1))v=c(v,choose(m,k)*f(k)/(m-k+1));1-sum(v)}}


Not really original as a solution. If any comment, please feel free !

Ungolfed :

function(m)
if(m==0){1}
else
v=c()
for(k in 0:(m-1))
v=c(v,choose(m,k)*f(k)/(m-k+1))

1-sum(v)

• I know this is a bit late now but you can save 3 bytes by changing the if/else statement order and use m>0 as well as looping over 1:m-1 instead. Sep 28, 2016 at 7:31

# GolfScript, 63 bytes

~:i.!+.[3i&(2i>]*i(,{i\-,{1$1$(=2$2$)=+*2/}%\;}/~\2i?.(*\--1?**


Using the same formula as my Python answer.

Test Script

61,{[.
~:i.!+.[3i&(2i>]*i(,{i\-,{1$1$(=2$2$)=+*2/}%\;}/~\2i?.(*\--1?**
]p}/


The apphb link will time out on this. If you don't have GolfScript installed locally, I recommend using the anarchy golf interpreter (use form, select GolfScript, paste, submit).

# Perl, 101 bytes

#!perl -p
@a=($_%4-1,$_<2)x$_; @a=map$_*($a[$_-1]+$a[$_+1])/2,0..@a-3for 2..$_;$_=!$_||$_/(4**$_-2**$_)*$a[1]  Counting the shebang as three, input is taken from stdin. Using the same formula as my Python answer. Sample Usage $ echo 60 | perl bernoulli.pl
-2.13999492572253e+034


# Ruby, 66 61 bytes

This is a Ruby version of my Python answer.

b=->m{s,p=0r,1;m.times{|k|a=m-k;s+=p*b[k]/-~a;p=p*a/-~k};1-s}


Since this uses Rational in its answers, I'm fairly sure this works up to 60, but I had trouble running even b[24], so I implemented the lookup table again for 86 81 80 bytes.

t=->m{s,p,r=0r,1,m>0?t[m-1]:[];m.times{|k|a=m-k;s+=p*r[k]/-~a;p=p*a/-~k};r<<1-s}


# J, 10 bytes

(%1-^@-)t:


Computes the nth Bernoulli number by finding the nth coefficient of the exponential generating function of x/(1 - e-x).

## Usage

If the input is given integer or floats as an argument, it will output a float. If given an extended integer, marked with a suffix x, it will output either an extended integer or a rational, two extended integers separated by r.

   f =: (%1-^@-)t:
f 1
0.5
f 1x
1r2
(,.f"0) i. 10x
0     1
1   1r2
2   1r6
3     0
4 _1r30
5     0
6  1r42
7     0
8 _1r30
9     0


## Explanation

(%1-^@-)t: Input: n
computes the coefficients of its egf
(      )   A monad that operates on x
-      Negate x
^@       Computes its exponential, e^-x
1-         Subtract it from 1
%           Divide x by it, x/(1 - e^-x)


# Actually, 46 45 bytes (non-competing)

I've been meaning to do a Seriously/Actually answer for months and now I can. As this uses commands that Seriously didn't have in November of 2015, it's non-competing. Golfing suggestions welcome. Try it online!

Edit: In February of 2017, there was an update to Actually that changed which function literals are which. Normally, this would simply be non-competing for any challenge written before February, but since this answer is already non-competing, I edited this answer anyway. Enjoy.

This uses the explicit definition of the Bernoulli numbers on Wikipedia.

;╖ur⌠;╝ur⌠;;0~ⁿ(╛█*╜(uⁿ*⌡MΣ╛uk⌡M┬i@;π;)♀\*@k▼


Ungolfing

;╖     Duplicate and save m in register 0.
ur     Range [0..m]
⌠      Start first for loop
;╝     Duplicate and save k in register 1.
ur     Range [0..k]
⌠      Start second for loop (as string).
;;     Duplicate v twice.
0~ⁿ    Push -1, and pow() to get (-1)**v.
(╛█    Rotate a duplicate v to TOS, push k, and binom(k, v).
*      Multiply (-1)**v by binom(k, v).
╜(uⁿ   Push m, rotate last duplicate v to TOS, increment, and pow() to get (v+1)**m.
*      (-1)**v * binom(k, v) * (v+1)**m
⌡      End second for loop (string turned to function).
MΣ     Map over range [0..v] and sum
╛u     Push k and increment (the denominator)
(Note: second for loop does numerators only as denominator only depends on k)
k      Push fraction in list form [numerator, denominator]
⌡      End first for loop
M      Map over range [0..k]
┬i@    Transpose all of the fractions, flatten and swap.
Stack: [denominators] [numerators]
;π     Duplicate and take product of denominators.
;)     Duplicate product and move to bottom of stack.
Stack: product [denominators] [numerators] product
♀\     For all items in denominators, integer divide product by item.
Return a list of these scaled-up denominators.
*      Dot product of numerators and the scaled-up denominators as new numerator.
(In effect, getting the fractions to the same denominator and summing them)
@k     Swap new numerator and product (new denominator) and turn into a list (fraction).
▼      Divide fraction by gcd(numerator, denominator) (Simplify fraction).

• Uses commands Seriously didn't have in November 2015? Man, this uses an entirely new language that didn't exist in November 2015! I'm so proud...
– user45941
Aug 26, 2016 at 7:04

# Axiom, 134 147 bytes

b(n:NNI):FRAC INT==(v:=[1/1];k:=1;repeat(k>n=>break;r:=1-reduce(+,[binomial(k,j)*v.(j+1)/(k-j+1)for j in 0..k-1]);v:=append(v,[r]);k:=k+1);v.(n+1))


ungolf and test

(23) -> b
(23)
b n ==
1
v := [-]
1
k := 1
repeat
if n < k
then break
else
binomial(k,j)v(j + 1)
r := 1 - reduce(+,[[--------------------- for j in 0..(k - 1)]])
k - j + 1
v := append(v,[r])
k := k + 1
v(n + 1)
Type: FunctionCalled b
(50) -> [[i,b(i)]  for i in [0,1,2,3,4,5,6,7,8,9,10]]
(50)
1     1              1            1              1             5
[[0,1],[1,-],[2,-],[3,0],[4,- --],[5,0],[6,--],[7,0],[8,- --],[9,0],[10,--]]
2     6             30           42             30            66
Type: List List Fraction Integer

(51) -> b 1000
(51)
-
18243104738661887254572640256857788879338336867042906052197158157641126_
2572624911158657472577321069709615489924627495522908087488299539455188_
7918567582241551668492697244184914012242579830955617098629924652251740_
9791915637226361428342780548971002281045465308441161372350696920220116_
2441791760680262602019620260255790058416539271332852806000966628467639_
0683434226380702951226108116666172815817157023611889303668166839919156_
3797683877845690114843122753427426880591799883780255338278664578660218_
5045895962670442011443630321460259486764674312436994856054301765557425_
1371150213401051058408679874766352952749178734973676859834707623881634_
6251471489942512878190574323531299070406930309477389251738705417680653_
1183648189451892725726445949589759600705334767585389769924857630972963_
9976364832442643512622073858780110731539833099817555775136008111170797_
6250597322951308884900670113339167641953793994512377610306198429310933_
1214632141683542607746641232089854815064629129596536997380608256428801_
9784909897301658268809203555030692846151917069465607257641149187197651_
0905515966840312411845543650593021402849221691341852819791233589301994_
1012291773441794027493574651881059432274494354092231954894280742068472_
7146192942133436054611475404867886313250114399681532753236429290625909_
3411000391368336312138915621701535954814084208794241665492294270773347_
6055878415765927582014214726584822236443691314366097570085473354584000_
9985915190584047337934331297339403392719579093995842312746836871169674_
9786460913411872527166990047126222109345933847358924230951718379883743_
2563465487604170316077418754242710065269818190591271690695446633836120_
3745255515267088218383996330164203403732365333352120338272021319718003_
5994220458994876460018350270385634117807768745161622933834063145505621_
9106004731529642292049578901
/
342999030
Type: Fraction Integer

(52) -> b 60
1215233140483755572040304994079820246041491
(52)  - -------------------------------------------
56786730
Type: Fraction Integer


# APL(NARS), 83 chars, 166 bytes

r←B w;i
r←,1⋄i←0x⋄w+←1
→3×⍳w≤i+←1⋄r←r,1-+/{(1+i-⍵)÷⍨(⍵!i)×r[⍵+1]}¨0..i-1⋄→2
r←r[i]


Input as integer output as big rational

  B 0
1
B 1
1r2
B 2
1r6
B 3
0
B 4
¯1r30
B 10
5r66
B 100
¯94598037819122125295227433069493721872702841533066936133385696204311395415197247711r33330
B 1000
¯1824310473866188725457264025685778887933833686704290605219715815764112625726249111586574725773210697096154899246
27495522908087488299539455188791856758224155166849269724418491401224257983095561709862992465225174097919156
37226361428342780548971002281045465308441161372350696920220116244179176068026260201962026025579005841653927
13328528060009666284676390683434226380702951226108116666172815817157023611889303668166839919156379768387784
56901148431227534274268805917998837802553382786645786602185045895962670442011443630321460259486764674312436
99485605430176555742513711502134010510584086798747663529527491787349736768598347076238816346251471489942512
87819057432353129907040693030947738925173870541768065311836481894518927257264459495897596007053347675853897
69924857630972963997636483244264351262207385878011073153983309981755577513600811117079762505973229513088849
00670113339167641953793994512377610306198429310933121463214168354260774664123208985481506462912959653699738
06082564288019784909897301658268809203555030692846151917069465607257641149187197651090551596684031241184554
36505930214028492216913418528197912335893019941012291773441794027493574651881059432274494354092231954894280
74206847271461929421334360546114754048678863132501143996815327532364292906259093411000391368336312138915621
70153595481408420879424166549229427077334760558784157659275820142147265848222364436913143660975700854733545
84000998591519058404733793433129733940339271957909399584231274683687116967497864609134118725271669900471262
22109345933847358924230951718379883743256346548760417031607741875424271006526981819059127169069544663383612
03745255515267088218383996330164203403732365333352120338272021319718003599422045899487646001835027038563411
78077687451616229338340631455056219106004731529642292049578901r342999030


import Data.Ratio
p=product
b m=sum[p[k-v+1..k]*(v+1)^m%(p[-v..0-1]*(k+1))|k<-[0..m],v<-[0..k]]


This implements the explicit definition of Bernoulli numbers outlined on the Wikipedia page.

## Perl 6, 83 bytes

my &B={$^m??1-[+] (^$m).map: {combinations($m,$_)*B($_)/($m+1-$_)}!!1};say B slurp;  A faster, 114-byte solution: my @b=1;for 1..+slurp() {@b.push: 1-[+] (^$^m).map: {([*] $m+1-$_..$m)*@b[$_]/($m+1-$_)/([*] 1..\$_)}};say @b[*-1];

• Your code for a code golf challenge needs to be as short as possible, even if it takes a few lifetimes of the universe to terminate for certain inputs.
– user45941
Aug 26, 2016 at 7:47

## Javascript, 168 bytes

function h(b,a){return a?h(a,b%a):b}for(var c=[],a=[],e=0,b,d,f,g;e<=k;)for(c[b=d=e]=1,a[e]=++e;d;)f=c[--d]*a[b]-(c[b]*=g=a[d]),r=h(f*=b,g=a[b]*=g),c[d]=f/r,a[--b]=g/r;


Set the 'k' variable to the Bernoulli number you want, and the result is c[0] over a[0]. (numerator & denominator)

## Sample Usage

k = 2;
console.log(c[0] + "/" + a[0]);


Not as small as the others, but the only one I've written that comes close. See https://marquisdegeek.com/code_ada99 for my other (non-golf) attempts.

# Axiom, 57 bytes

g(n)==factorial(n)*coefficient(taylor(t*%e^t/(%e^t-1)),n)


code for test and results

(18) -> [[i, g(i)]  for i in 0..29]
(18)
1      1                1              1                1
[[0,1], [1,-], [2,-], [3,0], [4,- --], [5,0], [6,--], [7,0], [8,- --],
2      6               30             42               30
5                  691               7                 3617
[9,0], [10,--], [11,0], [12,- ----], [13,0], [14,-], [15,0], [16,- ----],
66                 2730               6                  510
43867                 174611               854513
[17,0], [18,-----], [19,0], [20,- ------], [21,0], [22,------], [23,0],
798                    330                  138
236364091               8553103                 23749461029
[24,- ---------], [25,0], [26,-------], [27,0], [28,- -----------], [29,0]]
2730                    6                        870
Type: List List Expression Integer

(19) -> g 60
1215233140483755572040304994079820246041491
(19)  - -------------------------------------------
56786730
Type: Expression Integer


one has to note that the function is not the one someone wrote above but t*%e^t/(%e^t-1)) with %e Euler costant

# Pyth, 22 bytes

L?b-1sm*.cbdcyd-btdUb1


Try it online!

Defines a function that is called as y<number>, e.g. yQ.

L                      # y=lambda b:
?b                  1 # ... if b else 1
-1                  # 1 -
s                 #     sum(
m            Ub  #         map(lambda d: ... , range(b))
*.cbd           #           combinations(b, d) *
cyd        #             y(d) /      (float division)
-btd    #                    b - (d - 1)
`