1 Introduction

1.1 Motivation: from sums to integrals

Many problems in analysis and computation require estimating a finite sum of the form \[ S=\sum_{k=m}^{n} f(k), \] where \(f\) is a smooth function on an interval large enough to contain the integers in question. Direct evaluation can be expensive, while integrals are often easier to approximate. The Euler–Maclaurin summation formula provides a systematic bridge between these two worlds: it expresses \(S\) as an integral plus explicit correction terms that depend on derivatives of \(f\) at the endpoints.

1.2 Historical context and origin

The formula is attributed to Leonhard Euler and Colin Maclaurin. Euler developed techniques that effectively relate discrete and continuous behavior through expansions involving boundary data. Maclaurin refined and generalized such ideas, leading to the modern statement in which Bernoulli numbers appear naturally. The result has since become a standard tool in asymptotic analysis, numerical methods, and analytic estimates.

1.3 Basic notation and assumptions on the function

A typical statement uses a function \(f\) that is sufficiently differentiable on \([m,n]\). In many applications, one assumes \(f\) has derivatives up to some order \(2p\) that are continuous on the closed interval, and that higher derivatives behave reasonably (e.g., integrable or of controlled growth). Under such hypotheses, the formula yields an asymptotic series in powers of the endpoint spacing.

2 The Core Formula

2.1 Standard finite-interval version

A common finite-interval form (with step size 1) is \[ \sum_{k=m}^{n} f(k) = \int_{m}^{n} f(x)\,dx +\frac{f(m)+f(n)}{2} +\sum_{r=1}^{p-1}\frac{B_{2r}}{(2r)!}\Big(f^{(2r-1)}(n)-f^{(2r-1)}(m)\Big) +R_p, \] where \(B_{2r}\) are Bernoulli numbers and \(R_p\) is a remainder term depending on derivatives of order \(2p\) (or higher) of \(f\).

2.2 Infinite/limiting forms (when applicable)

If one or both endpoints extend to infinity, or if the sum is over a semi-infinite range, additional conditions are needed to ensure convergence and to control boundary behavior. Under decay or analyticity assumptions on \(f\), the formula can be applied with \(m\to -\infty\) or \(n\to +\infty\), causing boundary contributions to simplify or vanish.

2.3 Role of Bernoulli numbers and Bernoulli polynomials

Bernoulli numbers appear in the correction coefficients in the unshifted version. In shifted variants—where the sum is taken over \(k+a\) rather than \(k\)—Bernoulli polynomials \(B_r(a)\) replace the corresponding numbers. This is also how periodic Bernoulli functions enter the derivation: the coefficients encode the behavior of certain periodic functions related to \(B_r\).

2.4 Boundary terms from endpoint derivatives

The correction terms depend only on derivatives evaluated at the endpoints. For each even index \(2r\), one uses the \((2r-1)\)-st derivative at \(m\) and \(n\), combined with the difference \(f^{(2r-1)}(n)-f^{(2r-1)}(m)\). This “endpoint-only” structure reflects that the discrepancy between discrete sampling and continuous integration is largely concentrated near boundaries after repeated smoothing.

2.5 Error term and its interpretation

The remainder \(R_p\) can be written in several equivalent ways. A typical interpretation is that truncating after \(p-1\) correction terms leaves an error governed by the next higher derivative(s) of \(f\). In many asymptotic settings, \(R_p\) is bounded by a quantity involving \(\|f^{(2p)}\|\) and a factor related to the scale of the interval. In other settings, \(R_p\) is itself asymptotic to a term that reveals how quickly the approximation improves as the interval grows or as more terms are included.

3 Derivation and Intuition

3.1 Integral remainder viewpoint

A useful intuition is to compare the sum with an integral over the same range by “filling in” each unit interval \([k,k+1]\). One can write \[ \sum_{k=m}^{n-1} f(k) \approx \sum_{k=m}^{n-1}\int_{k}^{k+1} f(x)\,dx, \] then quantify the difference by repeatedly expanding \(f(x)\) around integer points and collecting the resulting boundary contributions. Each level of expansion produces additional derivative terms and ultimately yields the Bernoulli-number coefficients.

3.2 Connection to the periodic Bernoulli function

