原始ピタゴラス数を高速に求める
原始ピタゴラス数生成公式
$${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)