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][1] for details or [Section 6.1][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][3]. 

---  

Using the [summation formula][1] 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).  


[1]: https://link.springer.com/article/10.1007/s00211-018-0978-y 
[2]: https://arxiv.org/abs/1511.03247 
[3]: https://www.wolframcloud.com/obj/ipinelis/Published/1.nb