Return to Answer

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast and very accurately as follows. Let \begin{equation*} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation*} so that \begin{equation*} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation*} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation*} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation*} where \begin{equation*} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation*} \begin{equation*} F(x):=\ln\frac{x+1}{x+1+p} \end{equation*} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation*} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation*} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation*} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation*} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation*} H_{2/3;m,c}:=\sum_{k=0}^c f(k)-G_{m,c} \end{equation*} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation, $H_{2/3;m,c}$, of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec. Since $H_{2/3}=\frac{3}{2}+\frac{\pi }{2 \sqrt{3}}-\frac{3 \ln3}{2}$, we find that the actual difference $H_{2/3;m,c}-H_{2/3}$ is $-3.08\ldots\times10^{-1003}$.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time if the desired number $d$ of correct digits is large enough (such as $2000$ or more).

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast and very accurately as follows. Let \begin{equation*} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation*} so that \begin{equation*} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation*} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation*} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation*} where \begin{equation*} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation*} \begin{equation*} F(x):=\ln\frac{x+1}{x+1+p} \end{equation*} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation*} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) \\ =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation*} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation*} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation*} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation*} H_{2/3;m,c}:=\sum_{k=0}^c f(k)-G_{m,c} \end{equation*} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation, $H_{2/3;m,c}$, of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec. Since $H_{2/3}=\frac{3}{2}+\frac{\pi }{2 \sqrt{3}}-\frac{3 \ln3}{2}$, we find that the actual difference $H_{2/3;m,c}-H_{2/3}$ between the approximate and true values is $\approx-3.08\times10^{-1003}$. The corresponding Mathematica notebook can be viewed here.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time if the desired number $d$ of correct digits is large enough (such as $2000$ or more).

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast and very accurately as follows. Let \begin{equation} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation} so that \begin{equation} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation} where \begin{equation} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation} \begin{equation} F(x):=\ln\frac{x+1}{x+1+p} \end{equation} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation} \sum_{k=0}^c f(k)-G_{m,c} \\ (\approx 0.758981249114928906960[\text{...967 digits omitted...}]071910261571) \end{equation} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time if the desired number $d$ of correct digits is large enough (such as $2000$ or more).

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast and very accurately as follows. Let \begin{equation*} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation*} so that \begin{equation*} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation*} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation*} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation*} where \begin{equation*} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation*} \begin{equation*} F(x):=\ln\frac{x+1}{x+1+p} \end{equation*} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation*} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation*} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation*} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation*} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation*} H_{2/3;m,c}:=\sum_{k=0}^c f(k)-G_{m,c} \end{equation*} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation, $H_{2/3;m,c}$, of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec. Since $H_{2/3}=\frac{3}{2}+\frac{\pi }{2 \sqrt{3}}-\frac{3 \ln3}{2}$, we find that the actual difference $H_{2/3;m,c}-H_{2/3}$ is $-3.08\ldots\times10^{-1003}$.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time if the desired number $d$ of correct digits is large enough (such as $2000$ or more).

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast as follows. Let \begin{equation} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation} so that \begin{equation} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation} where \begin{equation} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation} \begin{equation} F(x):=\ln\frac{x+1}{x+1+p} \end{equation} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation} \sum_{k=0}^c f(k)-G_{m,c} \\ (\approx 0.758981249114928906960[\text{...967 digits omitted...}]071910261571) \end{equation} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time (and also computer memory usage) if the desired number $d$ of correct digits is very large (such as $64000$).

We can compute $H_p = p \sum_{k=1}^\infty \frac{1}{k (k + p)}$ very fast and very accurately as follows. Let \begin{equation} f(x):=\frac p{(x+1)(x+1 + p)}=\frac{1}{x+1}-\frac{1}{x+1+p}, \end{equation} so that \begin{equation} H_p=\sum_{k=0}^c f(k)+\sum_{k=1}^\infty f_c(k), \end{equation} where $c$ is a natural number and $f_c(x):=f(c+x)$.

If $c$ is large enough, then $\sum_{k=1}^\infty f_c(k)$ can be evaluated fast with high precision by the Euler--Maclaurin (EM) formula: \begin{equation} \sum_{k=1}^\infty f_c(k)=-G_{m,c}+R_{m,c}, \end{equation} where \begin{equation} G_{m,c}:= F(c) + \frac{f(c)}2 + \sum_{j=1}^{m-1}\frac{B_{2j}}{(2j)!}f^{(2j-1)}(c), \end{equation} \begin{equation} F(x):=\ln\frac{x+1}{x+1+p} \end{equation} (so that $F$ is the antiderivative of $f$ with $F(\infty-)=0$), $m$ is a natural number, the $B_k$'s are the Bernoulli numbers, and \begin{equation} |R_{m,c}|<\frac{2.02}{(2\pi)^{2m-1}}f^{(2m-2)}(c) =\frac{2.02(2m-2)!}{(2\pi)^{2m-1}}((c+1)^{1-2m}-((c+1+p)^{1-2m}) \end{equation} if $m\ge4$.

If we want to get $H(p)$ with $d$ correct decimal digits after the decimal point, then good choices for $m$ and $c$ are as follows: \begin{equation} m=\Big\lceil\frac d{2\log_{10}d}\Big\rceil,\quad c=\Big\lceil\frac1{\pi e}\,m\,10^{d/(2m)}\Big\rceil; \end{equation} see Section 6.1 for details or Section 6.1 for more details.

E.g., if we want to get $H(2/3)$ with $1000$ correct decimal digits after the decimal point, then we can choose $m=167$ and $c=19288$, so that the error $|R_{m,c}|$ of the corresponding approximation \begin{equation} \sum_{k=0}^c f(k)-G_{m,c} \\ (\approx 0.758981249114928906960[\text{...967 digits omitted...}]071910261571) \end{equation} of $H(2/3)$ is $<3\times10^{-1002}$. This approximation of $H(2/3)$ was computed by Mathematica in just about $0.11$ sec.

Using the summation formula alternative to the EM formula, we can parallelize calculations and thus further reduce the execution time if the desired number $d$ of correct digits is large enough (such as $2000$ or more).