前回のその２ではべき乗の和の剰余計算を３つの方法で計算しました。
今回は以下のような「べき乗の和の和」$SS(p,n) \pmod{m}$を考えます

\begin{align}
&S(p,n) = \sum_{k=1}^{n} k^p　\pmod{m}  \tag{1}\\
&SS(p, n) = \sum_{i=1}^{n}S(p,i) =\sum_{i=1}^{n}\sum_{k=1}^i k^p　\pmod{m} \tag{2} 
\end{align}

【例題１】$SS(5, 10^6) \pmod{(10^9+7)} $を求めよ

【コード#1】素直に計算

今回はすべての$1 \le i \le n$を求める必要があるので、ファウルハーバー係数を使った計算は効果的ではないので、素直に計算しました。(itertoolsのaccumulate利用)。

from itertools import accumulate
def SS1(p, n, m):
  return sum(list(accumulate([pow(k, p, m) for k in range(1, n+1)])))%m

p, n, m = 5, 10**6, 10**9+7
print(f"Answer: {SS1(p,n,m)}")
# 803608076

式(２)を変形してシグマを１重にする

式（２）を行列で書くと

SS(p, n) =  
\begin{pmatrix} 
 1^p \\
1^p & 2^p \\
1^p & 2^p & 3^p \\ 
1^p & 2^p & 3^p & 4^p \\
\cdots \\
1^p & 2^p & 3^p & 4^p &\cdots & n^p
\end{pmatrix}

これを各列で見ると

\begin{align}
SS(p, n) &= n \cdot 1^p + (n-1) \cdot 2^p + (n-2)\cdot 3^p \cdots 1 \cdot n^p\\
& = \sum_{k=1}^{n} (n+1-k) \cdot k^{p} \\
&= (n+1) \cdot S(p,n) - S(p+1, n) 　\tag{3}
\end{align}

となり結局$S(p,n)$を求める問題に帰着できます。

【例題２】$SS(10^3, 10^{14} \pmod{(10^9+7)} $を求めよ

式（３）を使ってSS(p,n)を求める

前回作った関数$S3$を使えば簡単に求めることができました。

def SS3(p, n, m):
  return ((n+1)*S3(p,n,m)-S3(p+1,n,m))%m

p, n, m = 10**3, 10**14, 10**9+7
print(f"Answer : {SS3(p, n, m)}")
# Answer : 419798333

（開発環境：Google Colab）

この考え方はProject Euler, Problem 487: Sums of Power Sumsを解くのに役に立ちます。