A Simpson–Type Decomposition of the Euler–Mascheroni Constant

1. Introduction

The Euler–Mascheroni constant is a classical object of analysis and number theory [1]. Let

$$H\left(N\right)=\sum _{k=1}^N{\displaystyle \frac{1}{k}}$$

(1)

denote the N-th harmonic number. Then the Euler–Mascheroni constant can be written succinctly as

$$\gamma =\underset{N\to \infty }{lim}\left(H\left(N\right)-log\left[N\right]\right)\approx 0.57721\dots .$$

(2)

Figure 1 illustrates the difference between the integral of $1/x$ and its approximation by the rectangular rule. The existence of $\gamma$ is usually established via asymptotic expansions derived from the Euler–-Maclaurin summation formula or related analytic techniques [2]. While this limit is familiar, most standard proofs rely on the Euler–Maclaurin summation formula, and it is natural to ask whether $\gamma$ can be understood using only elementary numerical ideas.

According to Havil [3], no proof of the irrationality of Euler’s constant $\gamma$ is currently known, and the problem remains open. Moreover, the classical limit defining $\gamma$ is known to converge very slowly, a fact that motivates the search for alternative approximations with provably faster convergence. In this paper, we study a class of rational approximations arising from Simpson–type quadrature formulas, and derive explicit error bounds that clarify the rate of convergence of the resulting sequences.

2. Analytic Approximation of the Logarithmic Integral

2.1. Simpson–Type Approximation of Logarithmic Increments

Definition 1. 

For integers $n\ge 1$, define

$$f\left(n\right){\int }_{2n-1}^{2n+1}{\displaystyle \frac{1}{x}}dx=log[2n+1]-log[2n-1].$$

(3)

This integral is approximated using Simpson’s rule with a step size of 1, as illustrated in Figure 2. Using the local Simpson weight

$$g\left(n\right)={\displaystyle \frac{1}{3}}\left({\displaystyle \frac{1}{2n-1}}+{\displaystyle \frac{4}{2n}}+{\displaystyle \frac{1}{2n+1}}\right),$$

(4)

we now define a global Simpson-regularized approximation of the harmonic sum. Notably, $g\left(n\right)$ is a finite linear combination of reciprocal integers and depends solely on the values of $1/x$ at the integer points.

2.2. Convexity and Comparison Inequalities

The local quadrature error is defined as follows:

$$d\left(n\right)=g\left(n\right)-f\left(n\right).$$

(5)

The function $1/x$ is strictly decreasing and convex on $(0,\infty )$. Let $p_n\left(x\right)$ denote a quadratic polynomial interpolating $1/x$ at the three points.

$$x=2n-1,2n,2n+1.$$

(6)

By construction,

$$g\left(n\right)={\int }_{2n-1}^{2n+1}{p}_n\left(x\right)dx.$$

(7)

As $1/x$ is convex, the interpolation polynomial $p_n\left(x\right)$ lies above $1/x$ on $[2n-1,2n+1]$. Thus, $p_n\left(x\right)-f\left(x\right)>0$ except at the midpoint, and integration yields $d\left(n\right)>0$. Consequently,

$$d\left(n\right)>0\mathrm{for}n\ge 1.$$

(8)

Moreover, convexity implies a comparison across adjacent intervals; in particular,

$$f\left(n\right)>g(n+1)\mathrm{for}n\ge 2.$$

(9)

It follows that

$$0<d\left(n\right)<g\left(n\right)-g(n+1)\mathrm{for}n\ge 2.$$

(10)

2.3. Telescoping Bounds and Error Convergence

Summing the aforementioned inequalities for $n=2,\dots ,N$ yields

$$0<\sum _{n=2}^{N}d\left(n\right)<\sum _{n=2}^N\{g\left(n\right)-g(n+1)\}.$$

(11)

Both bounds telescope:

$$0<\sum _{n=2}^{N}d\left(n\right)<g\left(2\right)-g(N+1)$$

(12)

As $g\left(n\right)\to 0$ as $n\to \infty$ and $d\left(n\right)>0$ for $n\ge 2$, the partial sums

$$\sum _{n=2}^{N}d\left(n\right)$$

(13)

form a monotone increasing sequence. Moreover, since

$$\sum _{n=2}^{N}d\left(n\right)<g\left(2\right)={\displaystyle \frac{23}{15}}\approx 1.5333\dots ,$$

(14)

the sequence of partial sums is bounded above and therefore converges.

