Puzzle DE Programming

パズルの説明

ピタゴラスの定理は、平面幾何学において直角三角形の辺を \(a, b, c \ (a + b \gt c)\) とすると、次式が成り立つという皆さんお馴染みの有名な定理です。

「ピタゴラス数」または「ピタゴラスの三つ組数 (pythagoras triple)」は、上式を満たす自然数の組 (a, b, c) のことで、とくに a, b, c が互いに素であるとき、(a, b, c) を「原始ピタゴラス数」といいます。たとえば、(3, 4, 5), (5, 12, 13), (8, 15, 17) などがあります。

参考 URL 『ピタゴラスの定理 - Wikipedia』によると、原始ピタゴラス数は次の方法で簡単に見つけることができるそうです。

• (\(a, b, c\)) = (\(m^2 - n^2, 2mn, m^2 + n^2\)) or (\(2mn, m^2 - n^2, m^2 + n^2\))
• ただし、\(m \gt n\), \(m\) と \(n\) は互いに素, \(m - n\) が奇数 (\(m + n\) が奇数でもよい)

本ページでは証明を割愛しますが、「ピタゴラス数 - 桃の天然水」の説明がわかりやくて参考になると思います。inamori さんに感謝いたします。

それでは問題です。

1. 三辺の合計値 (a + b + c) が k 以下の原始ピタゴラス数をすべて求めるプログラムを作ってください。
2. 三辺の合計値 (a + b + c) が k 以下のピタゴラス数をすべて求めるプログラムを作ってください。

本ページでは Python3 (ver 3.8.10) を使ってプログラムを作ることにします。

●解答１

プログラムは次のようになります。

リスト : 三辺の合計値が k 以下の原始ピタゴラス数を求める

def primitive_pythagoras(k):
    xs = []
    m = 2
    while m * m <= k // 2:
        for n in range(1, m):
            if 2 * m * (m + n) > k: break
            if (m - n) % 2 != 0 and gcd(m, n) == 1:
                a = m * m - n * n
                b = 2 * m * n
                c = m * m + n * n
                if a < b:
                    xs.append((a, b, c))
                else:
                    xs.append((b, a, c))
        m += 1
    return xs

三辺の合計値は次の式で求めることができます。

上式から三辺の合計値は必ず偶数になることがわかります。変数 \(m\) の値を \(\sqrt {\frac{k}{2}}\) とすると、\(n = 1\) のときに三辺の合計値は \(2m^2 + 2m = k + 2 \sqrt {\frac{k}{2}}\) となり、\(k\) よりも大きくなるので、\(m\) の上限値を \(\sqrt {\frac{k}{2}}\) とします。

while ループで m の値を 2 から \(\sqrt {\frac{k}{2}}\) まで増やしていき、for ループの中で原始ピタゴラス数を生成します。三辺の合計値が k を超えたならば break で for ループを脱出します。そうでなければ、m - n が奇数で、m と n が互いに素であることを確認します。そして、辺 a, b, c を計算して append() でリスト xs にタプル (a, b, c) を追加するだけです。

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

>>> from pythagoras import *
>>> primitive_pythagoras(100)
[(3, 4, 5), (5, 12, 13), (8, 15, 17), (7, 24, 25), (20, 21, 29), (9, 40, 41),
 (12, 35, 37)]
>>> primitive_pythagoras(500)
[(3, 4, 5), (5, 12, 13), (8, 15, 17), (7, 24, 25), (20, 21, 29), (9, 40, 41),
 (12, 35, 37), (11, 60, 61), (28, 45, 53), (33, 56, 65), (13, 84, 85), (16, 63, 65),
 (48, 55, 73), (39, 80, 89), (15, 112, 113), (36, 77, 85), (65, 72, 97), (17, 144, 145),
 (20, 99, 101), (60, 91, 109), (51, 140, 149), (19, 180, 181), (44, 117, 125),
 (88, 105, 137), (85, 132, 157), (57, 176, 185), (21, 220, 221), (24, 143, 145),
 (119, 120, 169), (95, 168, 193), (52, 165, 173), (104, 153, 185), (133, 156, 205),
 (28, 195, 197), (84, 187, 205)]

●解答２

三辺の合計値が k となる全てのピタゴラス数を求めることも簡単です。つきのリストを見てください。

リスト : 合計値が k となるピタゴラス数をすべて求める

def pythagoras(k):
    xs = []
    for a, b, c in primitive_pythagoras(k):
        p = k % (a + b + c)
        if p == 0:
            q = k // (a + b + c)
            xs.append((q * a, q * b, q * c))
    return xs

primitive_pythagoras(k) で原始ピタゴラス数を生成し、k が三辺の合計値 (a + b + c) で割り切れることを確認するだけです。とても簡単ですね。

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

>>> pythagoras(12)
[(3, 4, 5)]
>>> pythagoras(120)
[(30, 40, 50), (20, 48, 52), (24, 45, 51)]
>>> pythagoras(240)
[(60, 80, 100), (40, 96, 104), (48, 90, 102), (15, 112, 113)]

