An Algorithm for Calculating Terms of the Stirling’s Formula Remainder

1. Introduction

In this note, we study Stirling’s approximation to $N!$, a fundamental result in asymptotic analysis and number theory. Originally derived in the 18th century, the approximation provides an accurate estimate of the factorial function for large values of N. Over the centuries, various forms and refinements of this approximation have been developed. A particularly elegant and instructive derivation, based solely on elementary transformations of integrals and term-by-term integration of a series, was introduced by Marsaglia [1].

Our first objective is to briefly derive Stirling’s formula. We begin with the representation of $lnN!$ as a sum of logarithms:

$$\begin{array}{cc}\hfill lnN!=& ln1+ln2+...+lnN=\hfill \\ \hfill =& {\int }_1^{N}lnxdx+\sum _{k=1}^{N-1}{\displaystyle \frac{ln(k+1)-lnk}{2}}-\sum _{k=1}^{N-1}{D}_k,\hfill \end{array}$$

(1)

where each term $D_k$ measures the discrepancy between the trapezoidal rule and the exact integral of $lnx$ over the interval $[k,k+1]$. That is,

$$\begin{array}{cc}\hfill D_k=& {\int }_k^{k+1}lnxdx-{\displaystyle \frac{lnk+ln(k+1)}{2}}=\hfill \\ \hfill =& \left(xlnx-x\right){|}_k^{k+1}-{\displaystyle \frac{lnk+ln(k+1)}{2}}=\hfill \\ \hfill =& (k+1)ln(k+1)-klnk-1-{\displaystyle \frac{lnk+ln(k+1)}{2}}=\hfill \\ \hfill =& \left(k+{\displaystyle \frac{1}{2}}\right)ln\left(1+{\displaystyle \frac{1}{k}}\right)-1.\hfill \end{array}$$

(2)

The integral in (1) can be evaluated explicitly,

$$\begin{array}{cc}\hfill {\int }_1^{N}lnxdx=& NlnN-N+1,\mathrm{and}\sum _{k=1}^{N-1}{\displaystyle \frac{ln(k+1)-lnk}{2}}={\displaystyle \frac{lnN}{2}}\hfill \end{array}$$

is a telescoping sum, as the terms cancel pairwise. In this way,

$$\begin{array}{cc}\hfill lnN!=& NlnN-N+1+{\displaystyle \frac{lnN}{2}}-\sum _{k=1}^{N-1}{D}_k=\hfill \\ \hfill =& ln\left(N^N\right)-N+1+ln\left(N^{1/2}\right)-\sum _{k=1}^{N-1}{D}_k,\hfill \end{array}$$

which leads to

$$\begin{array}{c}\hfill N!=N^{N+1/2}{\mathrm{e}}^{1-N}exp\left(-\sum _{k=1}^{N-1}{D}_k\right),\end{array}$$

(3)

where $D_k$ is given in (2).

It is well known, from the Mercator series, that

$$\begin{array}{cc}\hfill ln(1+x)=& x-{\displaystyle \frac{x^2}{2}}+{\displaystyle \frac{x^3}{3}}-{\displaystyle \frac{x^4}{4}}+...=\sum _{i=1}^{+\infty }{(-1)}^{i+1}{\displaystyle \frac{x^i}{i}},-1<x\le 1.\hfill \end{array}$$

Substituting $x=1/k$, where $k=1,2,3,...$, yields

$$\begin{array}{cc}\hfill ln\left(1+{\displaystyle \frac{1}{k}}\right)=& {\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{2k^2}}+{\displaystyle \frac{1}{3k^3}}-{\displaystyle \frac{1}{4k^4}}+...=\sum _{i=1}^{+\infty }{(-1)}^{i+1}{\displaystyle \frac{1}{i{k}^i}}.\hfill \end{array}$$

Thus, $D_k$ in (2) can be rewritten as

$$\begin{array}{cc}\hfill D_k=& \left(k+{\displaystyle \frac{1}{2}}\right)ln\left(1+{\displaystyle \frac{1}{k}}\right)-1=\hfill \\ \hfill =& \left(1-{\displaystyle \frac{1}{2k}}+{\displaystyle \frac{1}{3k^2}}-{\displaystyle \frac{1}{4k^3}}+...\right)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{2k^2}}+{\displaystyle \frac{1}{3k^3}}-...\right)-1=\hfill \\ \hfill =& \left({\displaystyle \frac{1}{3}}-{\displaystyle \frac{1}{2.2}}\right){\displaystyle \frac{1}{k^2}}-\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{1}{2.3}}\right){\displaystyle \frac{1}{k^3}}+\left({\displaystyle \frac{1}{5}}-{\displaystyle \frac{1}{2.4}}\right){\displaystyle \frac{1}{k^4}}-...=\hfill \\ \hfill =& \sum _{i=2}^{+\infty }{(-1)}^i\left({\displaystyle \frac{1}{i+1}}-{\displaystyle \frac{1}{2i}}\right){\displaystyle \frac{1}{k^i}}=\sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{a_i}{k^i}},\hfill \end{array}$$

