How to caculate the summation of successive modular inverses?
The sum you are inquiring about is pretty much the right side of the Fermat quotient identity discovered by Eisenstein (proof here):
$$
-2q_p(2) = \sum_{i=1}^{(p-1)/2} 1/i \pmod{p}\,,
$$
where $q_p(a) = \frac{a^{p-1} - 1}{p}$.
Thus the sum can be computed as $-2 \frac{2^{p-1} - 1}{p} \bmod p$. Because computing $2^{p-1}$ is not feasible for very large $p$, and because one only requires the result divided by $p$ and then reduced again by $p$, computing $\frac{2^{p-1} - 1}{p}$ can be performed as $\frac{(2^{p-1} \bmod p^2) - 1}{p}$ without any precision loss. And to absorb that factor of $2$, one can compute instead $-\frac{(2^{p} \bmod p^2) - 1}{p} \bmod p$, which is exactly the formula you found.