今回は下記の「べき乗の和」を色々な方法で求めてみたいと思います。
\begin{align}
\sum_{k=1}^n k^p= 1^p+2^p+3^p+\cdots +n^p \tag{1}
\end{align}
ファウルハーバーの公式
Faulhaber's formula (Wikipedia)によるとファウルハーバーの公式は以下のようになります。
\begin{align}
\sum_{k=1}^n k^p=\frac{1}{p+1} \sum_{r=0}^{p} \binom{p+1}{r} B_r n^{p+1-r} \tag{2}
\end{align}
ここで$B_r$はベルヌーイ数です。これだけだと、なかなか分かりづらいですが、そのままコードにして$n$の多項式の係数を求めます。sympyのbinomialとbernoulliを使っています。
import sympy as sp
N = 6
for p in range(1,N+1):
f = [sp.binomial(p+1,r)*sp.bernoulli(r)/(p+1) for r in range(p+1)]
print(p, f)
# 1 [1/2, 1/2]
# 2 [1/3, 1/2, 1/6]
# 3 [1/4, 1/2, 1/4, 0]
# 4 [1/5, 1/2, 1/3, 0, -1/30]
# 5 [1/6, 1/2, 5/12, 0, -1/12, 0]
# 6 [1/7, 1/2, 1/2, 0, -1/6, 0, 1/42]
これによってお馴染みの２乗・３乗の和の公式はこの係数を使って表すことができます。
\begin{align}
\sum_{k=1}^n k^2=\frac1{3}n^3 + \frac1{2}n^2 + \frac1{6}n=\frac1{6}(n+1)(2n+1) \\
\sum_{k=1}^n k^3=\frac1{4}n^4 + \frac1{2}n^3 + \frac1{4}n^2=\lbrace \frac1{2}n(n+1)\rbrace^2 \\
\end{align}
このプログラムはsympyを使うと簡単に書けますが、ベルヌーイ数は自分で書くと結構たいへんです。
Wikipediaをよく見たらベルヌーイ数を使わないで行列をつかって係数を求める方法がありましたので紹介します。
行列でファウルハーバー係数を導出
結論から言うとパスカルの三角数の最後の１を除いて、$i$行$j$列の$(i+j)$が奇数の数にマイナスをつけた以下のような行列を$\bar{A_7}$とします。
\bar{A_7} =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
-1 & 2 & 0 & 0 & 0 & 0 & 0\\
1 & -3 & 3 & 0 & 0 & 0 & 0\\
-1 & 4 & -6 & 4 & 0 & 0 & 0\\
1 & -5 & 10 & -10 & 5 & 0 & 0 \\
-1 & 6 & -15 & 20 & -15 & 6 & 0\\
1 & -7 & 21 & -35 & 35 & -21 & 7
\end{pmatrix}
この逆行列を計算すると驚くことにファウルハーバー係数が現れます。
G_7 =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
\frac1{2} & \frac1{2} & 0 & 0 & 0 & 0 & 0\\
\frac1{6} & \frac1{2} & \frac1{3} & 0 & 0 & 0 & 0\\
0 & \frac1{4} & \frac1{2} & \frac1{4} & 0 & 0 & 0\\
-\frac1{30} & 0 & \frac1{3} & \frac1{2} & \frac1{5} & 0 & 0\\
0 & -\frac1{12} & 0 & \frac{5}{12} & \frac1{2} & \frac1{6} & 0\\
\frac1{42} & 0 & -\frac1{6} & 0 & \frac1{2} & \frac1{2} & \frac1{7}\\
\end{pmatrix}
この詳細を知りたい方はFaulhaber's formula#Matrix formを参照して下さい。
これをコードにします。sympyのbinomialとbernoulliを使っていないことに注目して下さい。ただし逆行列はsympyのMatrixの機能を使用してます。
N = 7
A, b = [], [1]
for j in range(1,N+1):
b = [1]+[(b[i]+b[i+1]) for i in range(len(b)-1)]+[1] # binomia coefficient (Pascal's triangle)
A += [[b[i]*(-1)**(i+1+j) for i in range(len(b)-1)] + [0]*(N+1-len(b))] # change sign, if (i+1+j) is odd
A = sp.Matrix(A) # list -> Matrix
G = sp.Matrix(A)**-1 # inverse Matrix
for i in range(1, N):
print(i, list(reversed(list(G[i,0:(i+1)]))))
# 1 [1/2, 1/2]
# 2 [1/3, 1/2, 1/6]
# 3 [1/4, 1/2, 1/4, 0]
# 4 [1/5, 1/2, 1/3, 0, -1/30]
# 5 [1/6, 1/2, 5/12, 0, -1/12, 0]
# 6 [1/7, 1/2, 1/2, 0, -1/6, 0, 1/42]
これでファウルハーバー係数は割と簡単に計算できることが分かりました。ただこれが威力を発揮するのは$n \ge 10^6$の場合ですが、この場合３乗だと結果が$\ge 10^{18}$となり余り実用的ではないです。
従って次回ではべき乗の和の剰余を求める方法を考えていきたいと思います。
（開発環境：Google Colab）