Definition 2 

(The constant $\delta$). The constant δ is naturally defined by

$$\delta \underset{N\to \infty }{lim}\sum _{n=2}^{N}d\left(n\right)=\underset{N\to \infty }{lim}\sum _{n=2}^N\{g\left(n\right)-f\left(n\right)\}.$$

(15)

3. Arithmetic Consequences of Simpson-Regularized Sums

We now assemble the local Simpson approximations introduced above to obtain a global regularized analogue of the harmonic sum.

3.1. Recovery of the Euler–Mascheroni Constant

Definition 3 

(Simpson-regularized harmonic sum). Let N be a positive integer.

(i) N odd. If N is odd, we partition the interval $[1,N]$ into subintervals of length two,

$$[1,3],[3,5],\dots ,[N-2,N],$$

(16)

and apply Simpson’s rule to each subinterval to approximate the integral

$${\int }_1^N{\displaystyle \frac{1}{x}}dx.$$

(17)

The resulting approximation is denoted by $h\left(N\right)$ and is given by

$$h\left(N\right)\sum _{k=1}^{{\displaystyle \frac{N-1}{2}}}g\left(k\right).$$

(18)

(ii) N even. If N is even, then $N-1$ is odd. We apply Simpson’s rule on the interval $[1,N-1]$ as above, and approximate the remaining interval $[N-1,N]$ by the trapezoidal rule. Namely, we define

$$h\left(N\right)h(N-1)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{N-1}}+{\displaystyle \frac{1}{N}}\right).$$

(19)

Lemma 1 

(Representation by integer reciprocals). The Simpson-regularized harmonic sum $h\left(N\right)$ can be written as a finite linear combination of the reciprocals of positive integers.

(i) N odd. If N is odd, then

$$h\left(N\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{3}}\sum _{n=2}^{N-1}{\displaystyle \frac{c_n}{n}}+{\displaystyle \frac{1}{3N}},$$

(20)

where the coefficients $c_n$ are given by

