幂之和和

令$f_k(n)$为前$n$个正整数的$k$次方之和，例如$f_2(10)=\sum_{n=1}^{10}n^2=385$。

再令$S_k(n)=\sum_{i=1}^nf_k(i)$，例如$S_4(100)=35375333830$。

对所有在$2\times 10^9$与$2\times10^9+2000$之间的素数$p$，求$S_{10000}(10^{12}) \bmod p$之和。

\solution 注意到 \begin{align*} S_k(n)&=\sum_{i=1}^n\sum_{j=1}^i j^k, \\ &=\sum_{i=1}^n (n-i+1)\cdot i^k, \\ &=(n+1)\cdot\sum_{i=1}^n i^k-\sum_{i=1}^n i^{k+1}, \\ &=(n+1)f_k(n)-f_{k+1}(n). \end{align*} 而 $$f_k(n)=\sum_{i=1}^ni^k=\frac{n^{k+1}}{k+1}+\frac{1}{2}n^k+\sum_{i=2}^{k}\frac{B_i}{i!}\cdot\frac{k!}{(k-i+1)!}\cdot n^{k-i+1}.$$ 其中$B_i$是第$i$个Bernoulli数。此外，不难验证当$k$与$p-1$互素时有$\sum_{i=0}^{p-1}i^k$模$p$余$0$。

因此若记$r=n \bmod p$，则有 $$f_k(n)=f_k(r)=(-1)^kf(p-1-r) \pmod p.$$ 接下来直接计算即可得结果$106650212746$。

本题难度:

解答

注意到 \begin{align*} S_k(n)&=\sum_{i=1}^n\sum_{j=1}^i j^k, \\ &=\sum_{i=1}^n (n-i+1)\cdot i^k, \\ &=(n+1)\cdot\sum_{i=1}^n i^k-\sum_{i=1}^n i^{k+1}, \\ &=(n+1)f_k(n)-f_{k+1}(n). \end{align*} 而 $$f_k(n)=\sum_{i=1}^ni^k=\frac{n^{k+1}}{k+1}+\frac{1}{2}n^k+\sum_{i=2}^{k}\frac{B_i}{i!}\cdot\frac{k!}{(k-i+1)!}\cdot n^{k-i+1}.$$ 其中$B_i$是第$i$个Bernoulli数。此外，不难验证当$k$与$p-1$互素时有$\sum_{i=0}^{p-1}i^k$模$p$余$0$。

因此若记$r=n \bmod p$，则有 $$f_k(n)=f_k(r)=(-1)^kf(p-1-r) \pmod p.$$ 接下来直接计算即可得结果$106650212746$。

注：因使用sympy作素性检验，以下代码为Python 3，代码中还打印了进度信息。

import sympy

K=10000
N=10**12
L=2000000000
R=2000002000

res=0
for p in range(L,R+1):
    if sympy.isprime(p): 
        r=N%p
        m=p-1-r
        a=b=0
        for i in range(1,m+1):
            x=pow(i,K,p)
            a=(a+x)%p
            b=(b+x*i)%p
        res+=(-a*(r+1)-b)%p
    print(f"p={p} checked, current result: {res}")

print(res)