where $a_i=1/(i+1)-1/\left(2i\right)$, for $i=2,3,4,...$, is a positive decreasing sequence of real numbers that converges to zero.

Thus,

$$\begin{array}{cc}\hfill \sum _{k=1}^{N-1}{D}_k=& \sum _{k=1}^{N-1}\sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{a_i}{k^i}}=\sum _{i=2}^{+\infty }{(-1)}^i{a}_i\left(\sum _{k=1}^{N-1}{\displaystyle \frac{1}{k^i}}\right).\hfill \end{array}$$

From the Riemann Zeta function (see [2] for details),

$$\begin{array}{cc}\hfill \zeta \left(x\right)=& \sum _{j=1}^{+\infty }{\displaystyle \frac{1}{j^x}},\hfill \end{array}$$

which is convergent for $x\in \mathbf{R}$ with $x>1$, it follows that

$$\begin{array}{cc}\hfill \sum _{k=1}^{N-1}{\displaystyle \frac{1}{k^i}}=& \zeta \left(i\right)-\sum _{k=N}^{+\infty }{\displaystyle \frac{1}{k^i}},\hfill \end{array}$$

and, therefore,

$$\begin{array}{cc}\hfill \sum _{k=1}^{N-1}{D}_k=& \sum _{i=2}^{+\infty }{(-1)}^i{a}_i\zeta \left(i\right)-\sum _{i=2}^{+\infty }{(-1)}^i{a}_i\left(\sum _{k=N}^{+\infty }{\displaystyle \frac{1}{k^i}}\right).\hfill \end{array}$$

From [2,3],

$$\begin{array}{cc}\hfill \sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{\zeta \left(i\right)}{2i}}=& {\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\zeta \left(2\right)}{2}}-{\displaystyle \frac{\zeta \left(3\right)}{3}}+{\displaystyle \frac{\zeta \left(4\right)}{4}}-{\displaystyle \frac{\zeta \left(5\right)}{5}}+...\right)={\displaystyle \frac{\gamma }{2}},\hfill \end{array}$$

where $\gamma$ is the Euler-Mascheroni constant.

Additionally, from [3,4], based on the Riemann Zeta function and Rational Zeta series, it follows that

$$\begin{array}{cc}\hfill \sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{\zeta \left(i\right)-1}{i+1}}=& {\displaystyle \frac{1}{2}}\left(\gamma +3-ln\sqrt{2\pi }\right)-ln2.\hfill \end{array}$$

Moreover, using the Mercator series evaluated at $x=1$,

$$\begin{array}{cc}\hfill \sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{1}{i+1}}=& ln2-{\displaystyle \frac{1}{2}},\hfill \end{array}$$

we arrive at the Suryanarayana formula (see [3] for example),

$$\begin{array}{cc}\hfill \sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{\zeta \left(i\right)}{i+1}}=& {\displaystyle \frac{1}{2}}\left(\gamma +3-ln\sqrt{2\pi }\right)-ln2+\left(ln2-{\displaystyle \frac{1}{2}}\right)={\displaystyle \frac{\gamma }{2}}+1-ln\sqrt{2\pi },\hfill \end{array}$$

and, consequently,

$$\begin{array}{cc}\hfill \sum _{i=2}^{+\infty }{(-1)}^i{a}_i\zeta \left(i\right)=& {\displaystyle \frac{\gamma }{2}}+1-ln\sqrt{2\pi }-{\displaystyle \frac{\gamma }{2}}=1-ln\sqrt{2\pi }.\hfill \end{array}$$

In this way,

$$\begin{array}{cc}\hfill \sum _{k=1}^{N-1}{D}_k=& 1-ln\sqrt{2\pi }-\sum _{i=2}^{+\infty }\sum _{k=N}^{+\infty }{(-1)}^i{\displaystyle \frac{a_i}{k^i}}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill exp\left(-\sum _{k=1}^{N-1}{D}_k\right)=& {\mathrm{e}}^{-1}\sqrt{2\pi }{\mathrm{e}}^{r_N},\hfill \end{array}$$

