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