The coefficients in the Euler–Maclaurin formula are closely linked to Bernoulli polynomials evaluated at fractional parts. In derivations, one introduces a periodic Bernoulli function—essentially a periodic version of \(B_r(\{x\})\)—whose Fourier series leads to the same constants \(B_{2r}\). The periodic function captures how the “sampling grid” deviates from uniform integration.

3.3 Using the Euler transform and finite differences

Another route uses finite differences. The forward difference operator \(\Delta f(k)=f(k+1)-f(k)\) and its iterates relate to derivatives through smoothness assumptions. By expressing the sum in terms of repeated differences and then connecting those differences to integrals of derivatives, one arrives at an expansion resembling Euler–Maclaurin. This perspective explains why higher derivatives appear with alternating parity: the structure of the difference operator naturally filters certain derivative orders.

3.4 Sketch of proof via contour/analytic techniques (overview)

For sufficiently analytic functions, complex analysis provides an alternative derivation. One constructs an integral representation of the sum using a kernel with known poles at integers (or a related meromorphic function). Shifting contours and summing residues yields an expansion whose coefficients match Bernoulli numbers. This method also offers a clear route to remainder estimates under analytic continuation and growth conditions, though the full technical details depend on the chosen kernel and function class.

4 Variants and Generalizations

4.1 Different normalization conventions

Different sources use distinct indexing, endpoint choices, and conventions for the remainder term. Some formulas include the integral from \(m\) to \(n+1\) rather than \(m\) to \(n\), or use \(\sum_{k=m}^{n-1}\) rather than \(\sum_{k=m}^{n}\). These changes shift the half-endpoint term and may alter the sign pattern in the derivative differences, while leaving the underlying coefficient structure governed by Bernoulli numbers.

4.2 Symmetric (midpoint) formulations

A symmetric form uses midpoint information. Instead of treating endpoints \(m\) and \(n\) asymmetrically, one may rewrite the sum in terms of derivatives evaluated at symmetric points around the interval center. This can be advantageous when \(f\) has special symmetry (e.g., even or odd behavior) that forces cancellation of certain terms.

4.3 Versions with shifts and starting indices

When the summand is evaluated at \(k+a\) for a fixed shift \(a\), Bernoulli polynomials enter: \[ \sum_{k=m}^{n} f(k+a) = \int_{m}^{n} f(x+a)\,dx +\text{(boundary corrections involving } B_r(a)\text{)} +\cdots. \] The shift changes the “fractional part” statistics of the periodic Bernoulli function, altering coefficients accordingly.

4.4 Application to alternating sums

Alternating sums such as \(\sum (-1)^k f(k)\) can be handled by incorporating a factor that changes the sampling pattern. Conceptually, one can apply Euler–Maclaurin to a transformed function (or relate the alternating sum to sums over even and odd indices). The resulting expansion typically uses Bernoulli numbers modified by the alternating behavior, sometimes leading to suppression of specific terms.

4.5 Multidimensional extensions (overview)

There are generalizations to sums over integer lattices, where corrections involve mixed derivatives and higher-dimensional analogues of Bernoulli data. In such settings, the structure depends on the geometry of the summation region (e.g., boxes versus more complex polytopes). The core idea remains: the discrete-to-continuum discrepancy is captured by boundary integrals and derivative corrections along faces, edges, and vertices.

5 Practical Computation

5.1 How to choose the number of correction terms

The formula yields an asymptotic expansion in many contexts. Practically, one selects \(p\) so that the computed correction terms improve the approximation without amplifying numerical errors. For analytic functions and large intervals, more terms often reduce the truncation error up to a point, after which the asymptotic nature may lead to deterioration.

5.2 Truncation error estimates

Remainder estimates depend on what is known about derivatives of \(f\). A typical strategy is:

1 Introduction

2 The Core Formula

3 Derivation and Intuition

For many smooth, slowly varying functions, including a moderate number of Bernoulli corrections provides rapid improvement.

5.3 Smoothness requirements and common sufficient conditions

Euler–Maclaurin requires enough derivatives for the terms included. For finite-interval applications, one commonly assumes \(f\in C^{2p}\) (continuously differentiable up to order \(2p\)) on \([m,n]\), with an integrable or bounded \(f^{(2p)}\). When the interval endpoints are moved outward (for asymptotics), additional decay or controlled growth ensures boundary contributions remain finite.

