原始ピタゴラス数を高速に求める

原始ピタゴラス数生成公式

$${m, n}$$ を次の 3 条件を満たす整数とする。  
(1) $${m>n>0}$$
(2) $${m−n}$$ は奇数
(3) $${gcd(m,n)=1}$$
このとき,$${(a,b,c) = (m^2−n^2,2mn,m^2+n^2)}$$ は,原始ピタゴラス数である。逆に,すべての原始ピタゴラス数は上の形で表すことができる。

原始ピタゴラス数生成公式を利用した python コード

# a, b の最大公約数(互除法利用)
def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

# 原始ピタゴラス数
def primitive_pythagorean_triple(k):
    ppt_list = []  # 記録するために空のListを用意
    m = 2  # (1) より m >= 2
    while m*m < k:  # c = m*m + n*n <= k より m*m <= k - n*n < k
        for n in range(1, m):  # (1) より
            if (m-n) % 2 == 1 and gcd(m, n) == 1:  # (2), (3) より
                c = m*m + n*n
                if c > k:
                    break

                a = m*m - n*n
                b = 2 * m * n
                if a < b:
                    ppt_list.append((a, b, c))  # 求まった値の組を記録
                else:
                    ppt_list.append((b, a, c))  # 求まった値の組を記録

        m += 1
    return sorted(ppt_list, key = lambda x: (x[2], x[1]))  # c, b で昇順にソート


if __name__ == '__main__':
    upper_limit = 1000  # 求める上限値

    ppt_list = primitive_pythagorean_triple(upper_limit)  # 原始ピタゴラス数のリスト
    
    for i, ppt in enumerate(ppt_list, 1):
        print(i, ' ', *ppt)

単純に求めるなら(低速)

三重ループのコード

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def primitive_pythagorean_triple(k):
    ppt_list = []
    for c in range(1, k + 1, 2):  # c は奇数なので、増分を2にする
        for b in range(1, c):  # b < c なので
            if gcd(c, b) != 1:  # b, c の順だと計算が1回多くなるので
                continue  # b, c が互いに素でないので、次の b へ

            a2 = c*c - b*b
            for a in range(c - b + 1, b):  # |b - c| < a < b なので
                if a*a == a2:
                    ppt_list.append((a, b, c))
                    break  # a が見つかったので、次の b へ
                elif a*a > a2:
                    break  # これ以上探しても無駄なので、次の b へ
                    
    return ppt_list


upper_limit = 1000
ppt_list = primitive_pythagorean_triple(upper_limit)
for i, ppt in enumerate(ppt_list, 1):
    print(i, ' ', *ppt)

二重ループ + 二分探索のコード

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

# 二分探索で a を探す
def binary_search(b, c):
    a2 = c*c - b*b
    left, right = c - b + 1, b - 1  # a を探す範囲
    while left <= right:  # 探す範囲を left = right になるまで狭めていく
        mid = (left + right) // 2  # 区間の中央値
        if mid*mid == a2:  # 見つかった。mid が a
            return mid
        elif mid*mid < a2:  # [left, mid] には a は存在しない
            left = mid + 1  # 探す範囲を [mid + 1, right] にする
        else:  # [mid, right] には a は存在しない
            right = mid - 1  # 探す範囲を [left, mid - 1] にする

    return False  # 見つからなかった

def primitive_pythagorean_triple(k):
    ppt_list = []
    for c in range(1, k + 1, 2):
        for b in range(1, c):
            a = binary_search(b, c)
            if a and gcd(b, a) == 1:
                ppt_list.append((a, b, c))
                
    return ppt_list


upper_limit = 1000
ppt_list = primitive_pythagorean_triple(upper_limit)
for i, ppt in enumerate(ppt_list, 1):
    print(i, ' ', *ppt)

二重ループ + is_integer() を使ったコード

float型の数値が整数(小数点以下が0)か判定する is_integer() メソッドを利用

from math import sqrt, gcd

def primitive_pythagorean_triple(k):
    ppt_list = []
    for c in range(1, k + 1, 2):
        for b in range(1, c):
            a = sqrt(c*c - b*b)  # 例えば sqrt(5^2 - 4^2) = 3.0 で float型
            if a.is_integer():  # 小数点以下が 0 か判定
                a = int(a)  # float型のままでは不都合なのでint型へ
                if a < b and gcd(b, a) == 1:
                    ppt_list.append((a, b, c))
                    
    return ppt_list


upper_limit = 1000
ppt_list = primitive_pythagorean_triple(upper_limit)
for i, ppt in enumerate(ppt_list, 1):
    print(i, ' ', *ppt)

この記事が気に入ったらサポートをしてみませんか?