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

An approximation of the Bernoulli numbers.

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

bernaasy

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'.

A sharp inclusion of the Bernoulli numbers.

Let Bn denote the Bernoulli numbers. If n is even and n ≥ 38 then

berinklusion

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.

A cute approximation of the Bernoulli numbers.

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)

bernapprox

 Equivalently this formula can be stated as

bernapprox

For example the standard approximation gives B(1000) ~ 0.53182..*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.

Asymptotic formulas

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

  1. LogB1(n) = (1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n
  2. 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)))
  3. 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.

previous Asymptotic inclusions and approximations for the Euler numbers.
next Asymptotic expansions of the factorial function.