スターリング数と冪和公式
著者:梅谷 武
語句:昇階乗,降階乗,第1種スターリング数,第2種スターリング数
スターリング数と冪和公式について述べる。
作成:2006-05-26
更新:2021-03-17
4.2節では三角数や四面体数等の重和が算術三角形を利用することで計算できることを示しました。ここでは、同じように冪和が算術三角形を利用することで計算できることを示します。
同じように有理数体上の多項式環
ℚ[X]において、正の整数
nに対して不定元
Xに対する
k次の昇階乗を
(4.3) | Xn = X(X+1)⋯(X+n-1) (n個の積)
|
|
n次の降階乗を
(4.4) | Xn = X(X-1)⋯(X-n+1) (n個の積)
|
|
と定義します。
0に対しては便宜上
X0 = X0 = 1と定義します。
算術三角形を階乗公式によって階乗で表現してみましょう。
重和を考えるときは算術三角形の横方向の数列を考えましたが、今度は縦方向の数列に着目してみます。例えば第3列
(C(i,3))i ∈ℕは、
| | i ∈ℕ
|
という数列になっています。この数列の第n項までの和は(P3)により
となりますから、
が得られます。この式は立方数の和の公式
と似ています。これは冪乗
i3と階乗
(i+1)(i+2)(i+3)が近い関係にあることによります。階乗の和の公式がすでにわかっていますので、冪乗と階乗との関係がわかれば、それを使って冪和の公式が得られそうです。そこで、まず冪乗と階乗の関係について考えていくことにしましょう。
Xnは
n次の多項式ですから、これを展開したときの
k次の係数を
と書くことにすると
と展開することができます。これにより自然数
n,k(n ≧ k ≧ 0)に対して数
を定義することができます。これを
第1種スターリング数だい1しゅすたーりんぐすう, Stirling cycle numberといいます。第1種スターリング数は負にはならず、
であることはすぐにわかります。また、
Xnと
Xnの展開式における係数は符号が異なるだけで
と書けることもわかります。第1種スターリング数は組み合わせと同じように加法公式が存在します。
任意の自然数
n,k(n>k>0)に対して次の式が成り立つ。
証明
Xn+1 = (X + n)Xn = X ⋅ Xn + n ⋅ Xnより、
となり
n>k>0なる
kについて係数を比較することによってわかる。■
この加法公式を使って、第1種スターリング数の三角形を作ってみましょう。作り方は三角形の細胞をパスカルの三角形と同じく
C(i,j)と書くことにすれば、
とします。言い換えれば、パスカルの算術三角形における
を
で置き換えたものです。まず、
とし、
n>0に対して
とすることによって第0行と第0列が埋まります。次に加法公式によって上の行の左から順番に
C(i,j) = C(i,j-1) + (i+j-1)C(i-1,j)
|
によって細胞を埋めていき、行末になったら下の行の左から同じように埋めていくという手順ですべての細胞を埋めることができます。
第1種スターリング数の三角形だい1しゅすたーりんぐすうのさんかくけいを一見して気が付くのは第1行が三角数になっていることと第1列が
と狭義の階乗になっていることでしょう。さらに底辺の和は次のようにちょうど狭義の階乗になっています。
これは第1種スターリング数を組み合わせ論的に解釈することによって比較的容易に証明することができますが、これについては章末の研究で詳しく考えることにします。
このようにしてできあがった第1種スターリング数の三角形が正しいかどうかを、多項式を展開した係数と比較することによって確認してみましょう。
from Poly1 import *
def RisingFactorial( n ):
if n == 0: return Poly1( [1] )
a = b = Poly1( [0,1] )
for i in range( n - 1 ):
b = b + Poly1( [1] )
a = a * b
return a
def FallingFactorial( n ):
if n == 0: return Poly1( [1] )
a = b = Poly1( [0,1] )
for i in range( n - 1 ):
b = b - Poly1( [1] )
a = a * b
return a
for i in range( 10 ):
print "i = ", i
print RisingFactorial( i )
print FallingFactorial( i )
i = 0
1
1
i = 1
X
X
i = 2
X^2 + X
X^2 - X
i = 3
X^3 + 3X^2 + 2X
X^3 - 3X^2 + 2X
i = 4
X^4 + 6X^3 + 11X^2 + 6X
X^4 - 6X^3 + 11X^2 - 6X
i = 5
X^5 + 10X^4 + 35X^3 + 50X^2 + 24X
X^5 - 10X^4 + 35X^3 - 50X^2 + 24X
i = 6
X^6 + 15X^5 + 85X^4 + 225X^3 + 274X^2 + 120X
X^6 - 15X^5 + 85X^4 - 225X^3 + 274X^2 - 120X
i = 7
X^7 + 21X^6 + 175X^5 + 735X^4 + 1624X^3 + 1764X^2 + 720X
X^7 - 21X^6 + 175X^5 - 735X^4 + 1624X^3 - 1764X^2 + 720X
i = 8
X^8 + 28X^7 + 322X^6 + 1960X^5 + 6769X^4 + 13132X^3 + 13068X^2 + 5040X
X^8 - 28X^7 + 322X^6 - 1960X^5 + 6769X^4 - 13132X^3 + 13068X^2 - 5040X
i = 9
X^9 + 36X^8 + 546X^7 + 4536X^6 + 22449X^5 + 67284X^4 + 118124X^3 + 109584X^2 + 40320X
X^9 - 36X^8 + 546X^7 - 4536X^6 + 22449X^5 - 67284X^4 + 118124X^3 - 109584X^2 + 40320X
require 'mathn'
require './poly1'
def RisingFactorial(n)
if n == 0
return Poly1([1])
end
a = b = Poly1([0,1])
for i in 0...(n-1)
b = b + Poly1([1])
a = a * b
end
return a
end
def FallingFactorial(n)
if n == 0
return Poly1([1])
end
a = b = Poly1([0,1])
for i in 0...(n-1)
b = b - Poly1([1])
a = a * b
end
return a
end
for i in 0...10
print "i = ", i, "\n"
print RisingFactorial(i), "\n"
print FallingFactorial(i), "\n"
end
i = 0
1
1
i = 1
X
X
i = 2
X^2 + X
X^2 - X
i = 3
X^3 + 3X^2 + 2X
X^3 - 3X^2 + 2X
i = 4
X^4 + 6X^3 + 11X^2 + 6X
X^4 - 6X^3 + 11X^2 - 6X
i = 5
X^5 + 10X^4 + 35X^3 + 50X^2 + 24X
X^5 - 10X^4 + 35X^3 - 50X^2 + 24X
i = 6
X^6 + 15X^5 + 85X^4 + 225X^3 + 274X^2 + 120X
X^6 - 15X^5 + 85X^4 - 225X^3 + 274X^2 - 120X
i = 7
X^7 + 21X^6 + 175X^5 + 735X^4 + 1624X^3 + 1764X^2 + 720X
X^7 - 21X^6 + 175X^5 - 735X^4 + 1624X^3 - 1764X^2 + 720X
i = 8
X^8 + 28X^7 + 322X^6 + 1960X^5 + 6769X^4 + 13132X^3 + 13068X^2 + 5040X
X^8 - 28X^7 + 322X^6 - 1960X^5 + 6769X^4 - 13132X^3 + 13068X^2 - 5040X
i = 9
X^9 + 36X^8 + 546X^7 + 4536X^6 + 22449X^5 + 67284X^4 + 118124X^3 + 109584X^2 + 40320X
X^9 - 36X^8 + 546X^7 - 4536X^6 + 22449X^5 - 67284X^4 + 118124X^3 - 109584X^2 + 40320X
第1種スターリング数によって、階乗を冪乗の多項式として表わすことができました。今度は逆に冪乗を階乗の多項式として表わすことを考えてみましょう。冪乗は
X2 = X(X-1) + X, X3 = X(X-1)(X-2) + 3X(X-1) + X
|
というように降階乗の多項式で表わすことができます。これを一般化して
により自然数
n,k(n≧k≧0)に対して数
を定義します。これを
第2種スターリング数だい2しゅすたーりんぐすう, Stirling subset numberといいます。第2種スターリング数が負にはならず、
であることはすぐにわかります。第2種スターリング数は次の加法公式を満たします。
任意の自然数
n,k(n>k>0)に対して次の式が成り立つ。
証明
まず、
{Xk | k∈ℕ}が
ℚ[X]において
ℚ上線形独立であることを注意しておく。すなわち、
ならばすべての
akは
0である。
より、
であるから、
となり
n>k>0なる
kについて係数を比較することによってわかる。■
この加法公式を使って第2種スターリング数の三角形を作ってみましょう。作り方は第1種スターリング数の三角形と同様に
とします。
n>0に対して
とし、加法公式
C(i,j) = C(i,j-1) + j C(i-1,j)
|
によって細胞を埋めていきます。
第2種スターリング数の三角形においても第1行が三角数になっています。また、第2列が
となっています。これも組み合わせ論的に解釈することによって容易に証明できますので、章末の研究で考えることにします。
冪乗の昇階乗に関する多項式は符号を考慮するだけで得られます。
第2種スターリング数の三角形が正しいかどうかを確認してみましょう。
from Poly1 import *
def RisingFactorial( n ):
if n == 0: return Poly1( [1] )
a = b = Poly1( [0,1] )
for i in range( n - 1 ):
b = b + Poly1( [1] )
a = a * b
return a
def FallingFactorial( n ):
if n == 0: return Poly1( [1] )
a = b = Poly1( [0,1] )
for i in range( n - 1 ):
b = b - Poly1( [1] )
a = a * b
return a
stir2 = [ [1,1],\
[1,3,1],\
[1,6,7,1],\
[1,10,25,15,1],\
[1,15,65,90,31,1],\
[1,21,140,350,301,63,1],\
[1,28,266,1050,1701,966,127,1],\
[1,36,462,2646,6951,7770,3025,255,1] ]
for n in range( 2, 10, 1 ):
coef = stir2[n-2]
i = n; a = b = Poly1([0])
sign = p = Poly1([1]); m = Poly1([-1])
for k in coef:
a += Poly1([k])*FallingFactorial( i )
b += sign * Poly1([k])*RisingFactorial( i )
if sign == p: sign = m
else: sign = p
i = i - 1
print "X^%d = " % n, a, ", ", b
X^2 = X^2 , X^2
X^3 = X^3 , X^3
X^4 = X^4 , X^4
X^5 = X^5 , X^5
X^6 = X^6 , X^6
X^7 = X^7 , X^7
X^8 = X^8 , X^8
X^9 = X^9 , X^9
require 'mathn'
require './poly1'
def RisingFactorial(n)
if n == 0
return Poly1([1])
end
a = b = Poly1([0,1])
for i in 0...(n-1)
b = b + Poly1([1])
a = a * b
end
return a
end
def FallingFactorial(n)
if n == 0
return Poly1([1])
end
a = b = Poly1([0,1])
for i in 0...(n-1)
b = b - Poly1([1])
a = a * b
end
return a
end
stir2 = [ [1,1],
[1,3,1],
[1,6,7,1],
[1,10,25,15,1],
[1,15,65,90,31,1],
[1,21,140,350,301,63,1],
[1,28,266,1050,1701,966,127,1],
[1,36,462,2646,6951,7770,3025,255,1] ]
for n in 2...10
coef = stir2[n-2]
i = n; a = b = Poly1([0])
sign = p = Poly1([1]); m = Poly1([-1])
for k in coef
a += Poly1([k]) * FallingFactorial(i)
b += sign * Poly1([k]) * RisingFactorial(i)
if sign == p
sign = m
else
sign = p
end
i = i - 1
end
printf "X^%d = %s, %s\n", n, a, b
end
X^2 = X^2 , X^2
X^3 = X^3 , X^3
X^4 = X^4 , X^4
X^5 = X^5 , X^5
X^6 = X^6 , X^6
X^7 = X^7 , X^7
X^8 = X^8 , X^8
X^9 = X^9 , X^9
さて、ここまでの準備によって
Pk(n)を
nの多項式として書き下すことができるようになりました。まず、冪乗を第2種スターリング数によって昇階乗に書き換えます。
昇階乗の指数が同じものの和をまとめ、パスカルの帰結(P3)を適用します。
| | k ∑ j=0 | (-1)k-j | | | | |
|
|
|
| | |
この昇階乗を第1種スターリング数によって冪乗に書き換えて、同じ次数の単項式の和と
して整理します。
| | |
| | k+1 ∑ t=0 | | | | nt |
|
|
|
ここで
t>j+1のとき
としていることに注意してください。これで
Pk(n)を
nの多項式として書くことができました。
であることを考慮してこれを命題としておきましょう。
自然数列の
k次の冪和は次のように表すことができる。
これからベルヌーイ数をスターリング数を使って求める公式が導かれます。
証明
ベルヌーイ数
Bkは
Pk(X)の一次の係数であるから、
この式に
を代入する。■
次にスターリング数による冪和公式があっているかどうかを確かめてみますが、計算の都合上、
と定義することにします。
from Poly1 import *
# calculate tables
tr1 = []
tr2 = []
for i in range( 10 ):
tr1.append( range( 0, 10, 1 ) )
tr2.append( range( 0, 10, 1 ) )
for j in range( 10 ):
tr1[0][j] = tr2[0][j] = 1
for i in range( 1, 10, 1 ):
for j in range( 1, 10, 1 ):
tr1[i][j] = tr1[i][j-1] + (i+j-1) * tr1[i-1][j]
tr2[i][j] = tr2[i][j-1] + j * tr2[i-1][j]
def s( n, k ):
if not isinstance( n, int ): raise TypeError
if not isinstance( k, int ): raise TypeError
if ( k < 0 ) or ( n < 0 ): raise TypeError
if ( k > 9 ) or ( n > 9 ): raise TypeError
if ( k > n ): return 0
else:
if ( n == 0 ) and ( k == 0 ): return 1
elif n == k: return 1
elif k == 0: return 0
else: return tr1[n-k][k]
def S( n, k ):
if not isinstance( n, int ): raise TypeError
if not isinstance( k, int ): raise TypeError
if ( k < 0 ) or ( n < 0 ): raise TypeError
if ( k > 9 ) or ( n > 9 ): raise TypeError
if ( k > n ): return 0
else:
if ( n == 0 ) and ( k == 0 ): return 1
elif n == k: return 1
elif k == 0: return 0
else: return tr2[n-k][k]
def coef( k, t ):
if not isinstance( k, int ): raise TypeError
if not isinstance( t, int ): raise TypeError
if ( k < 0 ) or ( t < 1 ): raise TypeError
else:
a = Fraction(0,1)
for j in range( t - 1, k + 1, 1 ):
b = ((-1)**(k-j)) * S(k,j) * s(j+1,t)
a += Fraction(b,j+1)
return a
for k in range( 1, 9, 1 ):
list = [ Fraction(0,1) ]
for t in range( 1, k + 2, 1 ):
list.append( coef( k, t ) )
p = Poly1( list )
print p
for n in [ 3, 5, 7 ]:
sum = 0
for i in range( 1, n + 1, 1 ):
sum += i**k
print "n =", n, sum, p.Horner(n)
1/2X^2 + 1/2X
n = 3 6 6
n = 5 15 15
n = 7 28 28
1/3X^3 + 1/2X^2 + 1/6X
n = 3 14 14
n = 5 55 55
n = 7 140 140
1/4X^4 + 1/2X^3 + 1/4X^2
n = 3 36 36
n = 5 225 225
n = 7 784 784
1/5X^5 + 1/2X^4 + 1/3X^3 - 1/30X
n = 3 98 98
n = 5 979 979
n = 7 4676 4676
1/6X^6 + 1/2X^5 + 5/12X^4 - 1/12X^2
n = 3 276 276
n = 5 4425 4425
n = 7 29008 29008
1/7X^7 + 1/2X^6 + 1/2X^5 - 1/6X^3 + 1/42X
n = 3 794 794
n = 5 20515 20515
n = 7 184820 184820
1/8X^8 + 1/2X^7 + 7/12X^6 - 7/24X^4 + 1/12X^2
n = 3 2316 2316
n = 5 96825 96825
n = 7 1200304 1200304
1/9X^9 + 1/2X^8 + 2/3X^7 - 7/15X^5 + 2/9X^3 - 1/30X
n = 3 6818 6818
n = 5 462979 462979
n = 7 7907396 7907396
require 'mathn'
require './poly1'
# calculate tables
$tr1 = Array.new(10) { Array.new(10) { 0 } }
$tr2 = Array.new(10) { Array.new(10) { 0 } }
for j in 0...10
$tr1[0][j] = 1; $tr2[0][j] = 1
end
for i in 1...10
for j in 1...10
$tr1[i][j] = $tr1[i][j-1] + (i+j-1) * $tr1[i-1][j]
$tr2[i][j] = $tr2[i][j-1] + j * $tr2[i-1][j]
end
end
def s(n, k)
if ( ! n.kind_of?(Integer) ) or ( ! k.kind_of?(Integer) )
raise TypeError
end
if ( k < 0 ) or ( n < 0 ) or ( k > 9 ) or ( n > 9 )
raise TypeError
end
if ( k > n )
return 0
else
if ( n == 0 ) and ( k == 0 )
return 1
elsif n == k
return 1
elsif k == 0
return 0
else
return $tr1[n-k][k]
end
end
end
def S(n, k)
if ( ! n.kind_of?(Integer) ) or ( ! k.kind_of?(Integer) )
raise TypeError
end
if ( k < 0 ) or ( n < 0 ) or ( k > 9 ) or ( n > 9 )
raise TypeError
end
if ( k > n )
return 0
else
if ( n == 0 ) and ( k == 0 )
return 1
elsif n == k
return 1
elsif k == 0
return 0
else
return $tr2[n-k][k]
end
end
end
def coef(k, t)
if ( ! k.kind_of?(Integer) ) or ( ! t.kind_of?(Integer) )
raise TypeError
end
if ( k < 0 ) or ( t < 1 )
raise TypeError
else
a = 0
for j in (t-1)..k
b = ((-1)**(k - j)) * S(k,j) * s(j+1,t)
a += b / (j + 1)
end
return a
end
end
for k in 1...9
list = [0]
for t in 1..(k+1)
list.push(coef(k, t))
end
p = Poly1(list)
print p, "\n"
for n in [3, 5, 7]
sum = 0
for i in 1..n
sum += i**k
end
printf "n = %d %d %d\n", n, sum, p.Horner(n)
end
end
1/2X^2 + 1/2X
n = 3 6 6
n = 5 15 15
n = 7 28 28
1/3X^3 + 1/2X^2 + 1/6X
n = 3 14 14
n = 5 55 55
n = 7 140 140
1/4X^4 + 1/2X^3 + 1/4X^2
n = 3 36 36
n = 5 225 225
n = 7 784 784
1/5X^5 + 1/2X^4 + 1/3X^3 - 1/30X
n = 3 98 98
n = 5 979 979
n = 7 4676 4676
1/6X^6 + 1/2X^5 + 5/12X^4 - 1/12X^2
n = 3 276 276
n = 5 4425 4425
n = 7 29008 29008
1/7X^7 + 1/2X^6 + 1/2X^5 - 1/6X^3 + 1/42X
n = 3 794 794
n = 5 20515 20515
n = 7 184820 184820
1/8X^8 + 1/2X^7 + 7/12X^6 - 7/24X^4 + 1/12X^2
n = 3 2316 2316
n = 5 96825 96825
n = 7 1200304 1200304
1/9X^9 + 1/2X^8 + 2/3X^7 - 7/15X^5 + 2/9X^3 - 1/30X
n = 3 6818 6818
n = 5 462979 462979
n = 7 7907396 7907396
数 学
昇階乗 しょうかいじょう, rising factorial
昇階乗、降階乗、第1種スターリング数、第2種スターリング数を表わすときに『コンピュータの数学』において採用された記号を使用している。
降階乗 こうかいじょう, falling factorial
第1種スターリング数 だい1しゅすたーりんぐすう, Stirling cycle number
第1種スターリング数の三角形 だい1しゅすたーりんぐすうのさんかくけい
関孝和の『括要算法』にも昇階乗を展開することにより得られた第1種スターリング数の三角形が書かれているが、関がこれと関・ベルヌーイ数との関係を知っていたかどうかは不明である。
第2種スターリング数 だい2しゅすたーりんぐすう, Stirling subset number