where $r_N$ is given by

$$\begin{array}{cc}\hfill r_N=& \sum _{i=2}^{+\infty }\sum _{k=N}^{+\infty }{(-1)}^i{\displaystyle \frac{a_i}{k^i}},\hfill \end{array}$$

(4)

and, as before, $a_i=1/(i+1)-1/\left(2i\right)$, for $i=2,3,4,...$, is a positive decreasing sequence of real numbers that converges to zero.

Finally, Eq. (3), known in the literature as Stirling’s formula, can be expressed as

$$\begin{array}{cc}\hfill N!=& \sqrt{2\pi }{N}^{N+1/2}{\mathrm{e}}^{-N}{\mathrm{e}}^{r_N},\hfill \end{array}$$

(5)

where the remainder $r_N$ is given in (4).

Eq (5) was first presented in [5], where it was also established that $11/2<r_N+ln\sqrt{2\pi }<1$. Since then, the term $r_N$ has been the subject of detailed study (see [6,7,8,9] for historical notes), particularly through works such as [10], who first provided sharp two-sided bounds for it,

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{12N+1}}<r_N<{\displaystyle \frac{1}{12N}}.\end{array}$$

Several years later, [11] improved the estimation of $r_N$,

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}<r_N<{\displaystyle \frac{1}{12N}}.\end{array}$$

More recently, [12] obtained the following result,

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}<r_N<{\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N(N+1)(N+2)}},\end{array}$$

[8] proved that

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{12N+c_1}}\le r_N<{\displaystyle \frac{1}{12N+c_2}},\end{array}$$

with the best possible constants $c_1=-12+1/(1-ln\sqrt{2\pi })\approx 0.336317$ and $c_2=0$, and [9] derived the following inequality:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{12N}}-{\displaystyle \frac{7}{360N(7N^2+2)}}-{\displaystyle \frac{13}{35280N^7}}<r_N<{\displaystyle \frac{1}{12N}}-{\displaystyle \frac{7}{360N(7N^2+2)}}.\end{array}$$

It is also important to highlight the works [1,7,13,14,15] who focused on obtaining explicit formulas for the coefficients of the asymptotic expansion of $r_N$.

The main goal of this work is to present an algorithm, relying only on elementary mathematical methods techniques, for efficiently calculating many terms of $r_N$ exactly. To the best of our knowledge, we understand that our algorithm represents an alternative to the procedures presented herein.

2. The Algorithm

Observe that $r_N$ in (4) can be rewritten as

$$\begin{array}{cc}\hfill r_N=& \sum _{k=N}^{+\infty }\sum _{i=2}^{+\infty }{(-1)}^i{\displaystyle \frac{a_i}{k^i}}=\sum _{k=N}^{+\infty }\left({\displaystyle \frac{a_2}{k^2}}-{\displaystyle \frac{a_3}{k^3}}+{\displaystyle \frac{a_4}{k^4}}-{\displaystyle \frac{a_5}{k^5}}+\dots \right),\hfill \end{array}$$

a numerical series whose general term forms a convergent alternating series for each N. Define

$$\begin{array}{cc}\hfill S_j\left(k\right)=& \sum _{i=2}^j{(-1)}^i{\displaystyle \frac{a_i}{k^i}},\hfill \end{array}$$

the sequence of partial sums of this series.

Step 1: Consider the polynomial $p_1$ defined by the following equation:

$$\begin{array}{cc}\hfill & S_4\left(k\right)-{\alpha }_1\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)={\displaystyle \frac{p_1\left(k\right)}{k^4(k+1)}},\hfill \end{array}$$

where ${\alpha }_1$ is a real parameter to be determined later. This leads to:

$$\begin{array}{cc}\hfill p_1\left(k\right)=& k^4(k+1)\left(S_4\left(k\right)-{\alpha }_1\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)\right)=\hfill \\ \hfill =& k^4(k+1)\left({\displaystyle \frac{a_2}{k^2}}-{\displaystyle \frac{a_3}{k^3}}+{\displaystyle \frac{a_4}{k^4}}-{\alpha }_1\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)\right)=\hfill \\ \hfill =& (a_2-{\alpha }_1)k^3+(a_2-a_3)k^2+(a_4-a_3)k+a_4.\hfill \end{array}$$