$$c_n=\left\{\begin{array}{cc}4,\hfill & \mathrm{if}n\mathrm{is}\mathrm{even},\hfill \\ 2,\hfill & \mathrm{if}n\mathrm{is}\mathrm{odd}.\hfill \end{array}\right.$$

(21)

(ii) N even. If N is even, then

$$h\left(N\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{3}}\sum _{n=2}^{N-1}{\displaystyle \frac{c_n}{n}}+{\displaystyle \frac{1}{6(N-1)}}+{\displaystyle \frac{1}{2N}}.$$

(22)

Lemma 2. 

The constant δ defined in Definition 2 can be written as

$$\delta =\underset{n\to \infty }{lim}\{h\left(n\right)-log\left[n\right]\}$$

(23)

Proof. 

By definition, $\delta ={\sum }_{n=2}^{\infty }d\left(n\right)$ with $d\left(n\right)=g\left(n\right)-f\left(n\right)$. Summing the exact logarithmic increments defined in equation (3),

$$\sum _{n=1}^{N}f\left(n\right)=log[2N+1].$$

(24)

yields a telescoping sum, from which the claim follows. □

Theorem 1 

(Limit formula for $H\left(n\right)-h\left(n\right)$).

$$\underset{n\to \infty }{lim}\{H\left(N\right)-h\left(N\right)\}={\displaystyle \frac{log\left[2\right]+1}{3}}=0.56438\dots .$$

(25)

Proof. 

Let us compare the coefficients of $1/n$ in the representations of $h\left(N\right)$ and $H\left(N\right)$.

For even n, the coefficient in $h\left(N\right)$ is $4/3$, hence

$${\displaystyle \frac{1}{n}}-{\displaystyle \frac{4}{3n}}=-{\displaystyle \frac{1}{3n}}.$$

(26)

For interior odd $n\ge 3$, the coefficient in $h\left(N\right)$ is $2/3$, and thus

$${\displaystyle \frac{1}{n}}-{\displaystyle \frac{2}{3n}}={\displaystyle \frac{1}{3n}}.$$

(27)

Therefore, up to boundary terms, we obtain

$$H\left(N\right)-h\left(N\right)={\displaystyle \frac{1}{3}}\left(\sum _{\begin{array}{c}3\le k\le N-2\\ k\mathrm{odd}\end{array}}{\displaystyle \frac{1}{k}}-\sum _{\begin{array}{c}2\le k\le N-1\\ k\mathrm{even}\end{array}}{\displaystyle \frac{1}{k}}\right)+{\displaystyle \frac{2}{3}}+O\left({\displaystyle \frac{1}{N}}\right).$$

(28)

The alternating harmonic series

$$\sum _{k=1}^{\infty }{(-1)}^{k-1}{\displaystyle \frac{1}{k}}$$

(29)

is known to converge to $log\left[2\right]$. Since the boundary contributions vanish as $N\to \infty$, we conclude that

$$\underset{n\to \infty }{lim}\{H\left(n\right)-h\left(n\right)\}={\displaystyle \frac{log\left[2\right]+1}{3}}.$$

(30)

□

Proposition 1 

(Recovery of the Euler–Mascheroni constant). From equations (12) and (28), we obtain

$$H\left(N\right)-log\left[N\right]={\displaystyle \frac{log\left[2\right]+1}{3}}+\delta +O\left({\displaystyle \frac{1}{N}}\right).$$

(31)

In particular,

$$\gamma =\underset{N\to \infty }{lim}\{H\left(N\right)-log\left[N\right]\}={\displaystyle \frac{log\left[2\right]+1}{3}}+\delta .$$

(32)

Equivalently, this yields the decomposition

$$\gamma ={\displaystyle \frac{log\left[2\right]+1}{3}}+\delta ,$$

(33)

which expresses the Euler–Mascheroni constant, as defined in equation (2), in the form of a sum of two constants.

4. Convergence of the Simpson-Regularized Harmonic Sum

4.1. Convergence of $h\left(N\right)$

Theorem 2 

(Convergence of $h\left(N\right)$). There exists a real constant δ and a function $r\left(N\right)$ such that

$$\left|h\left(N\right)-log\left[N\right]-\delta \right|\le r\left(N\right)$$

(34)

for all integers $N\ge 3$, where $r\left(N\right)$ is defined as follows.

(i) 

If N is odd and $N\ge 3$, then

$$r\left(N\right):={\displaystyle \frac{1}{15{(N-1)}^4}}.$$

(35)

(ii) 

If N is even and $N\ge 4$, then

$$r\left(N\right):={\displaystyle \frac{1}{15{(N-2)}^4}}+{\displaystyle \frac{1}{6{(N-1)}^3}}.$$

(36)

We define the error term

$$\epsilon \left(N\right):=h\left(N\right)-log\left[N\right]-\delta .$$

(37)

Then $\epsilon \left(N\right)$ satisfies

$$\left|\epsilon \right(N\left)\right|\le r\left(N\right).$$

(38)

In particular, $h\left(N\right)-log\left[N\right]\to \delta$ as $N\to \infty$.

Proof. 

We first treat the case where N is odd.

Let $f\left(x\right)=1/x$. On each subinterval $[2j-1,2j+1]$, application of Simpson’s rule with step size 1 yields a local error of the form

$$-{\displaystyle \frac{1}{90}}{f}^{\left(4\right)}\left({\xi }_j\right),{\xi }_j\in (2j-1,2j+1).$$

(39)

See, for example, [2].

Since

$$f^{\left(4\right)}\left(x\right)=-{\displaystyle \frac{24}{x^5}},$$

(40)

and $|f^{\left(4\right)}\left(x\right)|$ is positive and strictly decreasing for $x>0$, we obtain

$$|f^{\left(4\right)}\left({\xi }_j\right)|\le {\displaystyle \frac{24}{{(2j-1)}^5}}.$$

(41)

Hence the absolute error on each subinterval is bounded by

$${\displaystyle \frac{4}{15}}{(2j-1)}^{-5}.$$

(42)

Summing the tail of these local errors, we define

$$\delta :=\underset{N\to \infty }{lim}\left(h\left(N\right)-log\left[N\right]\right),$$

(43)

which exists since the error terms are positive and summable. Moreover,

$$\sum _{j\ge (N+1)/2}{(2j-1)}^{-5}\le {\int }_{N-1}^{\infty }{x}^{-5}dx={\displaystyle \frac{1}{4{(N-1)}^4}}.$$

(44)

Multiplying by the factor $4/15$ yields

$$\left|h\left(N\right)-log\left[N\right]-\delta \right|\le {\displaystyle \frac{1}{15{(N-1)}^4}},$$

(45)

which proves (i).

We next consider even N. By definition,

$$h\left(N\right)=h(N-1)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{N-1}}+{\displaystyle \frac{1}{N}}\right).$$

