**KEYWORDS: Bernoulli numbers, asymptotic inclusion,
asymptotic approximation, OEIS A000367, OEIS A027641.**

There is a standard asymptotic formula for the Bernoulli numbers. For example at Mathworld or at the Digital Library of Mathematical Functions (NIST) (§24.11) one can find the formula

The present author found two inclusions for the Bernoulli numbers which appear to be new. They are reported on this page. The bounds given comprise a simple but amazingly efficient approximation to the Bernoulli numbers which I will call the 'cute approximation'.

Let B_{n} denote the Bernoulli numbers. If n
is even and n ≥ 38 then

For example the inclusion predicts

0.5318704469415522033..*10^1770 < |B(1000)| < 0.5318704469415522039..*10^1770.

And indeed |B(1000)| = 0.5318704469415522036..*10^1770.

Note that the factorial function is not referenced in these formulas.The lower bound of the above inclusion can also be used as a convenient approximation for the Bernoulli numbers. It then takes the form (assuming n even)

Equivalently this formula can be stated as

For example the standard approximation gives B(1000) ~ 0.5318*2*..*10^1770,
which amounts to **4** valid decimal numbers, whereas the last
approximation gives for B(1000) an approximation with almost **18.2**
valid decimal digits.

More generally let us look at three asymptotic formulas for the logarithm of the Bernoulli numbers.

- LogB1(n) = (1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n
- LogB2(n) = (1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*

(1-ln((120*n^2+9)/(120*n^2-1))) - LogB3(n) = (1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*

(1-1/(12*n^2)*(1-1/(30*n^2)*(1-2/(7*n^2))))

We see that formula 2 as well as formula 3 are refinements of formula 1. We are now focusing on formula 3 in its exponential form exp(LogB3(n)).

What makes this asymptotic approximation especially useful -- besides being
a *much* better approximation than those given by the Digital Library of
Mathematical Functions (see §24.11) -- is that the error of the approximation
can be easily estimated.

In fact a convenient way to reason about the validity of an approximation formula is to give a lower bound for the number of exact decimal digits, i. e. to indicate the number of decimal digits which are guaranteed by the formula at least. In the case of exp(LogB3(n)) we have a good and simple way to express this bound:

EddE(n) = floor(3*log(3*n)) (valid for n ≥ 50).

For example for Bernoulli(31622776) this bound says that 55 decimal digits of exp(LogB3(n)) are guaranteed to be correct (the true value is 55.7).

To sum up: to compute an approximation to the Bernoulli numbers B(n) with n ≥ 50 and n even compute exp(logB3(n)) and retain floor(3*log(3*n)) decimal digits of the result.

def BernoulliAsympt(n) : if n < 50 : print "Value error, n has to be ≥ 50" return if is_odd(n): return 0 R = RealField(300) # to be increased for large values n nn = n*n LogB = R((1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*(1-1/(12*nn)* (1-1/(30*nn)*(1-2/(7*nn)))) B = (-1)^(1+n//2)*exp(LogB) Edd = floor(3*ln(3*n)) print SciFormat(B, Edd - 1) # see SciFormat return B for n in (49..52): print n, BernoulliAsympt(n) 49 → Value error, n has to be ≥ 50 50 → 7.50086674607696e24 51 → 0 52 → -5.0387781014810e26 BernoulliAsympt(31622776) -7.66922063003368519879820408820523505875084131626143035e198196563

Note that the function displays only valid digits as we used this formatting function: SciFormat. With Sage the asymptotic formulas for the Bernoulli numbers based on these formulas and the exact decimal digits (edd) of these approximations can be computed as follows.

Edd1(n) = -log10(abs(1-exp(LogB1(n))/abs(bernoulli(n))) Edd2(n) = -log10(abs(1-exp(LogB2(n))/abs(bernoulli(n))) Edd3(n) = -log10(abs(1-exp(LogB3(n))/abs(bernoulli(n))) EddE(n) = 3*ln(3*n) ANF = 50; END = 1000; STEP = 20 plot2 = list_plot([[n, Edd2(n)] for n in range(ANF,END,STEP)], color='red') plot3 = list_plot([[n, Edd3(n)] for n in range(ANF,END,STEP)], color='blue') plotE = list_plot([[n, EddE(n)] for n in range(ANF,END,STEP)], color='magenta') show(plot2 + plot3 + plotE)

The plot below compares the number of exact decimal digits of the approximation (formula 3, blue curve) with the number of exact decimal digits guaranteed by this formula, 3*ln(3*n) (magenta curve). For comparison the red curve shows the exact decimal digits of formula 2. Formula 1 (given by DLMF/NIST) gives a poor approximation not worth to be shown.

Restating exp(LogB3(n)) gives the following asymptotic expansion of the Bernoulli number B(n) for n>0 even:

Note: The first inclusion was given on Jan. 18 2007 by the present author. Here the announcement in the newsgroup de.sci.mathematik. Thanks to Charles R Greathouse IV who drew my attention to an error in a previous version.

Asymptotic inclusions and approximations for the Euler numbers.

Asymptotic expansions of the factorial function.