I need a good approximations for $H\_p$, for $p \in (0,1) \cap \mathbb{Q}$, the generalization of $H\_n=\sum\_{i=1}^n \frac{1}{i}$ to the real numbers.


I tried $H\_p = p \sum\_{k=1}^\infty \frac{1}{k (k + p)}$ and $H\_p = \sum\_{k=0}^\infty \sum\_{j=0}^\infty (-1)^j (1 + k)^{-2 - j} p^{1 + j}$ for $-\frac{1}{2}<p<\frac{1}{2}$, both converge very slow, just five digit of $H\_\frac{1}{2}$ take ~1,000,000 iterations.



Approximation for interpolation of harmonic numbers
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](https://link.springer.com/article/10.1007/s00211-018-0978-y) for details or [Section 6.1](https://link.springer.com/article/10.1007/s00211-018-0978-y) 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](https://www.wolframcloud.com/obj/ipinelis/Published/1.nb).




---


Using the [summation formula](https://link.springer.com/article/10.1007/s00211-018-0978-y) 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).