Now, ${\alpha }_1$ is chosen in order to reduce the degree of $p_1$ to the lowest possible value. In this case, ${\alpha }_1=a_2=1/12$ and, after some calculations, we obtain

$$\begin{array}{cc}\hfill p_1\left(k\right)=& (a_4-a_3)k+a_4={\displaystyle \frac{-1}{120}}k+{\displaystyle \frac{3}{40}}\hfill \end{array}$$

which is a polynomial of degree 1, that is, $deg\left(p_1\right)=1$. We denote by $\eta \left(p_1\right)=-1/120$ the coefficient associated with the highest degree term of the polynomial $p_1$.

In this way, the partial sum $S_4$ can be expressed as

$$\begin{array}{cc}\hfill S_4\left(k\right)=& {\displaystyle \frac{a_2}{k^2}}-{\displaystyle \frac{a_3}{k^3}}+{\displaystyle \frac{a_4}{k^4}}={\displaystyle \frac{1}{12}}\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)+{\displaystyle \frac{p_1\left(k\right)}{k^4(k+1)}}.\hfill \end{array}$$

As an important remark, $r_N$ in (4) can then be approximated by

$$\begin{array}{cc}\hfill r_N\approx & \sum _{k=N}^{+\infty }{S}_4\left(k\right)={\displaystyle \frac{1}{12}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_1\left(k\right)}{k^4(k+1)}}=\hfill \\ \hfill =& {\displaystyle \frac{1}{12N}}+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_1\left(k\right)}{k^4(k+1)}}.\hfill \end{array}$$

The sum above can be estimated using the integral test, comparing it to an improper integral.

Step 2: Similarly, we seek a polynomial $p_3$ that satisfies

$$\begin{array}{cc}\hfill & {\displaystyle \frac{p_1\left(k\right)}{k^4(k+1)}}-{\displaystyle \frac{a_5}{k^5}}+{\displaystyle \frac{a_6}{k^6}}-{\alpha }_3\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)={\displaystyle \frac{p_3\left(k\right)}{k^6{(k+1)}^3}},\hfill \end{array}$$

where ${\alpha }_3$ is a real parameter to be determined.

Thus,

$$\begin{array}{cc}\hfill p_3\left(k\right)=& k^2{(k+1)}^2{p}_1\left(k\right)-a_{5}k{(k+1)}^3+a_6{(k+1)}^3-{\alpha }_3{k}^3\left({(k+1)}^3-k^3\right).\hfill \end{array}$$

Again, ${\alpha }_3$ is chosen with the aim of reducing the degree of $p_3$ to the lowest possible value. In this case, ${\alpha }_3=\eta \left(p_1\right)/3=-1/360$.

Therefore, we can present $S_6\left(k\right)$ as

$$\begin{array}{cc}\hfill S_6\left(k\right)=& {\displaystyle \frac{1}{12}}\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)-{\displaystyle \frac{1}{360}}\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)+{\displaystyle \frac{p_3\left(k\right)}{k^6{(k+1)}^3}},\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill p_3\left(k\right)=& {\displaystyle \frac{1}{252}}{k}^3+{\displaystyle \frac{3}{56}}{k}^2+{\displaystyle \frac{47}{420}}k+{\displaystyle \frac{5}{84}}.\hfill \end{array}$$

Note that $deg\left(p_3\right)=3$ and $\eta \left(p_3\right)=1/252$.

Again, as a remark, $r_N$ in (4) is better approximated by

$$\begin{array}{cc}\hfill r_N\approx & \sum _{k=N}^{+\infty }{S}_6\left(k\right)={\displaystyle \frac{1}{12}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)-{\displaystyle \frac{1}{360}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)+\hfill \\ \hfill +& \sum _{k=N}^{+\infty }{\displaystyle \frac{p_3\left(k\right)}{k^6{(k+1)}^3}}={\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_3\left(k\right)}{k^6{(k+1)}^3}}.\hfill \end{array}$$

Step 3: Similarly, aiming for a better approximation of $r_N$, we write

$$\begin{array}{cc}\hfill & {\displaystyle \frac{p_3\left(k\right)}{k^6{(k+1)}^3}}-{\displaystyle \frac{a_7}{k^7}}+{\displaystyle \frac{a_8}{k^8}}-{\alpha }_5\left({\displaystyle \frac{1}{k^5}}-{\displaystyle \frac{1}{{(k+1)}^5}}\right)={\displaystyle \frac{p_5\left(k\right)}{k^8{(k+1)}^5}},\hfill \end{array}$$

