[Paper Review] An efficient algorithm for the computation of Bernoulli numbers

G. J. Fee, Simon Plouffe|arXiv (Cornell University)|Feb 11, 2007

Advanced Mathematical Identities5 references4 citations

TL;DR

This paper presents a highly efficient algorithm for computing large Bernoulli numbers B(n) using an asymptotic approximation derived from the Euler product formula for the Riemann zeta function. By leveraging ζ(n) ≈ 1 for large even n, the method computes B(n) directly via B(n) ≈ 2n!πn / (2n), enabling the calculation of B(5,000,000) with a numerator of 27,332,507 decimal digits—surpassing prior records.

ABSTRACT

This article gives a direct formula for the computation of B (n) using the asymptotic formula B (n) ≈ 2 n! π n 2 n where n is even and n ≫ 1. This is simply based on the fact that ζ (n) is very near 1 when n is large and since B (n) = 2 ζ(n)n! πn2n exactly. The formula chosen for the Zeta function is the one with prime numbers from the well-known Euler product for ζ (n). This algorithm is far better than the recurrence formula for the Bernoulli numbers even if each B(n) is computed individually. The author could compute B (750,000) in a few hours. The current record of computation is now (as of Feb. 2007) B (5,000,000) a number of (the numerator) of 27332507 decimal digits is also based on that idea. 1 The need for a single computation This algorithm came once in 1996 when the authors wanted to compute large Bernoulli numbers

Motivation & Objective

• To develop a faster alternative to the traditional recurrence-based computation of Bernoulli numbers for large indices.
• To enable the practical computation of extremely large Bernoulli numbers, such as B(5,000,000), which were previously infeasible with standard methods.
• To exploit the asymptotic behavior of the Riemann zeta function ζ(n) ≈ 1 for large even n to simplify and accelerate the calculation of B(n).
• To provide a direct, non-recursive algorithm that avoids the computational overhead of recurrence relations.
• To set a new record in computing Bernoulli numbers by achieving B(5,000,000) in a few hours using the proposed method.

Proposed method

• The algorithm uses the exact formula B(n) = 2ζ(n)n!πn / (2n) for even n, with ζ(n) approximated as 1 when n ≫ 1.
• It applies the Euler product representation of ζ(n) over primes to refine the asymptotic approximation for improved accuracy.
• The method computes B(n) directly without relying on recurrence relations, reducing computational complexity.
• The algorithm leverages the fact that ζ(n) approaches 1 rapidly for large even n, simplifying the evaluation of B(n).
• High-precision arithmetic is used to handle the massive size of the numerator of B(n), which reaches 27,332,507 decimal digits for B(5,000,000).
• The computation is optimized for individual B(n) values, making it suitable for targeted large-scale calculations.

Experimental results

Research questions

• RQ1 Can a direct asymptotic formula based on the zeta function's behavior for large n significantly outperform recurrence-based methods in computing Bernoulli numbers?
• RQ2 To what extent can the Euler product formula for ζ(n) improve the accuracy of the asymptotic approximation for large even n?
• RQ3 Is it feasible to compute Bernoulli numbers with millions of digits in their numerator using a non-recursive, direct algorithm?
• RQ4 What is the computational performance ceiling of such an algorithm in terms of the largest achievable B(n)?
• RQ5 Can this method set a new record for the largest computed Bernoulli number, surpassing previous benchmarks?

Key findings

• The algorithm computes B(750,000) in just a few hours, demonstrating high efficiency for large-scale computations.
• The method achieves a new record by computing B(5,000,000), whose numerator contains 27,332,507 decimal digits.
• The asymptotic approximation B(n) ≈ 2n!πn / (2n) is sufficiently accurate for large even n due to ζ(n) ≈ 1.
• The use of the Euler product for ζ(n) enhances the precision of the asymptotic formula without introducing significant computational overhead.
• The algorithm outperforms traditional recurrence-based methods even when computing individual B(n) values.
• The approach enables feasible computation of Bernoulli numbers far beyond the reach of standard recurrence techniques.