点を多角形の形に並べたとき、その総数を「多角数 (polygonal number)」といいます。三角形に配置したものを「三角数 (triangular number)」、四角形に配置したものを「四角数 (square number)」、五角形に配置したものを「五角数 (pentagonal number)」といいます。三角数、四角数、五角数を下図に示します。
１ ３ ６ 10 15
● ● ● ● ●
●● ●● ●● ●●
●●● ●●● ●●●
●●●● ●●●●
●●●●●
図 : 三角数
１ ４ ９ 16 25
● ●● ●●● ●●●● ●●●●●
●● ●●● ●●●● ●●●●●
●●● ●●●● ●●●●●
●●●● ●●●●●
●●●●●
図 : 四角数
１ ５ 12 22
● ● ● ●
● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ●
● ● ● ●
● ● ● ●
● ● ● ●
図 : 五角数
それでは問題です。
多角数の点の増分を表に示すと、次のようになります。
n 三角数 四角数 五角数
---+-----------------------------------------------------------
1 | 1 1 1
2 | 3 = 1+2 4 = 1+3 5 = 1+4
3 | 6 = 1+2+3 9 = 1+3+5 12 = 1+4+7
4 | 10 = 1+2+3+4 16 = 1+3+5+7 22 = 1+4+7+10
5 | 15 = 1+2+3+4+5 25 = 1+3+5+7+9 35 = 1+4+7+10+13
6 | 21 = 1+2+3+4+5+6 36 = 1+3+5+7+9+11 51 = 1+4+7+10+13+16
・・・・・・ ・・・・・・・ ・・・・・
n | n(n + 1) / 2 n^2 n(3n - 1) / 2
表を見ればお分かりのように、三角数は公差 1、四角数は公差 2、五角数は公差 3、p 角数は公差 p - 2 の等差数列の和になります。初項を a, 公差を d とすると、等差数列の和 Sn は次式で求めることができます。
Sn = n(2a + (n - 1)d) / 2
公式の導出はあとで説明します。
a = 1, d = p - 2 を代入して計算すると、多角数 Pp,n は次式で求めることができます。
Pp,n = ((p - 2)n^2 - (p - 4)n) / 2
この式を Python でプログラムすると、次のようになります。
リスト : 多角数
def polygonal_number(p, n):
return ((p - 2) * n * n - (p - 4) * n) / 2
for p in xrange(3, 9):
print p, ":",
for x in xrange(1, 17):
print polygonal_number(p, x),
print
それでは実行してみましょう。
3 : 1 3 6 10 15 21 28 36 45 55 66 78 91 105 120 136 4 : 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256 5 : 1 5 12 22 35 51 70 92 117 145 176 210 247 287 330 376 6 : 1 6 15 28 45 66 91 120 153 190 231 276 325 378 435 496 7 : 1 7 18 34 55 81 112 148 189 235 286 342 403 469 540 616 8 : 1 8 21 40 65 96 133 176 225 280 341 408 481 560 645 736
p 角数の初項は 1 で、第 2 項は p になります。参考 URL 1 によると、多角数にはいろいろな面白い性質があるようです。
次のように、一定の差で並んだ数列を「等差数列」といいます。
a, a + d, a + 2d, a + 3d, ..., a + (n - 1)d
a を「初項」、d を「公差」といいます。等差数列の一般項は次の式で表すことができます。
an = a + (n - 1)d
初項から an までの和 Sn は次の式で求めることができます。
Sn = n(2a + (n - 1)d) / 2
初項を 1, 公差 を 1 とすると、1 から n までの和は n(n + 1)/ 2 となります。
この公式は次のように導出することができます。
Sn = a + (a + d) + ,,, + (a + (n - 2)d) + (a + (n - 1)d) Sn = (a + (n - 1)d) + (a + (n - 2)d) + ... + (a + d) + a 足し算すると 2Sn = (2a + (n - 1)d) + (2a + (n - 1)d) + ... (2a + (n - 1)d) + (2a + (n - 1)d) 2Sn = n(2a + (n - 1)d) Sn = n(2a + (n - 1)d)/2
このように、右辺を逆順に並べ替えて足し算すると、2a + (n - 1)d が n 個並ぶことになります。あとは、これを 2 で割り算すればいいわけです。
整数 m が p 角数か判定するプログラムも簡単です。次式で n が整数値になれば、m は p 角数の第 n 項であることがわかります。
(p - 4) + √((p - 4)^2 + 8(p - 2)m)
n = -----------------------------------
2(p - 2)
リスト : 多角数の判定
def is_polygonal_number(p, m):
a = p - 2
b = p - 4
c = b * b + 8 * a * m
d = int(sqrt(c))
return d * d == c and (b + d) % (2 * a) == 0
print is_polygonal_number(3, 55)
print is_polygonal_number(3, 54)
print is_polygonal_number(4, 64)
print is_polygonal_number(4, 65)
print is_polygonal_number(5, 117)
print is_polygonal_number(5, 118)
p - 2 を変数 a に、p - 4 を変数 b に、ルートの中を計算して変数 c にセットします。sqrt() で c の平方根を求め、int() で整数値に変換して d にセットします。あとは、d * d が c に等しく、b + d が 2 * a で割り切れることを確認するだけです。
ただし、sqrt() は整数を浮動小数点数に変換してから平方根を計算するため、大きな整数を与えると誤差により正確な値を計算することはできません。簡単な例を示しましょう。
>>>> sqrt(999999999999999 ** 2) 999999999999999.0 >>>> sqrt(9999999999999999 ** 2) 1e+16
浮動小数点数の精度は 15 桁くらいしかないので、これよりも大きな数の判定に is_polygonal_number() を使うことはできません。大きな数に対応する場合は、ニュートン法などで整数の平方根を求めるようにするといいでしょう。このへんの話は、inamori さんの 数値計算による誤差 - 桃の天然水 が参考になると思います。inamori さんに感謝いたします。
それでは実行してみましょう。
True False True False True False
もちろん、n の値を求めることもできます。興味のある方はプログラムを修正してみてください。
Python でナイーブにプログラムすると次のようになります。
リスト : 平方三角数
def square_triangular_num(n):
m = int(sqrt(n))
for x in xrange(1, m + 1):
if is_polygonal_number(3, x * x):
print x * x,
print
square_triangular_num(1500000)
生成した四角数が三角数になるか is_polygonal_number() でチェックするだけです。実行結果は次のようになります。
1 36 1225 41616 1413721
ところで、参考 URL 3 に掲載されている平方三角数の公式 (一般項と漸化式) を使うと、もっと簡単にプログラムすることができます。
f(0) = 0
f(1) = 1
f(n) = 34 * f(n - 1) - f(n - 2) + 2
平方三角数の漸化式
漸化式をそのままプログラムすると次のようになります。
リスト : 平方三角数 (2)
def square_triangular(n):
if n == 0: return 0
if n == 1: return 1
return 34 * square_triangular(n - 1) - square_triangular(n - 2) + 2
for x in xrange(1, 9):
print square_triangular(x)
二重再帰なので、実行速度は遅いです。興味のある方は繰り返しに変換してください。
実行結果は次のようになります。
1 36 1225 41616 1413721 48024900 1631432881 55420693056
Python でナイーブにプログラムすると次のようになります。
リスト : 三角数で五角数の数
def triangular_pentagonal_num(n):
i = 1
x = 1
while x <= n:
if is_polygonal_number(3, x):
print x,
i += 1
x = polygonal_number(5, i)
print
triangular_pentagonal_num(10000000)
五角数を生成して、それが三角数になるか is_polygonal_number() でチェックするだけです。実行結果は次のようになります。
1 210 40755 7906276
ところで、オンライン整数列大辞典の数列 A069017 によると、三角数かつ五角数の漸化式は次のようになるそうです。
f(0) = 0
f(1) = 1
f(n) = 194 * f(n - 1) - f(n - 2) + 16
三角数かつ五角数の漸化式
これをプログラムすると次のようになります。
リスト : 三角数で五角数の数 (2)
def triangular_pentagonal(n):
if n == 0: return 0
if n == 1: return 1
return 194 * triangular_pentagonal(n - 1) - triangular_pentagonal(n - 2) + 16
for x in xrange(1, 9):
print triangular_pentagonal(x)
二重再帰なので、実行速度は遅いです。興味のある方は繰り返しに変換してください。
実行結果は次のようになります。
1 210 40755 7906276 1533776805 297544793910 57722156241751 11197800766105800
Python でナイーブにプログラムすると次のようになります。
リスト : 四角錐数
def square_pyramidal(n):
a = 0
for m in xrange(1, n + 1): a += m * m
return a
for x in xrange(1, 11):
print square_pyramidal(x),
print
print square_pyramidal(1000000)
1 5 14 30 55 91 140 204 285 385 333333833333500000
次の公式を使うともっと簡単に求めることができます。
四角錐数(n) = n(n + 1)(2n + 1) / 6
四角錐数の公式は次のように導出することができます。
S = 1 + 22 + 32 + ... + n2 U = 1 + 23 + 33 + ... + (n + 1)3 = 1 + (1 + 1)3 + (2 + 1)3 + ... + (n + 1)3 とする。(a + b)3 = a3 + 3a2b + 3ab2 + b3 なので、U は次のように変形できる。 U = (1 + 23 + ... + n3) ; 左辺に移項すると消える + 3(1 + 22 + ... + n2) ; 3S + 3(1 + 2 + ... + n) ; 3n(n + 1) / 2 + (n + 1) (n + 1)3 = 3S + 3n(n + 1)/2 + (n + 1) 両辺を 2 倍して式を変形する。 6S = 2(n + 1)3 - 3n(n + 1) - 2(n + 1) = (n + 1)(2n2 + 4n + 2 - 3n - 2) = (n + 1)(2n2 + n) = n(n + 1)(2n + 1) => S = n(n + 1)(2n + 1) / 6
Python でナイーブにプログラムすると次のようになります。
リスト : 三角錐数
def triangular_pyramidal(n):
a = 0
for m in xrange(1, n + 1): a += m * (m + 1) / 2
return a
for x in xrange(1, 11):
print triangular_pyramidal(x),
print
print triangular_pyramidal(1000000)
1 4 10 20 35 56 84 120 165 220 166667166667000000
次の公式を使うともっと簡単に求めることができます。
三角錐数(n) = n(n + 1)(n + 2) / 6
四角錐数の公式を使うと、三角錐数の公式は次のように求めることができます。
Σ(n(n + 1)/2) = (Σn2 + Σn) / 2
= (n(n + 1)(2n + 1) + 3n(n + 1)) / 12
= (2n3 + 6n2 + 4n) / 12
= n(n2 + 3n + 2) / 6
= n(n + 1)(n + 2) / 6