where ${\alpha }_5$ is, as before, a real parameter to be determined.

The polynomial $p_5$ can be rewritten as

$$\begin{array}{cc}\hfill p_5\left(k\right)=& k^2{(k+1)}^2{p}_3\left(k\right)-a_{7}k{(k+1)}^5+a_8{(k+1)}^5-{\alpha }_5{k}^3\left({(k+1)}^5-k^5\right).\hfill \end{array}$$

The parameter ${\alpha }_5$ is chosen with the aim of reducing the degree of $p_5$ to the lowest possible value. Thus, ${\alpha }_5=\eta \left(p_3\right)/5=1/1260$.

Then, $S_8\left(k\right)$ is given by

$$\begin{array}{cc}\hfill S_8\left(k\right)=& {\displaystyle \frac{1}{12}}\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)-{\displaystyle \frac{1}{360}}\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)+\hfill \\ \hfill +& {\displaystyle \frac{1}{1260}}\left({\displaystyle \frac{1}{k^5}}-{\displaystyle \frac{1}{{(k+1)}^5}}\right)+{\displaystyle \frac{p_5\left(k\right)}{k^8{(k+1)}^5}},\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill p_5\left(k\right)=& -{\displaystyle \frac{1}{240}}{k}^5+{\displaystyle \frac{29}{720}}{k}^4+{\displaystyle \frac{13}{72}}{k}^3+{\displaystyle \frac{5}{18}}{k}^2+{\displaystyle \frac{191}{1008}}k+{\displaystyle \frac{7}{144}}.\hfill \end{array}$$

Note that $deg\left(p_5\right)=5$ and $\eta \left(p_5\right)=-1/240$.

Therefore, a better approximation for $r_N$ in (4) is obtained:

$$\begin{array}{cc}\hfill r_N\approx & \sum _{k=N}^{+\infty }{S}_8\left(k\right)={\displaystyle \frac{1}{12}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)-{\displaystyle \frac{1}{360}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)+\hfill \\ \hfill +& {\displaystyle \frac{1}{1260}}\sum _{k=N}^{+\infty }\left({\displaystyle \frac{1}{k^5}}-{\displaystyle \frac{1}{{(k+1)}^5}}\right)+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_5\left(k\right)}{k^8{(k+1)}^5}}=\hfill \\ \hfill =& {\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}+{\displaystyle \frac{1}{1260N^5}}+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_5\left(k\right)}{k^8{(k+1)}^5}}.\hfill \end{array}$$

Step 4: Likewise, the polynomial $p_7$ is defined as

$$\begin{array}{cc}\hfill p_7\left(k\right)=& k^2{(k+1)}^2{p}_5\left(k\right)-a_{9}k{(k+1)}^7+a_{10}{(k+1)}^7-{\alpha }_7{k}^3\left({(k+1)}^7-k^7\right).\hfill \end{array}$$

Taking ${\alpha }_7=\eta \left(p_5\right)/7=-1/1680$, we obtain

$$\begin{array}{cc}\hfill p_7\left(k\right)=& {\displaystyle \frac{1}{132}}{k}^7+{\displaystyle \frac{7}{132}}{k}^6+{\displaystyle \frac{3349}{13860}}{k}^5+{\displaystyle \frac{8119}{13860}}{k}^4+{\displaystyle \frac{10891}{13860}}{k}^3+\hfill \\ \hfill +& {\displaystyle \frac{105}{176}}{k}^2+{\displaystyle \frac{479}{1980}}k+{\displaystyle \frac{9}{220}},\hfill \end{array}$$

where $deg\left(p_7\right)=7$ and $\eta \left(p_7\right)=1/132$.

Thus,

$$\begin{array}{cc}\hfill S_{10}\left(k\right)=& {\displaystyle \frac{1}{12}}\left({\displaystyle \frac{1}{k}}-{\displaystyle \frac{1}{k+1}}\right)-{\displaystyle \frac{1}{360}}\left({\displaystyle \frac{1}{k^3}}-{\displaystyle \frac{1}{{(k+1)}^3}}\right)+\hfill \\ \hfill +& {\displaystyle \frac{1}{1260}}\left({\displaystyle \frac{1}{k^5}}-{\displaystyle \frac{1}{{(k+1)}^5}}\right)-{\displaystyle \frac{1}{1680}}\left({\displaystyle \frac{1}{k^7}}-{\displaystyle \frac{1}{{(k+1)}^7}}\right)+{\displaystyle \frac{p_7\left(k\right)}{k^{10}{(k+1)}^7}},\hfill \end{array}$$

