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.
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).
