$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$?

• You say "Let $f$ be a rational polynomial." Do you mean a rational function, in view of your motivation? – David E Speyer Sep 9 '18 at 13:09
• Oops, fixed that. – Itai Bar-Natan Sep 9 '18 at 13:28
• Some expressions for the harmonic sum: $H_{p-1}\equiv ({2p\choose p}-2)/(4p)\equiv -p^2 B_{p-3}/3\equiv -p^2\zeta_p(3)\mod p^3$. I'm not sure if any of these are fast to compute, though. – Julian Rosen Sep 9 '18 at 15:23
• I don't have a copy with me right now, but maybe the chapter on $p$-adic $L$-functions in volume II of Henri Cohen's book on Analytic number theory might help, including the exercises. – EFinat-S Sep 9 '18 at 16:32
• @Julian: factorials can be computed in square root time using fast polynomial evaluation, so the binomial coefficient identity gives a $O(\sqrt{p})$ algorithm :) My guess is that this can be achieved for general rational functions. – Dror Speiser Sep 9 '18 at 18:55

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