Bernoulli numbers form a sequence of rational numbers that arise frequently in number theory and analysis, defined by the exponential generating function \frac{x}{e^x - 1} = \sum_{n=0}^\infty \frac{B_n x^n}{n!}, where B_n denotes the nth Bernoulli number.[1] Named after the Swiss mathematician Jakob Bernoulli (1654–1705), who introduced them in his 1713 posthumous work Ars Conjectandi to express formulas for the sums of powers of integers, such as \sum_{k=1}^m k^p = \frac{1}{p+1} \sum_{j=0}^p \binom{p+1}{j} B_j m^{p+1-j}.[2] These numbers were independently discovered around the same time by the Japanese mathematician Seki Takakazu (1642–1708), though Bernoulli's publication brought them prominence in Western mathematics.[3]The sequence begins with B_0 = 1, B_1 = -\frac{1}{2}, B_2 = \frac{1}{6}, B_3 = 0, B_4 = -\frac{1}{30}, and so on, exhibiting the property that B_n = 0 for all odd n > 1, while the even-indexed terms alternate in sign and grow rapidly in magnitude.[1] Leonhard Euler formalized the modern generating function in 1755 and extended the concept to Bernoulli polynomials B_n(x), which satisfy B_n(0) = B_n and appear in the Euler-Maclaurin summation formula, a key tool for approximating sums by integrals: \sum_{k=a}^b f(k) \approx \int_a^b f(x) \, dx + \frac{f(a) + f(b)}{2} + \sum_{k=1}^m \frac{B_{2k}}{(2k)!} (f^{(2k-1)}(b) - f^{(2k-1)}(a)) + R.[2] Bernoulli numbers also connect deeply to the Riemann zeta function, providing explicit values for even positive integers via \zeta(2n) = (-1)^{n+1} \frac{B_{2n} (2\pi)^{2n}}{2 (2n)!}, which Euler used to evaluate \zeta(2) = \frac{\pi^2}{6}.[3]Beyond these classical applications in series expansions of trigonometric functions like \cot x and \tan x, Bernoulli numbers play roles in modern contexts such as the von Staudt–Clausen theorem, which describes their denominators in terms of primes, and Kummer's congruence relating them to Fermat's Last Theorem for regular primes.[1] Their denominators are the product of the distinct primes p for which p-1 divides 2n, ensuring they are rational but not always integer-valued.[2] Overall, the Bernoulli numbers exemplify the intricate links between combinatorial identities, analytic continuation, and arithmetic properties in mathematics.
Fundamentals
Notation
The Bernoulli numbers are denoted by B_n, where n is a non-negative integer indexing the sequence of rational numbers, with the zeroth term defined as B_0 = 1. In the standard modern convention, the first Bernoulli number is B_1 = -\frac{1}{2}, a choice that aligns with applications in analysis and number theory, such as the Euler-Maclaurin formula.[4] This notation is prescribed by the National Institute of Standards and Technology (NIST) and adopted in most contemporary mathematical literature.[4]Historically, variations in sign conventions have existed, particularly for the first index. In Jacob Bernoulli's original presentation and some subsequent older texts, such as those by Édouard Lucas and Ernesto Cesàro, the convention set B_1 = +\frac{1}{2}, with the signs for higher odd indices adjusted accordingly.[5] To reconcile these differences and maintain consistency for even indices—which remain unchanged across conventions—the alternative sequence with B_1^* = +\frac{1}{2} is often denoted by B_n^*, while the modern B_n retains the negative sign for B_1.[1] This starred notation distinguishes the oldervariant without altering the values for n \neq 1.[1]Regarding indices, Bernoulli numbers with odd indices greater than 1 are zero in both conventions, so B_{2k+1} = 0 for all integers k \geq 1, emphasizing the focus on even-indexed terms B_{2k} in many contexts.[4] The numbers are named after the Swiss mathematician Jakob Bernoulli (1654–1705), who introduced them in his posthumously published work Ars Conjectandi (1713).[5]
Definitions
The Bernoulli numbers B_n form a sequence of rational numbers defined recursively by setting B_0 = 1 and, for each integer n > 1,\sum_{k=0}^{n-1} \binom{n}{k} B_k = 0.This relation allows the computation of each subsequent B_n from the previous ones, yielding B_n = -\frac{1}{n} \sum_{k=0}^{n-1} \binom{n}{k} B_k for n > 1, with the value of B_1 = -\frac{1}{2} fixed to satisfy the defining properties.[6] Once B_1 is specified, the recursion uniquely determines all higher B_n.[3]This recursive definition arises as the coefficients in the Taylor series expansion of the function \frac{x}{e^x - 1} around x = 0, given by\frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!}.The function \frac{x}{e^x - 1} appears naturally in contexts such as the Euler-Maclaurin formula and analytic number theory, motivating the choice of these coefficients.[6]The Bernoulli numbers considered here are the classical ones, referred to as Bernoulli numbers of the first kind, with B_1 = -\frac{1}{2}. In contrast, Bernoulli numbers of the second kind, denoted b_n, are defined by the generating function \frac{x}{\ln(1+x)} = \sum_{n=0}^\infty b_n \frac{x^n}{n!}, yielding b_1 = \frac{1}{2} and differing values for higher indices.[7]
History
Early history
The early history of Bernoulli numbers traces back to 17th-century efforts to find closed-form expressions for the sums of successive powers of positive integers, a problem that had intrigued mathematicians since antiquity. Blaise Pascal made significant strides in this direction around 1654 in his Traité du triangle arithmétique, where he developed methods using binomial coefficients to express such sums as polynomials, implicitly relying on patterns akin to finite differences for computing higher powers iteratively. These approaches provided practical tools but lacked a unified general framework, setting the stage for further algebraic investigations.Johann Faulhaber advanced these ideas in his 1631 publication Academia Algebrae, which included a section titled "Summae Potestatum" featuring extensive tables of formulas for sums of powers up to the 17th degree.[8] Faulhaber's work presented specific polynomial expressions, often in terms of triangular numbers N = n(n+1)/2, but offered no explicit general formula or proof, relying instead on recursive summation techniques and geometric analogies.[8] Modern historians, notably Donald Knuth in a 1993 analysis, have reconstructed Faulhaber's lost original manuscript and methods, revealing how his coefficients foreshadowed the structure of Bernoulli numbers through matrix-based derivations and reflective polynomials.[8]Independently, the Japanese mathematician Seki Takakazu (1642–1708) discovered similar coefficients around the same period, though his work remained less known in the West until later.[3]The numbers now known as Bernoulli numbers emerged implicitly in the work of Jakob Bernoulli, who, building on Faulhaber's tables and Pascal's binomial methods, sought a general solution for power sums in the late 17th century. In his posthumously published Ars Conjectandi (1713), Bernoulli included a dedicated chapter on "Summae Potestatum," listing formulas up to the 10th power and identifying recurring coefficients that unified the expressions as polynomials of degree one higher than the power.[9] These coefficients, appearing without a named sequence, represented the first systematic tabulation of what Euler would later formalize and name in the 18th century.
Euler's contributions
Leonhard Euler played a pivotal role in formalizing the theory of Bernoulli numbers during the 18th century, building on earlier groundwork by the Bernoulli family through his innovative approaches to series summation and function expansions. In his 1738 paper "Methodus generalis summandi progressiones," published in the Commentarii Academiae Scientiarum Petropolitanae, Euler introduced the exponential generating function for the Bernoulli polynomials, which yields the Bernoulli numbers as the coefficients when evaluated at zero. The function is expressed as\frac{x e^{y x}}{e^x - 1} = \sum_{n=0}^\infty B_n(y) \frac{x^n}{n!},where B_n(y) are the Bernoulli polynomials and B_n = B_n(0) are the Bernoulli numbers, satisfying\frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!}.This formulation provided a systematic way to generate the sequence and extended the numbers to polynomial forms useful in integral calculus and series manipulations.Euler computed the initial terms of the Bernoulli numbers using this generating function, obtaining B_0 = 1, B_1 = -\frac{1}{2}, B_2 = \frac{1}{6}, B_3 = 0, B_4 = -\frac{1}{30}, and noting the pattern that B_n = 0 for odd n > 1. These values allowed him to verify and expand upon known results for sums of powers. He recognized the central role of Bernoulli numbers in deriving explicit formulas for the sums of p-th powers of the first m natural numbers, expressing \sum_{k=1}^m k^p as a polynomial of degree p+1 in m, with coefficients involving the Bernoulli numbers up to B_p. For example, for p=2, the sum is \frac{m(m+1)(2m+1)}{6}, aligning with B_2 = \frac{1}{6}. This general method streamlined computations that had previously relied on ad hoc techniques.Euler further applied Bernoulli numbers to the Taylor series expansions of trigonometric functions, particularly the cotangent. In this work and related publications, he derived the expansion\frac{\pi \cot(\pi z)}{z} = \frac{1}{z^2} + \sum_{n=1}^\infty (-1)^n \frac{(2\pi)^{2n} B_{2n}}{(2n)!} z^{2n-2},which connects the numbers to infinite products for the sine function and aids in evaluating series like those for the Riemann zeta function at even integers. Similar uses appeared in expansions for the cosecant and other functions, highlighting the numbers' utility in analytic number theory.Euler's advancements emerged amid active correspondence with contemporaries, including Johann I Bernoulli and Christian Goldbach, where he discussed power sums and series as early as the 1730s; these exchanges, documented in letters from 1731 onward, informed his publications and spurred further developments in the field.[10]
Modern developments
In the 19th century, significant advancements in the arithmetical properties of Bernoulli numbers were made by Ernst Kummer, who linked them to the theory of cyclotomic fields and Fermat's Last Theorem. Kummer introduced the concepts of regular and irregular primes in 1850, defining an odd prime p as irregular if it divides the numerator of any Bernoulli number B_k for even k with $2 \leq k \leq p-3. This classification stemmed from Kummer's study of the divisibility of the class number of the p-th cyclotomic field by p, where irregular primes obstruct ideal factorization in these fields, impacting proofs of Fermat's Last Theorem for certain exponents. Subsequent work by mathematicians like Jensen in 1915 disproved Kummer's conjecture of finitely many irregular primes, showing that there are infinitely many such primes.[11]In the 20th century, Bernoulli numbers gained deeper ties to analytic number theory through their explicit relation to the Riemann zeta function, \zeta(2n) = (-1)^{n+1} \frac{B_{2n} (2\pi)^{2n}}{2 (2n)!} for positive integers n, which connects them to the distribution of prime numbers.[12] This relation underpins the explicit formulas in prime number theory, such as von Mangoldt's formula, where sums over the non-trivial zeros of the zeta function describe oscillations in the prime counting function; the locations of these zeros are conjectured by the Riemann hypothesis to lie on the critical line \Re(s) = 1/2, offering a tantalizing analytic bridge from Bernoulli numbers to unresolved questions about prime spacing.[13]Post-2000 developments have leveraged arithmetic geometry for high-precision computations of Bernoulli numbers, particularly using properties of modular forms to evaluate them modulo primes efficiently. In his 2007 book, William Stein outlined algorithms that compute generalized Bernoulli numbers via Hecke operators on modular forms, enabling practical calculations for large indices by exploiting the arithmetic of modular Jacobians.[14] Building on this, David Harvey's 2008 multimodular algorithm achieved a record computation of B_{10^8} using parallel processing and modular reductions, reducing time complexity to near-optimal levels and facilitating applications in number-theoretic software like SageMath.[15] These methods highlight the interplay between analytic formulas and geometric structures, surpassing earlier recursive approaches in scalability.Emerging applications in modern physics have revealed Bernoulli numbers in the partition functions of topological string theory, where they encode higher-genus contributions to amplitudes on Calabi-Yau manifolds. In 2004, Dunne and Schubert derived convolution identities for products of Bernoulli numbers from quantum field theory computations, linking them to perturbative expansions in topological strings.[16] Similarly, in two-dimensional topological string models, the partition function at higher genera is proportional to Bernoulli numbers, providing exact non-perturbative insights into enumerative invariants and BPS states. These connections underscore the numbers' role in unifying algebraic geometry with physical dualities.
Generating Functions and Explicit Formulas
Exponential generating function
The exponential generating function for the Bernoulli numbers B_n is defined by the equation\frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!},valid in a neighborhood of x = 0. This analytic expression provides a closed-form representation that bridges the recursive definition of the Bernoulli numbers to their role in series expansions and summation formulas.[1]To derive this generating function using a geometric series expansion, begin with the identity\frac{1}{e^x - 1} = \frac{e^{-x}}{1 - e^{-x}} = e^{-x} \sum_{k=0}^\infty e^{-k x} = \sum_{k=1}^\infty e^{-k x},where the infinite sum follows from the geometric series formula \sum_{k=0}^\infty r^k = 1/(1 - r) with r = e^{-x} and |r| < 1 for \Re(x) > 0. Multiplying both sides by x yields\frac{x}{e^x - 1} = x \sum_{k=1}^\infty e^{-k x}.Expanding each exponential term e^{-k x} = \sum_{n=0}^\infty \frac{(-k x)^n}{n!} and collecting coefficients leads to the power series form.[17]A key property of this generating function arises from its connection to the cotangent function through partial fraction decompositions. Specifically, the function \pi \cot(\pi z) admits the partial fraction expansion \pi \cot(\pi z) = 1/z + \sum_{m=1}^\infty \left[ 1/(z - m) + 1/(z + m) \right], and differentiating or integrating this leads to relations involving \frac{x}{e^x - 1}. In particular, \frac{x}{e^x - 1} = \frac{x}{2} \coth\left(\frac{x}{2}\right) - \frac{x}{2}, linking the generating function to hyperbolic functions and, by analytic continuation, to the trigonometric cotangent series \sum_{n=0}^\infty (-1)^n \frac{2^{2n} B_{2n} x^{2n}}{(2n)!} = x \cot x for even indices.[1]The power series converges for |x| < 2\pi, determined by the distance to the nearest singularity of \frac{x}{e^x - 1} at x = 2\pi i, where e^x = 1. Beyond this radius, the function admits analytic continuation as a meromorphic function in the complex plane, with simple poles at x = 2\pi i k for nonzero integers k, allowing extension of Bernoulli number properties to broader analytic contexts.[1]
Explicit summation formula
The explicit summation formula for Bernoulli numbers provides a closed-form expression that allows direct computation without recursion, derived from the relation to sums of powers and the generating function \frac{x}{e^x - 1}. One standard form, known as Worpitzky's formula, expresses the nth Bernoulli number B_n as a double sum involving binomial coefficients and powers:B_n = \sum_{k=0}^n \frac{1}{k+1} \sum_{j=0}^k (-1)^j \binom{k}{j} j^n.This formula arises from inverting Faulhaber's expression for the sum of mth powers of the first n integers, \sum_{j=1}^n j^m = \frac{1}{m+1} \sum_{k=0}^m \binom{m+1}{k} B_k n^{m+1-k}, by considering the case n = k+1 and substituting the expression for powers in terms of finite differences, leading to the inner alternating sum.[18]An alternative explicit representation involves Stirling numbers of the second kind S(n,k), which count the number of ways to partition n objects into k non-empty subsets. The formula isB_n = \sum_{k=1}^n (-1)^k \frac{k! \, S(n,k)}{k+1}.This follows briefly from the identity expressing powers as falling factorials via Stirling numbers, j^n = \sum_{k=0}^n S(n,k) (j)_k, where (j)_k = j(j-1)\cdots(j-k+1), and summing over j=1 to some limit; the inversion using the Bernoulli relation for sums of falling factorials yields the expression after collecting coefficients.[7]Euler's early work on Bernoulli numbers included attempts to compute them explicitly by solving finite systems from the power sum formula, calculating the first several terms up to B_{10} in the 1730s to evaluate series for the zeta function at even integers, though he relied on recursive methods rather than fully closed forms like those above.[3]These summation formulas, while explicit, face practical limitations for large n due to the alternating signs in the inner sums, which produce terms of rapidly growing magnitude (exponentially with n) that largely cancel to yield the small rational values of B_n, necessitating high-precision arithmetic to avoid numerical instability.[18]
Relations to Analytic Functions
Connection with Riemann zeta function
The connection between Bernoulli numbers and the Riemann zeta function arises prominently at positive even integers, providing explicit evaluations that link these fundamental objects in number theory. The Riemann zeta function, defined for \operatorname{Re}(s) > 1 by \zeta(s) = \sum_{k=1}^\infty k^{-s}, takes a closed form at even positive integers s = 2n involving the $2n-th Bernoulli number:\zeta(2n) = (-1)^{n+1} \frac{B_{2n} (2\pi)^{2n}}{2 (2n)!}, \quad n = 1, 2, 3, \dotsThis formula, first derived by Leonhard Euler, expresses \zeta(2n) in terms of \pi and Bernoulli numbers, enabling precise computations such as \zeta(2) = \pi^2/6 and \zeta(4) = \pi^4/90.[1][19]Euler established this relation in the 18th century, beginning with his 1734 solution to the Basel problem, which evaluated \zeta(2) using the infinite product representation of the sine function, \sin(\pi x)/(\pi x) = \prod_{k=1}^\infty (1 - x^2/k^2). By taking the logarithmic derivative, Euler obtained a partial fraction expansion for \pi \cot(\pi x) whose coefficients involved sums like \zeta(2n), allowing him to equate power series expansions and isolate the zeta values. He extended this approach in 1740 to derive the general formula for \zeta(2n) in terms of Bernoulli numbers, which appear naturally in the Taylor series for \cot(\pi x) or related trigonometric functions.[20][21][22]The formula has key implications for the properties of even-indexed Bernoulli numbers. It reveals that the signs of B_{2n} alternate as (-1)^{n-1} for n \geq 1 (e.g., B_2 = 1/6 > 0, B_4 = -1/30 < 0, B_6 = 1/42 > 0), ensuring \zeta(2n) > 0 despite the alternating factor (-1)^{n+1}. Furthermore, since \zeta(2n) \to 1 as n \to \infty, the magnitude of the even Bernoulli numbers grows asymptotically as |B_{2n}| \sim 2 (2n)! / (2\pi)^{2n}, reflecting their rapid super-exponential increase, which is crucial for convergence in series expansions involving these numbers.[1][23]This connection can also be approached through the ordinary generating function for Bernoulli numbers, whose relation to trigonometric functions yields the zeta evaluations upon suitable differentiation and evaluation.[1]
Integral representations
The Bernoulli numbers possess several integral representations that extend their definition beyond algebraic or series forms, facilitating both theoretical analysis and numerical computation.A fundamental contourintegral representation is given byB_n = \frac{n!}{2\pi i} \oint_{|z|=r} \frac{z^{-n}}{e^z - 1} \, dz,where the positively oriented circular contour encloses the origin with radius $0 < r < 2\pi, avoiding the poles of the integrand at multiples of $2\pi i. This expression derives from applying the residue theorem to extract the coefficient of z^n/n! in the generating function \frac{z}{e^z - 1}.[1]For even indices, a real integral representation, often attributed to Riesz, isB_{2n} = (-1)^{n+1} 4n \int_0^\infty \frac{t^{2n-1}}{e^{2\pi t} - 1} \, dt, \quad n \geq 1.This absolutely convergent integral provides an explicit path for evaluating even Bernoulli numbers and connects to broader analytic continuations in number theory.[24]The Bernoulli numbers are linked to the Gamma function through the analytic continuation of the Riemann zeta function to negative arguments. For integers n \geq 2,B_n = -n \, \zeta(1 - n),where \zeta(s) denotes the Riemann zeta function. Substituting into the functional equation of the zeta function,\zeta(s) = 2^s \pi^{s-1} \sin\left( \frac{\pi s}{2} \right) \Gamma(1 - s) \, \zeta(1 - s),evaluated at s = 1 - n, yieldsB_n = -n \cdot 2^{1-n} \pi^{-n} \sin\left( \frac{\pi (1 - n)}{2} \right) (n-1)! \, \zeta(n).Here, \Gamma(n) = (n-1)! explicitly incorporates the Gamma function, offering a bridge to positive arguments of \zeta(s) and enabling computations via known values such as \zeta(2k) for even positive integers.[19]Fourier series representations emerge from the periodic extensions of functions related to Bernoulli polynomials, such as the linear function x on [-\pi, \pi]. The periodic extension with period $2\pi has the Fourier seriesx = 2 \sum_{k=1}^\infty \frac{(-1)^{k+1}}{k} \sin(kx), \quad -\pi < x < \pi.Extensions to higher powers x^m on the same interval yield Fourier series whose coefficients involve Bernoulli numbers, as seen in the expansions of periodic Bernoulli functions \tilde{B}_m(x), providing integral-free but series-based expressions tied to trigonometric sums.[25]These contour integrals are utilized in computational algorithms for Bernoulli numbers by deforming the integration path to enclose additional poles at z = 2\pi i k (k \in \mathbb{Z} \setminus \{0\}), allowing evaluation via residue calculus at those points. This approach yields recurrence relations and asymptotic estimates for large n, improving efficiency over direct series summation.[26]
Computation Methods
Recursive computation
The Bernoulli numbers B_m for m \geq 1 can be computed recursively using the defining relation derived from their generating function \frac{x}{e^x - 1} = \sum_{m=0}^\infty B_m \frac{x^m}{m!}, which leads to the equation \sum_{k=0}^m \binom{m+1}{k} B_k = 0, with B_0 = 1. Solving for the highest-index term gives the explicit recursion B_m = -\frac{1}{m+1} \sum_{k=0}^{m-1} \binom{m+1}{k} B_k.[2]To compute the first few Bernoulli numbers step by step, begin with B_0 = 1. For m=1, the sum is \binom{2}{0} B_0 = 1, so B_1 = -\frac{1}{2} (1) = -\frac{1}{2}. For m=2, the sum up to k=1 is \binom{3}{0} B_0 + \binom{3}{1} B_1 = 1 + 3 \left(-\frac{1}{2}\right) = -\frac{1}{2}, so B_2 = -\frac{1}{3} \left(-\frac{1}{2}\right) = \frac{1}{6}. For m=3, the sum up to k=2 is \binom{4}{0} B_0 + \binom{4}{1} B_1 + \binom{4}{2} B_2 = 1 + 4 \left(-\frac{1}{2}\right) + 6 \left(\frac{1}{6}\right) = 0, so B_3 = -\frac{1}{4} (0) = 0. For m=4, the sum up to k=3 is \binom{5}{0} B_0 + \binom{5}{1} B_1 + \binom{5}{2} B_2 + \binom{5}{3} B_3 = 1 + 5 \left(-\frac{1}{2}\right) + 10 \left(\frac{1}{6}\right) + 10 (0) = \frac{1}{6}, so B_4 = -\frac{1}{5} \left(\frac{1}{6}\right) = -\frac{1}{30}.[2]This recursive algorithm requires computing binomial coefficients and performing arithmetic on the previous Bernoulli numbers at each step, resulting in a time complexity of O(n^2) to compute the first n numbers, as each of the n steps involves a sum of up to O(n) terms.For practical implementation, especially to compute up to moderate orders such as n \leq 10, it is essential to represent the Bernoulli numbers as exact rational fractions to avoid accumulation of floating-point errors and numerical instability, which can arise from repeated divisions and summations in floating-point arithmetic.[2]
Efficient numerical algorithms
Efficient numerical algorithms for computing Bernoulli numbers have advanced significantly since the early 2000s, enabling calculations for indices up to millions or more with practical efficiency. These methods address the exponential growth in the size of Bernoulli numbers, which have approximately n \log n bits for index n, by leveraging modular arithmetic, fast transforms, and specialized identities rather than direct recursion.[15]A seminal approach is the multimodular algorithm developed by David Harvey in 2010, which computes individual Bernoulli numbers B_n by evaluating them modulo numerous small primes and reconstructing the result via the Chinese Remainder Theorem. This method exploits the connection between Bernoulli numbers and values of the Riemann zeta function at even integers, B_{2k} = (-1)^{k+1} \frac{2 (2k)! \zeta(2k)}{(2\pi)^{2k}}, but performs the core computations in finite fields to avoid high-precision floating-point arithmetic. The fast Fourier transform (FFT) accelerates the evaluation of power sums modulo primes, achieving a time complexity of O(n^2 (\log n)^{2+o(1)}) bit operations. Using a parallel implementation, this algorithm set a record by computing B_{10^8}, demonstrating feasibility for n up to at least $10^6 on standard hardware.The Von Staudt–Clausen theorem plays a supporting role in these computations by providing explicit congruences that facilitate modular reductions, particularly for determining denominators and adjusting values modulo primes where p-1 divides the index. This allows efficient handling of the rational structure of Bernoulli numbers during reconstruction, reducing the number of moduli needed for accuracy.For approximations modulo a prime p, p-adic methods offer further optimization, especially when higher powers p^s are required. Harvey's p-adic multimodular extension, introduced in 2012, generalizes Voronoi's congruence to compute B_{2n} \mod p^s using a recursive formula involving binomial coefficients and power sums over residues modulo p. By choosing parameters like s \approx n^{2/3}, it achieves subquadratic complexity O(n^{4/3 + o(1)}), making it suitable for large-scale modular computations in number-theoretic applications.[27]In the 2020s, advancements have focused on parallelization and formal verification to scale these algorithms further. For instance, a 2025 verified implementation of Harvey's multimodular method in Isabelle/HOL emphasizes the inherent parallelism of the modular power sum evaluations, which can be distributed across arithmetic progressions of residues, enabling efficient computation on multicore systems with minimal memory overhead. This work highlights ongoing efforts to bridge theoretical efficiency with practical, verified software, addressing gaps in earlier implementations.[28]
Applications in Analysis
Sums of integer powers
One of the most prominent applications of Bernoulli numbers arises in evaluating the sums of the p-th powers of the first m positive integers, expressed through Faulhaber's formula. This formula provides a closed-form polynomial expression in terms of m, with coefficients involving Bernoulli numbers. Specifically,\sum_{k=1}^m k^p = \frac{1}{p+1} \sum_{j=0}^p (-1)^j \binom{p+1}{j} B_j m^{p+1-j},where the Bernoulli numbers B_j (with the convention B_1 = -\frac{1}{2}) ensure the right-hand side is a polynomial of degree p+1.[29]The derivation of Faulhaber's formula can be obtained using the exponential generating function for Bernoulli numbers, \frac{x}{e^x - 1} = \sum_{j=0}^\infty B_j \frac{x^j}{j!}, which encodes the necessary summation properties through operator methods or finite differences.[30] Alternatively, it follows from the calculus of finite differences, where the forward difference operator \Delta f(x) = f(x+1) - f(x) applied to the monomial f(x) = x^p leads to the antiderivative-like sum via Bernoulli polynomials, ultimately yielding the explicit form with Bernoulli numbers.[3]For small values of p, the formula simplifies to well-known polynomial expressions. For p=1,\sum_{k=1}^m k = \frac{m(m+1)}{2},which matches the general formula using B_0 = 1 and B_1 = -\frac{1}{2}. For p=2,\sum_{k=1}^m k^2 = \frac{m(m+1)(2m+1)}{6},incorporating B_2 = \frac{1}{6}. For p=3,\sum_{k=1}^m k^3 = \left( \frac{m(m+1)}{2} \right)^2,where the vanishing B_3 = 0 eliminates higher terms beyond the quartic. These examples illustrate how Bernoulli numbers systematically determine the coefficients, transforming the sum into an explicit polynomial.[3]Historically, the quest for such formulas dates to the early 17th century, when Johann Faulhaber computed explicit expressions for sums up to the 23rd power in his 1631 work Academia Algebræ, laying foundational techniques for power summation without reference to Bernoulli numbers. Jacob Bernoulli later unified these results in 1713 using his newly defined numbers, providing the general form and marking a key advancement in early calculus by enabling precise evaluations essential for integral approximations and series developments.[31]
Asymptotic expansions
The Euler-Maclaurin formula establishes a connection between definite integrals and sums, providing an asymptotic expansion where Bernoulli numbers serve as coefficients in the higher-order correction terms.[32] This formula, developed independently by Leonhard Euler and Colin Maclaurin in the 18th century, approximates the integral of a smooth function f(x) over an interval [a, b] by a midpoint rule enhanced with terms involving derivatives of f at the endpoints.[33] The expansion is particularly useful for deriving asymptotic behaviors in analysis and numerical computations.The symmetric form of the Euler-Maclaurin formula is given by\int_a^b f(x) \, dx \approx (b - a) f\left( \frac{a + b}{2} \right) + \sum_{k=1}^m \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(b) - f^{(2k-1)}(a) \right) + R_m,where B_{2k} are the even-indexed Bernoulli numbers, f^{(2k-1)} denotes the (2k-1)-th derivative, and R_m is the remainder term.[32] Only even powers of Bernoulli numbers appear because the odd-indexed ones (beyond B_1) vanish, reflecting the formula's structure derived from Bernoulli polynomials.[34] This approximation improves accuracy by incorporating successively higher derivatives, with the series converging asymptotically for analytic functions.A prominent application arises in deriving Stirling's approximation for the factorial n!, obtained by applying the Euler-Maclaurin formula to the sum \sum_{k=1}^n \log k \approx \int_1^n \log x \, dx plus corrections.[35] The resulting asymptotic series for \log n! is \log n! \sim n \log n - n + \frac{1}{2} \log (2\pi n) + \sum_{k=1}^\infty \frac{B_{2k}}{2k(2k-1) n^{2k-1}}, where the Bernoulli numbers dictate the subleading terms, providing increasingly precise estimates for large n.[36]The error term R_m in the expansion is typically bounded by an integral involving the (2m+1)-th derivative of f and the next Bernoulli polynomial, but its magnitude grows rapidly due to the factorial growth of |B_{2m}| for large m, limiting practical truncation to low orders unless f is highly smooth.[13] Analysis of this remainder highlights the asymptotic nature of the series, where higher Bernoulli numbers amplify errors beyond optimal truncation points.[37]In physics, the Euler-Maclaurin formula finds application in approximating discrete sums from lattice discretizations of path integrals in quantum mechanics and statistical mechanics, facilitating the transition to continuum limits for partition functions and effective actions.[38] For instance, it corrects Riemann sums in path integral formulations to yield accurate semiclassical approximations.[39]
Series expansions of functions
Bernoulli numbers play a central role in the Laurent and Taylor series expansions of trigonometric and hyperbolic functions, providing explicit coefficients that connect these expansions to number-theoretic properties. These series are particularly useful in analysis for representing functions in terms of powers of the variable, with convergence in disks avoiding the poles of the functions.[40]The Taylor series expansion of the tangent function around zero is given by\tan z = \sum_{n=1}^{\infty} \frac{(-1)^{n-1} 2^{2n} (2^{2n} - 1) B_{2n}}{(2n)!} z^{2n-1},valid for |z| < \pi/2. Here, the coefficients involve even-indexed Bernoulli numbers B_{2n}, which alternate in sign starting with B_2 = 1/6 > 0. This form arises from the relation between the tangent numbers and Bernoulli numbers, where the tangent numbers T_{2n-1} = (-1)^{n-1} 2^{2n} (2^{2n} - 1) B_{2n} appear when expressing the series as \tan z = \sum_{n=1}^{\infty} T_{2n-1} \frac{z^{2n-1}}{(2n-1)!}.[40]A closely related expansion is that of the cotangent function, which has a Laurent series due to its pole at zero:\cot z = \frac{1}{z} + \sum_{n=1}^{\infty} \frac{(-1)^n 2^{2n} B_{2n}}{(2n)!} z^{2n-1},converging for $0 < |z| < \pi. The principal term $1/z reflects the simple pole, and the remaining power series mirrors the structure of the tangent expansion but with adjusted signs. Scaling the argument yields the expansion for the scaled cotangent:\pi \cot(\pi z) = \frac{1}{z} + \sum_{n=1}^{\infty} (-1)^n \frac{(2\pi)^{2n} B_{2n}}{(2n)!} z^{2n-1},valid in $0 < |z| < 1. This series is derived from the generating function for Bernoulli numbers and connects to the Riemann zeta function for evaluating the coefficients, as B_{2n} = (-1)^{n-1} \frac{2 (2n)!}{(2\pi)^{2n}} \zeta(2n).[41]Hyperbolic analogs follow similarly, replacing trigonometric functions with their hyperbolic counterparts, which lack the alternating signs due to the absence of imaginary units. The Laurent series for the hyperbolic cotangent is\coth z = \frac{1}{z} + \sum_{n=1}^{\infty} \frac{2^{2n} B_{2n}}{(2n)!} z^{2n-1},converging for $0 < |z| < \pi i. This expansion is obtained by substituting z \to i z in the cotangent series and adjusting for the hyperbolic identity \coth z = i \cot(i z). The coefficients remain tied to Bernoulli numbers, facilitating computations in hyperbolic geometry and special functions.[42]In complex analysis, these series expansions underpin partial fraction decompositions of meromorphic functions. Notably, the expansion of \pi \cot(\pi z) equates to the infinite partial fraction sum \sum_{k=-\infty}^{\infty} \frac{1}{z + k}, providing a residue-based representation useful for contour integrals and summing series over integers. This connection highlights the Bernoulli numbers' role in bridging power series and residue calculus.[41]
Combinatorial Connections
Links to Stirling numbers
Bernoulli numbers exhibit a direct connection to Stirling numbers of the second kind through the explicit summation formulaB_n = \sum_{k=0}^n (-1)^k \frac{k!}{k+1} S(n,k),where S(n,k) denotes the Stirling number of the second kind, which counts the number of ways to partition n distinct objects into k non-empty unlabeled subsets.[43] This identity provides a combinatorial avenue for computing Bernoulli numbers, leveraging the recursive or explicit evaluations of S(n,k), and underscores the role of these numbers in bridging analysis and combinatorics.[43]Stirling numbers of the first kind, denoted by the signed s(n,k), which appear in the expansion of rising factorials and count permutations by cycle structure, connect to Bernoulli numbers via inversion relations inherent to the Stirling pair. A key identity is\sum_{i=k}^n \frac{S(n-1, i-1) s(i, k)}{i} = \frac{1}{n} \binom{n}{k} B_{n-k},illustrating how products of Stirling numbers from both kinds yield Bernoulli numbers scaled by binomial coefficients.[44] This relation arises from generating function manipulations and highlights the orthogonal structure between the two kinds of Stirling numbers in relation to Bernoulli sequences.[44]These links manifest prominently in the expression for sums of integer powers, where Stirling numbers of the second kind convert powers to falling factorials, whose sums are then captured by Bernoulli polynomials, providing an integrated framework for Faulhaber's formula.[7]
Relations to other combinatorial sequences
Bernoulli numbers exhibit connections to various combinatorial sequences beyond Stirling numbers, notably through identities involving Eulerian numbers, Worpitzky numbers, Pascal's triangle entries, and counts of alternating permutations.Worpitzky's identity establishes a fundamental link between powers, binomial coefficients, and Eulerian numbers A(n,k), which count permutations of n elements with exactly k-1 ascents:x^n = \sum_{k=0}^n A(n,k) \binom{x + k}{n}.This representation highlights the combinatorial role of Eulerian numbers in expressing monomials via binomial terms. Bernoulli numbers relate to this framework through generating functions for Eulerian polynomials and analogous identities for generalized Bernoulli numbers B_{n,k}. Specifically, a type-B variant of Worpitzky's identity is(2x + 1)^n = \sum_{k=0}^n B_{n,k} \binom{x + n - k}{n},where the B_{n,k} are type B Eulerian numbers. These connections underscore how Bernoulli numbers generalize the structure seen in Eulerian counts.The binomial coefficients in Pascal's triangle provide another combinatorial tie to Bernoulli numbers via their recursive definition. For m \geq 1, the Bernoulli numbers satisfy\sum_{k=0}^m \binom{m+1}{k} B_k = 0,with B_0 = 1. This equation interprets the Bernoulli numbers as weights that make the weighted sum of the (m+1)-th row of Pascal's triangle equal to zero, akin to a balancing condition or finite difference operator applied to the power basis. Such relations emphasize the role of binomial coefficients in generating and constraining the Bernoulli sequence combinatorially.[3]Alternating permutations, which are permutations \sigma of $$ satisfying \sigma(1) < \sigma(2) > \sigma(3) < \cdots (up-down type), connect directly to even-indexed Bernoulli numbers through Euler zigzag numbers. The number U_{2n-1} of such alternating permutations on $2n-1 elements equalsU_{2n-1} = \frac{2^{2n} (2^{2n} - 1)}{2n} |B_{2n}|.For example, when n=2, U_3 = 2, matching \frac{16 \cdot 15}{4} \cdot \frac{1}{30} = 2. This formula provides a combinatorial interpretation of |B_{2n}| as scaling the count of these permutations, linking analytic properties of Bernoulli numbers to permutation statistics. Eulerian numbers relate indirectly here, as they enumerate all permutations by ascent count, while alternating permutations represent extremal cases with maximally alternating ascent-descent patterns.[45]
Advanced Arithmetical Properties
Von Staudt–Clausen theorem
The von Staudt–Clausen theorem provides a precise description of the fractional part of even-indexed Bernoulli numbers B_{2n} for n \geq 1. Specifically, it states thatB_{2n} + \sum_{\substack{p \text{ prime} \\ (p-1) \mid 2n}} \frac{1}{p} \in \mathbb{Z},where the sum is taken over all primes p such that p-1 divides $2n. This theorem was independently discovered in 1840 by Karl Georg Christian von Staudt and Thomas Clausen, with von Staudt providing the full result on the denominators and Clausen contributing to the fractional part characterization.[46][3]A proof of the theorem can be sketched using the Hurwitz zeta function \zeta(s, a) = \sum_{k=0}^\infty (k+a)^{-s} for \Re(s) > 1 and its analytic continuation, which relates to Bernoulli polynomials via \zeta(1-m, a) = -\frac{B_m(a)}{m} for positive integers m. The connection arises from the partial fraction expansion of \pi \cot(\pi z), whose coefficients involve sums over residues that yield the Bernoulli numbers through the Euler-Maclaurin summation formula applied to the zeta function. By examining the poles and residues at integer points, or equivalently, the behavior of the Hurwitz zeta at even negative integers, one derives the required sum over reciprocals of primes dividing the denominator after clearing fractional parts via localization at each prime. This approach confirms the integer nature of the expression by showing that the denominator divides the product of relevant primes.[47]The theorem implies that the denominator of B_{2n} (in lowest terms) is exactly the product of all distinct primes p such that p-1 divides $2n, making it square-free and always divisible by 6 for n \geq 1. For example, when n=1, $2n=2, the primes are p=2 and p=3 (since $1 \mid 2 and $2 \mid 2), so B_2 + \frac{1}{2} + \frac{1}{3} = 1 \in \mathbb{Z} and the denominator of B_2 = \frac{1}{6} is $2 \times 3 = 6. Similarly, for n=2, $2n=4, the primes are p=2,3,5, yielding denominator 30 for B_4 = -\frac{1}{30}. This structure highlights the arithmetic purity of Bernoulli numbers, with no higher prime powers in the denominators.[46][47]
Congruences and p-adic aspects
Bernoulli numbers satisfy several important modular congruences that reveal their arithmetic structure, particularly in relation to prime moduli. Kummer congruences provide a foundational set of relations for these numbers. Specifically, for an odd prime p and integers m, n > 1 with m \equiv n \pmod{p-1}, the normalized Bernoulli numbers satisfy \frac{B_m}{m} \equiv \frac{B_n}{n} \pmod{p}. These congruences extend to higher powers, stating that \frac{B_{m + k(p-1)}}{m + k(p-1)} \equiv \frac{B_n}{n} \pmod{p} for suitable k, demonstrating a form of periodicity modulo p. Kummer's work established these as key tools for analyzing the p-adic properties of Bernoulli numbers.[48]In the p-adic setting, these congruences underpin the continuity of Bernoulli numbers. The sequence \frac{B_n}{n} for even n congruent to 0 modulo p-1 converges p-adically, allowing interpolation via the p-adic zetafunction introduced by Kubota and Leopoldt. This function extends the Riemann zeta function p-adically and interpolates values related to Bernoulli numbers at negative integers, ensuring \zeta_p(1 - n) = -\frac{B_n}{n} for n not divisible by p-1. For irregular primes—those where p divides the numerator of \frac{B_k}{k} for some even k = 2, 4, \dots, p-3—the p-adic valuation of these numerators can exceed 1, leading to B_n \equiv 0 \pmod{p^r} for higher r in specific arithmetic progressions. Such irregularity measures the deviation from p-adic continuity and relates to zeros of the p-adic zetafunction.[48][49]These arithmetical properties have profound applications, notably in the study of Fermat's Last Theorem. Kummer's criterion links the regularity of an odd prime p to the class number h of the cyclotomic field \mathbb{Q}(\zeta_p): p is regular if p \nmid h, which is equivalent to p not dividing the numerator of \frac{B_k}{k} for even k = 2, 4, \dots, p-3. For regular primes, Kummer proved that there are no non-trivial solutions to x^p + y^p = z^p in integers, establishing Fermat's Last Theorem for infinitely many exponents. Irregular primes, tied to higher p-adic divisibility of Bernoulli numerators, obstruct this proof but inform the distribution of solvable cases.[50]
Vanishing of odd Bernoulli numbers
The Bernoulli numbers B_n for odd indices n = 2k+1 with k \geq 1 vanish, i.e., B_{2k+1} = 0.[51] This property, except for the case n=1 where B_1 = -\frac{1}{2}, arises directly from the symmetry of their generating function.[51]The exponential generating function for the Bernoulli numbers is given by\frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!},valid for |x| < 2\pi.[51] To see the vanishing, consider the modified functionh(x) = \frac{x}{e^x - 1} + \frac{1}{2} x.This function is even, as h(-x) = h(x), because substituting -x yields\frac{-x}{e^{-x} - 1} + \frac{1}{2} (-x) = \frac{x e^x}{e^x - 1} - \frac{1}{2} x = x + \frac{x}{e^x - 1} - \frac{1}{2} x = h(x).[51] The evenness of h(x) implies that its Taylor series expansion around x = 0 contains only even powers of x. Since h(x) differs from the generating function by the linear term \frac{1}{2} x, and given that B_1 = -\frac{1}{2} makes the coefficient of x in h(x) zero, all higher odd-powered coefficients in the expansion of \frac{x}{e^x - 1} must vanish.[51]This vanishing has significant implications in analytic number theory, particularly for the Riemann zeta function \zeta(s). The values at negative integers are related by \zeta(1 - k) = -\frac{B_k}{k} for integers k \geq 2.[52] For k = 2m+1 odd and greater than 1, this yields \zeta(-2m) = 0 for positive integers m, corresponding to the trivial zeros of \zeta(s) at the negative even integers s = -2, -4, -6, \dots.[52]Historically, Leonhard Euler introduced the modern generating function approach to Bernoulli numbers in the mid-18th century, resolving earlier puzzles about their patterns observed by Jakob Bernoulli in sums of powers.[3] Euler's analysis highlighted the exceptional non-vanishing of B_1 and the systematic zeroing of higher odd indices, providing clarity to their otherwise irregular appearance in series expansions.[3]
Generalizations and Extensions
Generalized Bernoulli numbers
Generalized Bernoulli numbers B_{n,\chi} extend the classical Bernoulli numbers to Dirichlet characters \chi modulo a conductor f, playing a key role in analytic number theory, particularly for evaluating L-functions at integer points. They are defined byB_{n,\chi} = f^{n-1} \sum_{a=1}^f \chi(a) B_n\left( \frac{a}{f} \right),where B_n(x) are the Bernoulli polynomials, or equivalently via the generating function\sum_{n=0}^\infty B_{n,\chi} \frac{t^n}{n!} = \frac{t e^{x t}}{e^t - 1} \sum_{a=1}^f \overline{\chi}(a) e^{2\pi i a x / f}for a suitable x. When \chi is the trivial character, they reduce to the ordinary Bernoulli numbers up to scaling.[53]These numbers satisfy properties analogous to the classical case, such as vanishing for odd n > 1 under the trivial character, but more generally relate to the arithmetic of the character. A prominent application is in the class number formula for imaginary quadratic fields, where for non-principal characters, h(-d) = \frac{w \sqrt{d}}{2\pi} |L(1,\chi_d)| with L(1,\chi_d) = -\frac{\pi}{f \sqrt{f}} B_{1,\chi_d} for odd primitive \chi_d. They also appear in the von Staudt–Clausen theorem generalizations and Kummer congruences for irregular primes.[54] Higher-order generalizations, such as Bernoulli numbers of order l, are defined via \left( \frac{t}{e^t - 1} \right)^l = \sum_{n=0}^\infty B_n^{(l)} \frac{t^n}{n!}, reducing to classical for l=1.[55]
Bernoulli polynomials
The Bernoulli polynomials B_n(x) form a sequence of monic polynomials of degree n that extend the concept of Bernoulli numbers to a variable argument, with the generating function given by\frac{t e^{x t}}{e^t - 1} = \sum_{n=0}^\infty B_n(x) \frac{t^n}{n!}.This generating function encapsulates their periodic and analytic properties, distinguishing them from the constant Bernoulli numbers evaluated at x = 0. They are defined explicitly byB_n(x) = \sum_{k=0}^n \binom{n}{k} B_k x^{n-k},ensuring B_n(0) = B_n.[56]A fundamental property is the differencerecurrence relationB_n(x+1) - B_n(x) = n x^{n-1},which holds for all nonnegative integers n and real x, reflecting their role in difference equations and summation techniques. This relation arises directly from differentiating the generating function with respect to x and comparing coefficients. For illustration, the first few polynomials are B_0(x) = 1, B_1(x) = x - \frac{1}{2}, and B_2(x) = x^2 - x + \frac{1}{6}. The periodic extensions \widetilde{B}_n(x) = B_n(\{x\}), where \{x\} is the fractional part of x, further extend this utility in Fourier analysis and distribution theory.[56]On the unit interval [0,1], the Bernoulli polynomials possess a Fourier series representation, which connects them to exponential sums and zeta functions. For n \geq 2,B_n(x) = -\frac{n!}{(2\pi i)^n} \sum_{\substack{k=-\infty \\ k \neq 0}}^\infty \frac{e^{2\pi i k x}}{k^n},with convergence holding uniformly for n > 1 and pointwise otherwise, excluding endpoints for odd n \geq 1. Equivalent real forms separate into cosine series for even n and sine series for odd n, aiding evaluations in harmonic analysis. For example, when n=2,B_2(x) = \frac{1}{\pi^2} \sum_{k=1}^\infty \frac{\cos(2\pi k x)}{k^2},illustrating the series' explicit form. These expansions are instrumental in proving asymptotic behaviors and irregularities in the distribution of primes via the Riemann zeta function.[56]Bernoulli polynomials exhibit notable integral properties, such as\int_0^1 B_n(x) \, dx = 0for all n \geq 1, which follows from their construction via successive integration starting from B_0(x) = 1 and the choice of integration constants to enforce this orthogonality. This vanishing integral underscores their utility in approximating integrals by sums, particularly in the Euler-Maclaurin formula, where the polynomials (or their periodic extensions) provide the higher-order correction terms that bridge discrete sums and continuous integrals without deriving the full expansion here.[57]Furthermore, Bernoulli polynomials connect to multiple zeta values through evaluations at non-positive integers; specifically, multiple zeta functions \zeta(s_1, -m_2, -m_3) for suitable arguments can be expressed as linear combinations involving Bernoulli polynomials like B_3(s_1 - 2) and products of lower-degree terms, enabling analytic continuations and explicit computations in number theory. Variants such as generalized Bernoulli numbers arise from evaluating these polynomials at rational points, for instance relating to B_n(1-x) under symmetry relations.[58][56]
Other Representations and Views
Tree and algorithmic representations
The Seidel triangle provides a tabular method for computing Bernoulli numbers, analogous to Pascal's triangle but incorporating alternating sums to generate the sequence. Introduced by Ludwig Seidel in 1877, this construction begins with initial entries related to Genocchi numbers or simple integers and builds subsequent rows through a recursive relation where each entry is the sum of entries from the previous row, often with sign alternations to reflect the boustrophedon (ox-plowing) pattern. Specifically, for the Euler-Seidel matrix with entries a_k^n, the relation is a_k^n = a_{k-1}^n + a_{k-1}^{n+1} for n > 0, k > 1, with initial conditions such as a_0^n = B_n for n \geq 0 (adjusted for conventions), leading to rows whose alternating partial sums yield the Bernoulli numbers B_n. For instance, the alternating sum of the n-th row produces B_n / n! up to sign conventions, enabling systematic computation without direct use of the standard recursive definition involving binomial coefficients.[59]This triangle's structure facilitates visualization of the dependencies in Bernoulli number calculations, highlighting how higher-order terms emerge from lower ones through linear combinations. Seidel's original work, published in the Sitzungsberichte der Mathematisch-Physikalischen Classe der Königl. Bayer. Akademie der Wissenschaften, detailed this "Treppen-Schema" (staircase scheme) as an efficient unfolding of the numbers from elementary components, predating modern algorithmic perspectives but proving foundational for subsequent combinatorial analyses.[59]A binary tree representation offers another graphical and recursive framework for generating Bernoulli numbers, mirroring the recursive nature of their definition through branched partial sums. As described by Woon in 1998, the tree is constructed using left and right operators O_L and O_R, where O_L shifts factorial denominators and O_R introduces a factor related to 2! with sign flip, starting from the seed $1/2!. Each node represents a partial sum term of the form \pm 1/(a! b! \cdots), and the n-th Bernoulli number is obtained as B_n = n! times the sum of terms in the (n-1)-th level, providing a tree depth proportional to n with $2^{n-1} leaves. This structure visualizes the expansion of the generating function.[60]
Combinatorial interpretations via permutations
The combinatorial interpretations of Bernoulli numbers via permutations primarily arise through their connections to alternating (or up-down) permutations, where the counts of such permutations yield the Euler numbers, which in turn relate directly to the Bernoulli numbers via explicit formulas.In 1879, Désiré André established that the number of alternating permutations of the set {1, 2, ..., m}—permutations π satisfying π(1) < π(2) > π(3) < π(4) > ⋯ or the down-up variant π(1) > π(2) < π(3) > π(4) < ⋯—is given by the Euler number E_m, with the exponential generating function ∑{m≥0} E_m x^m / m! = sec x + tan x.[61] The even-indexed Euler numbers E{2n} are known as secant numbers, while the odd-indexed ones E_{2n-1} (starting from n=1) are the tangent numbers T_n, which count the alternating permutations of 2n-1 elements. These tangent numbers provide a direct link to the Bernoulli numbers through the formulaT_n = \frac{2^{2n} (2^{2n} - 1) |B_{2n}|}{2n},where B_{2n} denotes the 2n-th Bernoulli number.[62] Rearranging yields|B_{2n}| = \frac{2n \cdot T_n}{2^{2n} (2^{2n} - 1)} = \frac{2n \cdot E_{2n-1}}{2^{2n} (2^{2n} - 1)},expressing the magnitude of even Bernoulli numbers in terms of the count of alternating permutations. Similar relations hold for secant numbers, stemming from comparing the Taylor series expansions of tan x and sec x with the known series involving Bernoulli numbers.Entringer numbers offer a refinement of this interpretation by classifying alternating permutations according to their initial descent. Specifically, the Entringer number E(r, s) counts the number of down-up alternating permutations of {1, 2, ..., r+1} that begin with the value s+1. The row sums ∑s E(r, s) = E{r+1} recover the Euler numbers, forming the Seidel-Entringer-Arnold triangle. In a seminal 1966 paper, Richard Entringer provided a combinatorial interpretation of the Bernoulli numbers themselves using these refined counts of alternating permutations, linking the rational values of B_m to signed or weighted enumerations within the triangle.[63] For example, certain alternating sums over entries in specific rows or diagonals of the Entringer triangle yield expressions involving Bernoulli numbers, extending the unsigned counts of André to the signed rational structure of the Bernoulli sequence.[64]Bijective proofs of these enumerations often employ tree structures to establish equivalences. For instance, there are natural bijections between alternating permutations of m elements and certain classes of increasing trees, such as even trees with m nodes, both enumerated by E_m. This bijection, constructed via recursive insertion or growth rules on the trees corresponding to the permutation's up-down pattern, provides a combinatorial verification of André's theorem and, by extension, the relations to Bernoulli numbers. Such tree-based bijections facilitate proofs of identities linking Euler and Bernoulli numbers without relying on generating functions, by directly mapping permutation cycles or descents to tree branches.[65]Alternating permutations also connect briefly to Eulerian numbers, which count permutations by the number of ascents; generating functions for certain statistics on alternating permutations incorporate Eulerian numbers as coefficients in expansions related to the Euler-Bernoulli identities.[61]