M.Hiroi's Home Page

Puzzle DE Programming

数字のパズル : 多角数 (polygonal number)

Copyright (C) 2015-2022 Makoto Hiroi
All rights reserved.

パズルの説明

点を多角形の形に並べたとき、その総数を多角数 (polygonal number) といいます。三角形に配置したものを三角数 (triangular number)、四角形に配置したものを四角数 (square number)、五角形に配置したものを五角数 (pentagonal number) といいます。

三角数、四角数、五角数を下図に示します。 

１    ３      ６        10          15
●    ●      ●        ●          ●
     ●●    ●●      ●●        ●●
            ●●●    ●●●      ●●●
                     ●●●●    ●●●●
                                ●●●●●

            図 : 三角数

１   ４      ９        16          25
●  ●●   ●●●   ●●●●   ●●●●●
    ●●   ●●●   ●●●●   ●●●●●
           ●●●   ●●●●   ●●●●●
                    ●●●●   ●●●●●
                               ●●●●●

            図 : 四角数

１       ５             12                      22
●       ●             ●                      ●
      ●    ●       ●    ●               ●      ●
                 ●            ●        ●             ●
       ●  ●         ●              ●      ●           ●
                   ●    ●  ●            ●     ●
                                       ●            ●  ●
                     ● ● ●               ●
                                         ●     ● ●  ●

                                           ● ● ● ●

            図 : 五角数

それでは問題です。

  1. n 番目の m 角数を求めるプログラムを作ってください。
  2. 自然数 n が m 角数か判定するプログラムを作ってください。
  3. 三角数かつ四角数である数を平方三角数 (square triangular number) といいます。
    1,500,000 以下の平方三角数をすべて求めてください。
  4. 10,000,000 以下の三角数かつ五角数である数をすべて求めてください。
  5. 四角数を 1 から順番に足した数列を考えることができます。これを四角錐数 (square pyramidal number) といいます。
    1,000,000 番目の四角錐数を求めてください。
  6. 三角数を 1 から順番に足した数列を考えることができます。これを三角錐数 (triangular pyramidal number) といいます。
    1,000,000 番目の三角錐数を求めてください。

●解答１

多角数の点の増分を表に示すと、次のようになります。 

 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 とすると、等差数列の和 \(S_n\) は次式で求めることができます。

\( S_n = \dfrac{n(2a + (n - 1)d)}{2} \)

公式の導出はあとで説明します。

a = 1, d = p - 2 を代入して計算すると、多角数 \(P_{p,n}\) は次式で求めることができます。

\( P_{p,n} = \dfrac{(p - 2)n^2 - (p - 4)n}{2} \)

この式を Python3 (ver 3.8.10) でプログラムすると、次のようになります。 

リスト : 多角数 (polynum.py)

def polygonal_number(p, n):
    return ((p - 2) * n * n - (p - 4) * n) // 2

def test1():
    for p in range(3, 9):
        print(p, ":", end=' ')
        for x in range(1, 17):
            print(polygonal_number(p, x), end=' ')
        print()

それでは実行してみましょう。 

>>> from polynum import *
>>> test1()
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 『多角数 - Wikipedia』によると、多角数にはいろいろな面白い性質があるようです。

●等差数列の和

次のように、一定の差で並んだ数列を「等差数列」といいます。

\(a, \ a + d, \ a + 2d, \ a + 3d, \ \cdots, \ a + (n - 1)d, \ \cdots\)

a を「初項」、d を「公差」といいます。等差数列の一般項は次の式で表すことができます。

\(a_n = a + (n-1)d\)

初項から \(a_n\) までの和 \(S_n\) は次の式で求めることができます。

\(S_n = \displaystyle \sum_{i=1}^{n} a_i = \frac{n(2a + (n - 1)d)}{2}\)

初項を 1, 公差 を 1 とすると、1 から n までの和は \(\frac{n(n + 1)}{2}\) となります。

この公式は次のように導出することができます。

\(\begin{eqnarray} S_n &=& a &+& (a + d) &+& \cdots \ &+& (a + (n - 2)d) &+& (a + (n - 1)d) \\ S_n &=& (a + (n - 1)d) &+& (a + (n - 2)d) &+& \cdots \ &+& (a + d) &+& a \end{eqnarray}\)

