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

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 was computed by Mathematica in just about $0.11$ sec.  


[1]: https://link.springer.com/article/10.1007/s00211-018-0978-y 
[2]: https://arxiv.org/abs/1511.03247