5.4 Handling endpoints and singular behavior

If \(f\) has singularities at or near endpoints, the standard finite-interval formula may fail or require modification (e.g., splitting the integral away from singular points). In computations, it is typical to:

  • isolate regions where derivatives blow up,
  • use local expansions matched to the singular behavior,
  • apply Euler–Maclaurin only on subintervals where smoothness holds.

6 Applications

6.1 Asymptotic expansions of partial sums

For sequences \(S(N)=\sum_{k=1}^{N} f(k)\) with \(N\to\infty\), Euler–Maclaurin converts \(S(N)\) into a combination of:

  • an integral term growing like the leading order behavior,
  • a series of inverse powers (through derivatives and Bernoulli numbers),
  • a controlled remainder.

This is a standard mechanism for deriving asymptotic expansions of partial sums.

6.2 Numerical summation and acceleration

In numerical analysis, Euler–Maclaurin can accelerate convergence of series. For example, slowly decaying terms can be integrated, while correction terms account for the discrete mismatch. The result is an approximation that can outperform direct summation, particularly when high derivatives can be computed reliably.

6.3 Asymptotics in special functions

Many special functions admit definitions through sums that are amenable to Euler–Maclaurin expansion. By transforming the sum representation, one obtains asymptotic expansions of functions such as logarithms of gamma-type functions, related zeta-like quantities, or other series-defined objects. The formula thus provides a systematic route from discrete definitions to continuous asymptotic behavior.

6.4 Analytic number theory (typical use cases, non-controversial)

In analytic number theory, Euler–Maclaurin is often used to estimate sums over arithmetic functions by replacing them with integrals plus correction terms. Even when the summand itself is not elementary, one may obtain approximate forms by inserting smooth weight functions and then applying the formula to the resulting smooth expressions. This yields asymptotic estimates with explicit error terms, a key ingredient in many analytic arguments.

7 Worked Examples

7.1 Polynomial test functions

For \(f(x)=x^p\), the sum can be computed exactly by Faulhaber’s theorem, and Euler–Maclaurin reproduces these polynomial identities after a finite number of correction terms. Because derivatives of a polynomial eventually vanish, the remainder term becomes zero once the differentiation order exceeds the polynomial degree. This makes the formula a convenient check of coefficients and signs.

7.2 Exponential and trigonometric functions

For \(f(x)=e^{ax}\) (with suitable \(a\)) or \(f(x)=\sin(ax),\cos(ax)\), derivatives remain simple multiples of the original function. Euler–Maclaurin then produces an expansion whose terms involve the same exponential or trigonometric factors evaluated at endpoints. Such examples illustrate how quickly correction terms shrink when \(f\) varies smoothly.

7.3 Rational functions with controlled derivatives

Consider functions like \(f(x)=1/(x+b)\) or \(f(x)=1/(x^2+b^2)\) with \(b>0\). While the function is smooth on \([m,n]\) away from singularities, derivatives can grow in magnitude as \(x\) approaches a pole. Choosing \(m\) and \(n\) to keep denominators safely bounded allows Euler–Maclaurin to yield accurate estimates. The example demonstrates the importance of smoothness and derivative control near endpoints.

7.4 Comparing exact sums vs. Euler–Maclaurin approximations

A standard comparison approach is:

  • compute the sum exactly for moderate \(m,n\) (or using high-precision arithmetic),
  • compute the Euler–Maclaurin approximation for increasing \(p\),
  • observe convergence trends and how the remainder estimate correlates with observed error.

Such comparisons highlight that the series is typically asymptotic: improving \(p\) decreases error until an optimal truncation point is reached, after which numerical or asymptotic effects dominate.

8 Relations to Other Results

Euler–Maclaurin can be viewed as a systematic refinement of the trapezoidal rule. The leading correction \(\tfrac{f(m)+f(n)}{2}\) resembles the end-point averaging inherent in trapezoidal integration. Higher-order derivative terms correspond to error corrections for the trapezoidal rule when the integrand is smooth, thereby linking discrete quadrature error analysis with Bernoulli-number expansions.

8.2 Relation to Poisson summation (conceptual comparison)

Poisson summation relates sums over integer points to sums of Fourier transforms. While Euler–Maclaurin focuses on local smoothness and endpoint derivatives to approximate a sum by an integral plus corrections, Poisson summation captures global frequency information through Fourier analysis. Both methods connect discrete and continuous behavior, but they emphasize different aspects: one is derivative-driven and local, the other is transform-driven and spectral.