足し算すると

\(\begin{eqnarray} 2S_n &=& (2a + (n - 1)d) + (2a + (n - 1)d) + \cdots + (2a + (n - 1)d) + (2a + (n - 1)d) \\ 2S_n &=& n(2a + (n - 1)d) \\ S_n &=& \dfrac{n(2a + (n - 1)d)}{2} \end{eqnarray}\)

このように、右辺を逆順に並べ替えて足し算すると、\(2a + (n - 1)d\) が n 個並ぶことになります。あとは、これを 2 で割り算すればいいわけです。

●解答２

整数 m が p 角数か判定するプログラムも簡単です。次式で n が整数値になれば、m は p 角数の第 n 項であることがわかります。

\( n = \dfrac{(p - 4) + \sqrt{(p - 4)^2 + 8(p - 2)m}}{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

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 さんに感謝いたします。

それでは実行してみましょう。 

>>> is_polygonal_number(3, 55)
True
>>> is_polygonal_number(3, 54)
False
>>> is_polygonal_number(4, 64)
True
>>> is_polygonal_number(4, 65)
False
>>> is_polygonal_number(5, 117)
True
>>> is_polygonal_number(5, 118)
False

もちろん、n の値を求めることもできます。興味のある方はプログラムを修正してみてください。

●解答３

Python でナイーブにプログラムすると次のようになります。 

リスト : 平方三角数

def square_triangular_num(n):
    m = int(sqrt(n))
    for x in range(1, m + 1):
        if is_polygonal_number(3, x * x):
            print(x * x, end=' ')
    print()

square_triangular_num(1500000)

生成した四角数が三角数になるか is_polygonal_number() でチェックするだけです。実行結果は次のようになります。 

>>> square_triangular_num(1500000)
1 36 1225 41616 1413721

ところで、参考 URL 『平方三角数 - Wikipedia』に掲載されている平方三角数の公式 (一般項と漸化式) を使うと、もっと簡単にプログラムすることができます。

\(\begin{array}{l} f(0) = 0 \\ f(1) = 1 \\ f(n) = 34 \times f(n - 1) - f(n - 2) + 2 \end{array}\)

平方三角数の漸化式

漸化式をそのままプログラムすると次のようになります。 

リスト : 平方三角数 (2)

def square_triangular(n):
    if n < 2: return n
    return 34 * square_triangular(n - 1) - square_triangular(n - 2) + 2

二重再帰なので、実行速度は遅いです。興味のある方は繰り返しに変換してください。

実行結果は次のようになります。 

>>> for x in range(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, end=' ')
        i += 1
        x = polygonal_number(5, i)
    print()

五角数を生成して、それが三角数になるか is_polygonal_number() でチェックするだけです。実行結果は次のようになります。 

>>> triangular_pentagonal_num(10000000)
1 210 40755 7906276

ところで、『オンライン整数列大辞典の数列 A069017』によると、三角数かつ五角数の漸化式は次のようになるそうです。

\(\begin{array}{l} f(0) = 0 \\ f(1) = 1 \\ f(n) = 194 \times f(n - 1) - f(n - 2) + 16 \end{array}\)

三角数かつ五角数の漸化式

これをプログラムすると次のようになります。 

リスト : 三角数で五角数の数 (2)

def triangular_pentagonal(n):
    if n < 2: return n
    return 194 * triangular_pentagonal(n - 1) - triangular_pentagonal(n - 2) + 16

二重再帰なので、実行速度は遅いです。興味のある方は繰り返しに変換してください。

実行結果は次のようになります。 

>>> for x in range(1, 9): print(triangular_pentagonal(x))
...
1
210
40755
7906276
1533776805
297544793910
57722156241751
11197800766105800

●解答５

Python でナイーブにプログラムすると次のようになります。 

リスト : 四角錐数

def square_pyramidal(n):
    a = 0
    for m in range(1, n + 1): a += m * m
    return a

>>> for x in range(1, 11): print(square_pyramidal(x))
...
1
5
14
30
55
91
140
204
285
385
>>> square_pyramidal(1000000)
333333833333500000

次の公式を使うともっと簡単に求めることができます。

四角錐数\((n) = \dfrac{n(n + 1)(2n + 1)}{6}\)

四角錐数の公式は次のように導出することができます。

\(\begin{array}{l} S = 1 + 2^{2} + 3^{2} + \cdots + n^{2} \\ U = 1 + 2^{3} + 3^{3} + \cdots + (n + 1)^{3} = 1 + (1 + 1)^{3} + (2 + 1)^{3} + \cdots + (n + 1)^{3} \end{array}\)

とする。\((a + b)^{3} = a^{3} + 3a^{2}b + 3ab^{2} + b^{3}\) なので、\(U\) は次のように変形できる。

\(\begin{array}{l} U = (1 + 2^{3} + \cdots + n^{3}) + 3(1 + 2^{2} + \cdots + n^{2}) + 3(1 + 2 + \cdots + n) + (n + 1) \\ U - (1 + 2^{3} + \cdots + n^{3}) = 3(1 + 2^{2} + \cdots + n^{2}) + 3(1 + 2 + \cdots + n) + (n + 1) \\ (n + 1)^{3} = 3S + \dfrac{3n(n + 1)}{2} + (n + 1) \\ 6S = 2(n + 1)^3 - 3n(n + 1) - 2(n + 1) \\ 6S = (n + 1)(2(n + 1)^2 - 3n - 2) \\ 6S = (n + 1)(2n^2 + 4n + 2 - 3n - 2) \\ 6S = (n + 1)(2n^2 + n) \\ S = \dfrac{n(n + 1)(2n + 1)}{6} \end{array}\)

●解答６

Python でナイーブにプログラムすると次のようになります。 

リスト : 三角錐数

def triangular_pyramidal(n):
    a = 0
    for m in range(1, n + 1): a += m * (m + 1) // 2
    return a

>>> for x in range(1, 11): print(triangular_pyramidal(x))
...
1
4
10
20
35
56
84
120
165
220
>>> triangular_pyramidal(1000000)
166667166667000000

次の公式を使うともっと簡単に求めることができます。

三角錐数\((n) = \dfrac{n(n + 1)(n + 2)}{6}\)

四角錐数の公式を使うと、三角錐数の公式は次のように求めることができます。

\(\begin{eqnarray} \displaystyle \sum_{k=1}^n \dfrac{k(k + 1)}{2} &=& \dfrac{\sum_{k=1}^n k^2 + \sum_{k=1}^n k}{2} \\ &=& \dfrac{n(n + 1)(2n + 1) + 3n(n + 1)}{12} \\ &=& \dfrac{2n^3 + 6n^2 + 4n}{12} \\ &=& \dfrac{n(n^2 + 3n + 2)}{6} \\ &=& \dfrac{n(n + 1)(n + 2)}{6} \end{eqnarray}\)

●参考 URL

  1. 多角数 - Wikipedia
  2. 三角数 - Wikipedia
  3. 平方三角数 - Wikipedia

●プログラムリスト

#
# polynum.py : 多角数
#
#              Copyright (C) 2015-2022 Makoto Hiroi
#
from math import sqrt

def polygonal_number(p, n):
    return ((p - 2) * n * n - (p - 4) * n) // 2

def test1():
    for p in range(3, 9):
        print(p, ":", end=' ')
        for x in range(1, 17):
            print(polygonal_number(p, x), end=' ')
        print()

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

def square_triangular_num(n):
    m = int(sqrt(n))
    for x in range(1, m + 1):
        if is_polygonal_number(3, x * x):
            print(x * x, end=' ')
    print()

def square_triangular(n):
    if n < 2: return n
    return 34 * square_triangular(n - 1) - square_triangular(n - 2) + 2

def triangular_pentagonal_num(n):
    i = 1
    x = 1
    while x <= n:
        if is_polygonal_number(3, x):
            print(x, end=' ')
        i += 1
        x = polygonal_number(5, i)
    print()

def triangular_pentagonal(n):
    if n < 2: return n
    return 194 * triangular_pentagonal(n - 1) - triangular_pentagonal(n - 2) + 16

def square_pyramidal(n):
    a = 0
    for m in range(1, n + 1): a += m * m
    return a

def triangular_pyramidal(n):
    a = 0
    for m in range(1, n + 1): a += m * (m + 1) // 2
    return a
初版 2015 年 12 月 20 日
改訂 2022 年 10 月 22 日