$p$-adic sums of $p$ terms My question is inspired by this riddle: Let $p \geq 5$ be prime, and let
$$ 1 + \frac 1 2 + \frac 1 3 + \dots + \frac 1 {p-1} = \frac a b $$
where $a/b$ is the fraction expressed in lowest terms. Show that $p^2$ divides $a$.
Another way to express this claim is that we are giving a $p$-adic approximation to the sum
$$ \sum_{k=0}^{p-1} \frac {1} {k+1} = \sum_{k=1}^p \frac {1} {k} = \frac {1} {p} + E $$
where $|E|_p \leq p^{-2}$. This suggests a more general question: Let $f$ be a rational function with coefficients in $\mathbb {Q}_p$. Is there an algorithm faster than direct summation for estimating (being close in the $p$-adic metric) the sum $\sum _{k=0} ^{p-1} f (k)$? Say, something that's polynomial in $\log p$? In particular, is there a way of estimating the harmonic sum as in this riddle with error $\sim p^3$?
 A: This is not a complete answer, but the following result might useful. Basically it says that if the poles of $f(x)$ are negative integers, then the time to compute $\sum_{k=0}^{p-1}f(k)$ is bounded by a constant times the time to compute binomial coefficients ${kp\choose p}$ modulo powers of $p$.
Claim: Suppose $f(x)\in\mathbb{Q}(x)$ has poles contained in the negative integers, and let $n\in\mathbb{Z}$ be given. Then there is a polynomial in $p$, $p^{-1}$, and binomial coefficients ${kp\choose p}$, whose coefficients are rational numbers independent of $p$, which is congruent to $\sum_{m=0}^{p-1}f(m)$ modulo $p^n$ for all sufficiently large $p$.
Proof sketch: Using a partial fraction decomposition, it suffices to prove the claim for $f(x)=x^b$, $b\geq 0$, and for $f(x)=(x+a)^{-b}$, $a,b\geq 1$. For $f(x)=x^b$, we can use the formula for power sums, and the sum is a polynomial in $p$. For $f(x)=(x+a)^{-b}$, the sum agrees, up to a rational function of $p$ (which can be expanded in a Laurent series), with
$$
\sum_{j=1}^{p-1}\frac{1}{j^b}.
$$
Write $e_{i,p}$ for the $i$-th elementary symmetric function evaluated at $1^{-1},2^{-1},\ldots,(p-1)^{-1}$. For every integer $j$, we have
$$
\begin{align*}
{jp\choose p}&=j\sum_{i=0}^p (j-1)^i p^i e_{i,p}\\
&\equiv j\sum_{i=0}^{n+b-1}(j-1)^i p^i e_{i,p} \mod p^{n+b}.
\end{align*}
$$
This is a linear system whose coefficient matrix is Vandermonde, hence we can solve for the terms $p^i e_{i,p}$ modulo $p^{n+b}$ as a linear combination of binomial coefficients ${jp\choose p}$, and the coefficients are independent of $p$. For $i\leq b$, we get an expression for $e_{i,p}$ modulo $p^n$ as a linear combination of ${jp\choose p}$. Finally, we can use Newton's formula to expression $\sum_{j=1}^{p-1}\frac{1}{j^b}$ as a polynomial (with coefficients independent of $p$) in the $e_{i,p}$, $1\leq i\leq b$.

I wrote a Sage script to compute these expressions. A random example: for all $p$ sufficiently large,
$$
\sum_{n=0}^{p-1}\frac{n}{(n+1)^2(n+2)}=
-2 \, p^{2} - \frac{1}{36} \, {\left(\frac{{3 p \choose p}}{p} - \frac{6 \, {2 p \choose p}}{p} + \frac{9}{p}\right)}^{2} + 2 \, p + \frac{{3 p \choose p}}{3 \, p^{2}} - \frac{{2 p \choose p}}{p^{2}}+O(p^3).
$$

I'll mention that such sums can also be written in terms of $p$-adic zeta values. My software does more optimization here, so the expressions tend to be simpler. For example:
$$
\sum_{n=0}^{p-1}\frac{n}{(n+1)^2(n+2)}=-p^{-2}+2p-2p\zeta_p(3)-2p^{2}+2p^{3}-4p^{3}\zeta_p(5)-2p^4 + O(p^{5}).
$$