8.3 Relation to Faulhaber-type formulas (overview)

Faulhaber formulas express \(\sum_{k=1}^{N} k^p\) as polynomials in \(N\). Euler–Maclaurin provides a general framework that, for polynomial \(f\), reduces to an exact finite expression. In this sense, Faulhaber-type identities can be obtained as special cases of Euler–Maclaurin when the remainder term vanishes.

9 Bernoulli Material

9.1 Bernoulli numbers: definitions and properties

Bernoulli numbers \(B_n\) are a sequence of rational numbers defined by generating functions and satisfy many algebraic identities. Key properties include:

  • \(B_{2n+1}=0\) for \(n\ge 1\),
  • their appearance as coefficients in expansions of \(\frac{x}{e^x-1}\),
  • relationships that enable simplification of Euler–Maclaurin terms to even indices only.

9.2 Bernoulli polynomials

Bernoulli polynomials \(B_n(x)\) generalize Bernoulli numbers via evaluation at \(x\). They satisfy symmetry and differentiation identities, and they naturally encode endpoint effects in shifted Euler–Maclaurin formulas. Setting \(x=0\) recovers the Bernoulli numbers: \(B_n(0)=B_n\).

9.3 Useful identities and recursion relations

Bernoulli numbers and polynomials satisfy recursion relations and difference formulas, such as relationships between \(B_n(x+1)-B_n(x)\) and lower-degree polynomials. These identities are used to transform correction terms across different conventions and to connect coefficients appearing in various variants of the summation formula.

9.4 Generating functions (compact overview)

Generating functions provide compact characterizations. For Bernoulli polynomials, one uses a generating series whose coefficients are \(B_n(x)\), and for Bernoulli numbers one obtains the special case at \(x=0\). In derivations of Euler–Maclaurin, these generating functions link Bernoulli data to periodic kernels and Fourier expansions.

10 Error Bounds and Conditions

10.1 Remainder expressed with a derivative integral

A common remainder representation involves an integral over the interval of a derivative of \(f\) multiplied by a bounded Bernoulli-related kernel. This form makes clear that smoothness of \(f\) of order \(2p\) (and suitable integrability of that derivative) directly controls the magnitude of the remainder.

10.2 Big-O and asymptotic remainder conventions

Different authors express remainders using either a big-\(O\) bound or an asymptotic equivalence. For large parameter regimes—such as sums up to \(N\)—one often writes \(R_p=O(N^{-2p+1})\) under appropriate decay assumptions. Alternatively, one may record that the remainder is asymptotic to a constant multiple of a higher derivative scale, reflecting the asymptotic nature of the correction series.

10.3 Optimal smoothness and decay assumptions (overview)

The “right” assumptions depend on the intended limit. On finite intervals, continuity of derivatives up to order \(2p\) is often enough for a valid expansion with a controlled remainder. In limits where endpoints grow, additional decay (or analyticity) can ensure that boundary corrections stabilize and that the remainder becomes asymptotically small. Without such conditions, the boundary contributions may dominate or the correction series may fail to approximate the sum reliably.

11 Implementation Notes

11.1 Symbolic vs. numerical use

In symbolic computations, Euler–Maclaurin can produce explicit expansions in terms of Bernoulli numbers and derivatives, which may then be simplified algebraically. In numerical work, one must evaluate derivatives at endpoints accurately and guard against cancellation, especially when endpoint derivatives are large but alternately signed.

11.2 Common pitfalls in truncation and scaling

Two typical pitfalls are:

  • truncating beyond the range where the asymptotic series is beneficial,
  • using inconsistent scaling when formulas are written for step size 1 but the problem effectively uses a different spacing.

Another issue arises when the function changes rapidly near endpoints: the derivatives required for high-order corrections can become numerically unstable.

11.3 Stability considerations for computer algebra systems

Computer algebra systems can symbolically compute derivatives and Bernoulli numbers exactly, but numerical evaluation may introduce rounding errors. A stable implementation often:

  • selects a truncation order based on remainder estimates,
  • uses high precision for endpoint derivatives,
  • checks sensitivity by comparing partial sums with different truncation depths.