Wolstenholme's Theorem, Vandermonde's Convolution Identity and so much more

We will take $H_n = 1^{-1} + 2^{-1} + \cdots + n^{-1}$ to be the sum of the multiplicative inverses of $1,2,\dots,n$ in $\mathbb{Z}_p$ . The value of the remainder computed this way will be equal to the remainder computed with non-integer values of $H_n$ because no multiple of $p$ ever appears in the denominator of any term and $(p - 1)!$ clears all the denominators that arise from $H_n$ .

Let $s$ denote the value of the sum in question. To compute the remainder of $s \div p$ , we can ignore the terms corresponding to $n = 1,2,\dots,\frac 12(p - 1)$ since, for these terms, $\binom{2p - 2n}{p - n} \equiv \frac {(2p - 2n)(2p - 2n - 1) \cdots p \cdots (p - n + 1)}{(p - n)!} \equiv 0 \pmod p.$ A convenient way to sum the terms corresponding to $n = \frac 12(p + 1), \dots, p - 1$ is to reverse the convolution of $H_n \cdot 4^n$ and $\binom{2n}n$ . Put $m = p - n$ . The terms corresponding to $n = 1,\dots,\frac 12(p - 1)$ are the terms we can ignore, and they are exactly those terms corresponding to $m = \frac 12(p + 1), \dots, p - 1$ . Omitting these terms shows $s \equiv (p - 1)! \sum_{m = 1}^{\frac 12(p - 1)} H_{p - m} \cdot 4^{p - m} \cdot \binom{2m}m \pmod p.$ By Wilson's theorem, $(p - 1)! \equiv -1 \pmod p,$ and by Fermat's little theorem, $4^{p - m} \equiv 4 \cdot 4^{-m} \pmod p,$ so we have $s \equiv -4 \sum_{m = 1}^{\frac 12(p - 1)} H_{p - m} \cdot 4^{-m} \cdot \binom{2m}m \pmod p. \tag{1}$

By Wolstenholme's theorem, $H_{p - 1} \equiv 0 \pmod p,$ and for $1 \leq m \leq p - 1$ , $H_{p - m - 1} \equiv H_{p - m} - (p - m)^{-1} \pmod p.$ Also, $\binom 21 \equiv 2 \pmod p,$ and one can show that $\binom {2m + 2}{m + 1} \equiv (2m + 2)(2m + 1)(m + 1)^{-2} \binom {2m}m \pmod p.$ Using the previous four congruences, we can compute $H_{p - m}$ and $\binom {2m}m$ for each term on the right side of $(1)$ recursively. To find multiplicative inverses, we can apply Fermat's little theorem again, which states $a^{-1} \equiv a^{p - 2} \pmod p, \qquad a = 1,2,\dots,p - 1.$ The script below computes the summands one by one and takes the residues of the partial sums. Using this script, we obtain the result $s \equiv -1 \sum_{m = 1}^{\frac 12(p - 1)} H_{p - m} 4^{-m} \binom {2m}m \equiv \boxed{32761} \pmod p.$

p = 2**16 + 1
inverses = [None]       #store multiplicative inverses
for k in range(1,p):    #use Fermat's little theorem to compute inverses
    inverses.append(pow(k,p - 2,p))
sum = 0                 #store partial sums
H = 0                   #initial value of H_m
expon = inverses[4]     #initial value of 4^-m 
binom = 2               #initial value of binom(2n,n)
for m in range(1,int((p + 1)/2)):
    term = -4 * (H * expon * binom) % p
    sum = (sum + term) % p
    H = (H - inverses[p - m]) % p
    expon = (expon * inverses[4]) % p
    binom = (binom * 2*(m + 1) * (2*m + 1) *inverses[m + 1]**2) % p
print(sum)