and then

$$\begin{array}{cc}\hfill r_N\approx & \sum _{k=N}^{+\infty }{S}_{10}\left(k\right)={\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}+{\displaystyle \frac{1}{1260N^5}}-{\displaystyle \frac{1}{1680N^7}}+\sum _{k=N}^{+\infty }{\displaystyle \frac{p_7\left(k\right)}{k^{10}{(k+1)}^7}},\hfill \end{array}$$

where $deg\left(p_7\right)=7$.

In summary, the following algorithm is presented.

Step 5: Algorithm

$$\begin{array}{cc}\hfill & {\alpha }_1=1/12\hfill \\ \hfill & p_1\left(k\right)=-(1/120)k+3/40=\eta \left(p_1\right)k+3/40\hfill \\ \hfill & \mathrm{For}l=1,2,3,\dots \hfill \\ \hfill & \begin{array}{cc}\hfill & p_{2l+1}\left(k\right)=k^2{(k+1)}^2{p}_{2l-1}\left(k\right)-a_{2l+3}k{(k+1)}^{2l+1}+a_{2l+4}{(k+1)}^{2l+1}\hfill \\ \hfill & -{\alpha }_{2l+1}{k}^3\left({(k+1)}^{2l+1}-k^{2l+1}\right)\hfill \\ \hfill & {\alpha }_{2l+1}=\eta \left(p_{2l-1}\right)/(2l+1)\hfill \end{array}\hfill \end{array}$$

Remark 1.

The sequence $\left({\alpha }_{2l+1}\right)$, for $l=0,1,2,3,...$, defined above, is crucial because it leads to

$$\begin{array}{c}\hfill r_N=\sum _{l=0}^{+\infty }{(-1)}^l{\displaystyle \frac{{\alpha }_{2l+1}}{N^{2l+1}}}.\end{array}$$

Furthermore, based on computational tests, it was observed that, for $l=0,1,2,...$,

$$\begin{array}{c}\hfill deg\left(p_{2l+1}\right)=2l+1.\end{array}$$

The exact calculations of $p_{2l+1}\left(k\right)$ and ${\alpha }_{2l+1}$, for $l=1,2,3,...$, were carried out using symbolic computations in the Mathematica software [16], using the script that follows:

To illustrate our algorithm, we calculate ${\alpha }_1,{\alpha }_3,{\alpha }_5,...,{\alpha }_{39}$ using the Mathematica script provided above (CPU time approximately $0.1$ seconds, executed in the trial version of Mathematica Online) and then, $r_N$ in (4) can be approximated by

$$\begin{array}{cc}\hfill r_N=& \sum _{l=0}^{+\infty }{(-1)}^l{\displaystyle \frac{{\alpha }_{2l+1}}{N^{2l+1}}}\approx {\displaystyle \frac{1}{12N}}-{\displaystyle \frac{1}{360N^3}}+{\displaystyle \frac{1}{1260N^5}}-{\displaystyle \frac{1}{1680N^7}}+{\displaystyle \frac{1}{1188N^9}}-\hfill \\ \hfill -& {\displaystyle \frac{691}{360360N^{11}}}+{\displaystyle \frac{1}{156N^{13}}}-{\displaystyle \frac{3617}{122400N^{15}}}+{\displaystyle \frac{43867}{244188N^{17}}}-{\displaystyle \frac{174611}{125400N^{19}}}+\hfill \\ \hfill +& {\displaystyle \frac{77683}{5796N^{21}}}-{\displaystyle \frac{236364091}{1506960N^{23}}}+{\displaystyle \frac{657931}{300N^{25}}}-{\displaystyle \frac{3392780147}{93960N^{27}}}+{\displaystyle \frac{1723168255201}{2492028N^{29}}}-\hfill \\ \hfill -& {\displaystyle \frac{7709321041217}{505920N^{31}}}+{\displaystyle \frac{151628697551}{396N^{33}}}-{\displaystyle \frac{26315271553053477373}{2418179400N^{35}}}+\hfill \\ \hfill +& {\displaystyle \frac{154210205991661}{444N^{37}}}-{\displaystyle \frac{261082718496449122051}{21106800N^{39}}}.\hfill \end{array}$$

As a test, the CPU time to calculate the first 300 terms of $r_N$ was approximately 200 seconds. Remark 1 have been verified in this case.