Let $H_n = 1 + \half + \cdots + {1\over n}$ be the $n$-th partial sum
of the harmonic series. A classical result of Wolstenholme states
that, if $p > 3$ is prime, the numerator of $H_{p-1}$ is divisible by
$p^2$. Here we consider, for a given prime $p$, the set $J_p$ of $n$
for which $p$ divides the numerator of $H_n$. This set $J_p$ had been
previously determined for $p = 2,3,5,7$. One of our results is that
$J_{11}$ contains exactly $638$ integers, the largest of which is a
number of $31$ decimal digits. We determine $J_p$ for all $p < 550$
with three exceptions: $83$, $127$ and $397$.
¶ The computation is based on a new $p$-adically convergent formula for
the quantity $H_{pn} - H_n/p$. We describe a probabilistic model
for the sets $J_p$, based on branching processes. The model predicts
that $|J_p| = O(p^2 (\log\log p)^{2+\eps} )$, and that there are
infinitely many $p$ with $|J_p| \ge p^2 (\log\log p)^2$. This
strengthens an earlier conjecture of Eswarathasan and Levine that
$|J_p|$ is finite for all $p$. Another prediction of the model is
that there will be infinitely many pairs $(n,p)$ for which $p^3$
divides the numerator of $H_n$, but only finitely many for which $p^4$
divides $H_n$.
¶ It has been conjectured that there are infinitely many $p$ for which
$|J_p| = 3$. We give a probabilistic argument that suggests that such
primes have a density $1/e$ in the set of all primes, and
experimentally confirm this by a determination of all such $p \le
10^5$.