(46)

The first term $h(N-1)$ is estimated by part (i). For the remaining interval $[N-1,N]$, the trapezoidal rule gives an error of the form

$$-{\displaystyle \frac{1}{12}}{f}^{\prime \prime }\left({\eta }_N\right),{\eta }_N\in (N-1,N),$$

(47)

and since

$$f^{\prime \prime }\left(x\right)={\displaystyle \frac{2}{x^3}},$$

(48)

its absolute value is bounded by $1/\left(6{(N-1)}^3\right)$. Combining these two estimates proves (ii). □

Definition 4 

(A rational sequence associated with $\delta$). Let us define a sequence $a\left(N\right)$ by

$$a\left(N\right)2h\left(N\right)-h\left(N^2\right).$$

(49)

Since each $h\left(N\right)$ is rational, each $a\left(N\right)$ is rational.

Lemma 3. 

The constant δ satisfies

$$\delta =\underset{n\to \infty }{lim}a\left(N\right)\approx 0.01283\dots .$$

(50)

Proof. 

By Theorem 2, there exists a real constant $\delta$ and a function $r\left(N\right)\ge 0$ such that

$$\left|h\left(N\right)-log\left[N\right]-\delta \right|\le r\left(N\right)$$

(51)

for all integers $N\ge 3$, where $r\left(N\right)$ is defined explicitly according to the parity of N. Hence we may write

$$h\left(N\right)=log\left[N\right]+\delta +\epsilon \left(N\right),$$

(52)

with an error term $\epsilon \left(N\right)$ satisfying

$$\left|\epsilon \right(N\left)\right|\le r\left(N\right).$$

(53)

Substituting this representation into equation (49) yields

$$a\left(N\right)=\delta +2\epsilon \left(N\right)-\epsilon \left(N^2\right).$$

(54)

Therefore,

$$|a\left(N\right)-\delta |=|2\epsilon \left(N\right)-\epsilon \left(N^2\right)|\le 2r\left(N\right)+r\left(N^2\right).$$

(55)

(i) Suppose that N is odd. Then $N^2$ is also odd, and by the definition of $r\left(N\right)$ we have

$$r\left(N\right)={\displaystyle \frac{1}{15{(N-1)}^4}},r\left(N^2\right)={\displaystyle \frac{1}{15{(N^2-1)}^4}}.$$

Hence

$$|a\left(N\right)-\delta |\le {\displaystyle \frac{2}{15{(N-1)}^4}}+{\displaystyle \frac{1}{15{(N^2-1)}^4}}\le {\displaystyle \frac{1}{5{(N-1)}^4}},$$

(56)

which proves (i).

(ii) Suppose that N is even. Then $N^2$ is odd. Using the even-N bound for $r\left(N\right)$ and the odd-N bound for $r\left(N^2\right)$, we obtain

$$|a\left(N\right)-\delta |\le {\displaystyle \frac{2}{15{(N-2)}^4}}+{\displaystyle \frac{1}{6{(N-1)}^3}}+{\displaystyle \frac{1}{15{(N^2-1)}^4}},$$

(57)

which proves (ii). □

The numerical values in Table 1 illustrate the rapid convergence of $a\left(N\right)$ toward the constant $\delta$. This behavior can be explained by the cancellation of the dominant $O(1/n)$ error terms in the asymptotic expansion of $h\left(n\right)$.

Remark 1. 

The slow convergence traditionally observed in numerical approximations to the Euler–Mascheroni constant γ [3] can be traced to the $O(1/n)$ contribution associated with the alternating harmonic series introduced in equation (29). In the present framework, this contribution is isolated explicitly, and removing it leads to the constant δ, whose numerical values exhibit markedly faster convergence (see Table 1).

5. Conclusions

This study demonstrated that the convergence of the harmonic series minus the logarithm can be established using solely Simpson’s rule and the convexity of the function $1/x$. This provides a simple, purely numerical, and analytic alternative to classical proofs based on the Euler–Maclaurin formula.

In addition, we introduced the constant $\delta$, defined as the limit of a rational sequence arising from Simpson–type harmonic approximations. Combined with the identity established in equation (33)

From a numerical perspective, this formulation explains the accelerated convergence of the associated rational sequence, which results from the systematic cancellation of leading error terms.

More broadly, the approach illustrates how elementary quadrature methods can shed new light on classical constants and suggests that similar ideas may be fruitfully applied to higher-order Newton–Cotes formulas [4].