もうひとつ、クールな方法を紹介します。この方法は nineties さんのブログ「ブートストラッピングでコンパイラを作る日記」の Project Euler (Problem 3～10) を参考にさせていただきました。nineties さんに感謝いたします。

式 \(a^2 + b^2 = c^2\) と \(a + b + c = k\) を使って変数 c を削除して因数分解すると、次のようになります。

上式より \(k - a\) と \(k - b\) は \(\frac{k^2}{2}\) の約数であることがわかります。つまり、\(\frac{k^2}{2}\) の約数を求め、条件 \(a \lt b \lt k\) を満たすものを探せばいいわけです。プログラムは次のようになります。

リスト : 合計値が k となるピタゴラス数をすべて求める (2)

def pythagoras1(k):
    xs = divisor(k * k // 2)
    ys = []
    i = 0
    j = len(xs) - 1
    while i <= j:
        if xs[j] < k:
            a = k - xs[j]
            b = k - xs[i]
            ys.append((a, b, k - (a + b)))
        i += 1
        j -= 1
    return ys

def test(func, k):
    s = time.time()
    print(func(k))
    print(time.time() - s)

拙作のページ「約数」で作成した関数 divisor() は約数をリストに昇順に格納して返します。約数の一つを xs[i] (i = 0, 1, 2, ...) とし、もう一つの約数を xs[j] (j = len(xs) - 1, len(xs) - 2, ...) とします。i が j よりも大きくなったら、すべての組み合わせを調べたので while ループを終了します。xs[j] が k よりも小さい場合は条件を満たしているので、a, b, c を計算してリスト ys に追加します。

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

>>> test(pythagoras, 10000000)
[(2000000, 3750000, 4250000), (2187500, 3600000, 4212500), (234375, 4880000, 4885625)]
1.2602739334106445
>>> test(pythagoras1, 10000000)
[(234375, 4880000, 4885625), (2000000, 3750000, 4250000), (2187500, 3600000, 4212500)]
0.0004527568817138672

実行環境 : Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz

k が大きな値の場合、約数を高速に求めることができれば pythagoras1() のほうが速くなるようです。興味のある方はいろいろ試してみてください。

●参考 URL

1. ピタゴラスの定理 - Wikipedia
2. ピタゴラス数 - 桃の天然水, (inamori さん)
3. Project Euler (Problem 3～10) - ブートストラッピングでコンパイラを作る日記, (nineties さん)

●プログラムリスト

#
# pythagoras.py : ピタゴラス数
#
#                 Copyright (C) 2016-2022 Makoto Hiroi
#
import time

# 素因数分解
def factor_sub(n, m):
    c = 0
    while n % m == 0:
        c += 1
        n //= m
    return c, n

def factorization(n):
    buff = []
    c, m = factor_sub(n, 2)
    if c > 0: buff.append((2, c))
    c, m = factor_sub(m, 3)
    if c > 0: buff.append((3, c))
    x = 5
    while m >= x * x:
        c, m = factor_sub(m, x)
        if c > 0: buff.append((x, c))
        if x % 6 == 5:
            x += 2
        else:
            x += 4
    if m > 1: buff.append((m, 1))
    return buff

# p^q の約数を求める
def divisor_sub(p, q):
    a = []
    for i in range(0, q + 1):
        a.append(p ** i)
    return a

# 約数を求める
def divisor(n):
    xs = factorization(n)
    ys = divisor_sub(xs[0][0], xs[0][1])
    for p, q in xs[1:]:
        ys = [x * y for x in divisor_sub(p, q) for y in ys]
    return sorted(ys)

# 最大公約数
def gcd(a, b):
    while b > 0:
        a, b = b, a % b
    return a

# 原始ピタゴラス数を求める
def primitive_pythagoras(k):
    xs = []
    m = 2
    while m * m <= k // 2:
        for n in range(1, m):
            if 2 * m * (m + n) > k: break
            if (m - n) % 2 != 0 and gcd(m, n) == 1:
                a = m * m - n * n
                b = 2 * m * n
                c = m * m + n * n
                if a < b:
                    xs.append((a, b, c))
                else:
                    xs.append((b, a, c))
        m += 1
    return xs

# ピタゴラス数を求める
def pythagoras(k):
    xs = []
    for a, b, c in primitive_pythagoras(k):
        p = k % (a + b + c)
        if p == 0:
            q = k // (a + b + c)
            xs.append((q * a, q * b, q * c))
    return xs

# 別解
def pythagoras1(k):
    xs = divisor(k * k // 2)
    ys = []
    i = 0
    j = len(xs) - 1
    while i <= j:
        if xs[j] < k:
            a = k - xs[j]
            b = k - xs[i]
            ys.append((a, b, k - (a + b)))
        i += 1
        j -= 1
    return ys

def test(func, k):
    s = time.time()
    print(func(k))
    print(time.time() - s)