ピタゴラス数を求めよう！
a^2 + b^2 = c^2を満たす正の整数の組(a, b, c)のことをピタゴラス数といいます．
本記事では，ピタゴラス数を見つけたり，一定の範囲のピタゴラス数を数えたりして遊びます．
ピタゴラス数とは
a^2 + b^2 = c^2を満たす正の整数の組(a, b, c)のことをピタゴラス数いいます．
例えば，(3, 4, 5)は3^2 + 4^2 = 9 + 16 = 25 = 5^2 であるためピタゴラス数です．
また，(a, b, c)の最大公約数が1であるピタゴラス数のことを，原始ピタゴラス数といいます．
例えば，先ほどの(3, 4, 5)は原始ピタゴラス数ですし，(5, 12, 13)も原始ピタゴラス数です．逆に原始ピタゴラス数でないものの例として(6, 8, 10)や(9, 12, 15)などが挙げられます．ピタゴラス数を正の整数倍したものもピタゴラス数となりますので，原始ピタゴラス数でないピタゴラス数は，なんらかの原始ピタゴラス数を正の整数倍したものとなります．先ほどの(6, 8, 10)は(3, 4, 5)に2を掛けたものとなります．
ちなみに，a^2 + b^2 = c^2といえば，ピタゴラスの定理(三平方の定理)を連想しますね．ピタゴラスの定理が成り立つ直角三角形のうち，辺の長さが全て整数のものの組がピタゴラス数であると考えることもできます．そのように考えると，次のことが直ちに分かります．
命題．
a = bとなるピタゴラス数は存在しない．
ピタゴラス数を見つけよう！
さて，a^2 + b^2 = c^2を満たす正の整数にはどのようなものがあるか，探してみましょう！人力で探すのも面白いですが，ここではプログラムを書いてみましょう．
念のため探索する範囲を定めておきましょう．
問題．
正の整数の入力nに対して，a < b < c \leqq nを満たすピタゴラス数を出力するプログラムを作成せよ．
まずは愚直にfor文を回してみます．以下のコードの言語はPythonです．
# 正の整数nを入力
n = int(input('n:'))
# cを1からnまで１ずつ増やす．
for c in range(1,n+1):
# bを1からc-1まで1ずつ増やす．(b=cとなることは無い)
for b in range(1,c):
# aを1からb-1まで1ずつ増やす．(a=bとなることは無い)
for a in range(1,b):
# ピタゴラス数であれば出力
if a**2 + b**2 == c**2:
print((a, b, c))n=100のときの出力
n:100
(3, 4, 5)
(6, 8, 10)
(5, 12, 13)
(9, 12, 15)
(8, 15, 17)
(12, 16, 20)
(15, 20, 25)
(7, 24, 25)
(10, 24, 26)
(20, 21, 29)
(18, 24, 30)
(16, 30, 34)
(21, 28, 35)
(12, 35, 37)
(15, 36, 39)
(24, 32, 40)
(9, 40, 41)
(27, 36, 45)
(30, 40, 50)
(14, 48, 50)
(24, 45, 51)
(20, 48, 52)
(28, 45, 53)
(33, 44, 55)
(40, 42, 58)
(36, 48, 60)
(11, 60, 61)
(39, 52, 65)
(33, 56, 65)
(25, 60, 65)
(16, 63, 65)
(32, 60, 68)
(42, 56, 70)
(48, 55, 73)
(24, 70, 74)
(45, 60, 75)
(21, 72, 75)
(30, 72, 78)
(48, 64, 80)
(18, 80, 82)
(51, 68, 85)
(40, 75, 85)
(36, 77, 85)
(13, 84, 85)
(60, 63, 87)
(39, 80, 89)
(54, 72, 90)
(35, 84, 91)
(57, 76, 95)
(65, 72, 97)
(60, 80, 100)
(28, 96, 100)この結果を眺めているだけでも面白い発見があります．
ラマヌジャンが「1729」を，2通りの2つの立方数の和で表せる最小の数(1729 = 1^3 + 12^3 = 9^3 + 10^3)であると即答した逸話は有名です．これをスケールダウンさせて「2通りの2つの平方数の和で表せる最小の平方数」を考えてみると，プログラムの出力から，25^2 = 625であることが分かります(25^2 = 15^2 + 20^2 = 7^2 + 24^2)．
また，cの値が同一の原始ピタゴラス数の中で，cの値が最小のものが(33, 56, 65)と(16, 63, 65)であることが分かります．
眺めているだけでも発見があって楽しいですね！
ピタゴラス数を数えよう！
nを大きくした場合，全てのピタゴラス数を出力するのも大変ですので，少し問題を変えてみます．
問題．
正の整数の入力nに対して，a < b < c \leqq nを満たすピタゴラス数がいくつあるか出力するプログラムを作成せよ．
ピタゴラス数の出力から，ピタグラス数の個数(組数)の出力の問題に変えました．数学よりもプログラミングに寄った問題となっております．
こちらの問題もfor文で愚直に回してみましょう．
※高校生やプログラミング初心者を想定してコメントを書いています．コードの動きをただ説明しているだけのダメなコメントの書き方ですが，どうか温かい目で見てください．
# 正の整数nを入力
n = int(input('n:'))
# 変数countで個数を数える
count = 0
# cを1からnまで１ずつ増やす．
for c in range(1,n+1):
# bを1からc-1まで1ずつ増やす．(b=cとなることは無い)
for b in range(1,c):
# aを1からb-1まで1ずつ増やす．(a=bとなることは無い)
for a in range(1,b):
# ピタグラス数であれば個数をプラス1
if a**2 + b**2 == c**2:
count += 1
# 個数を出力
print(count)入出力例は以下の通りです．
n:100
52
さてここで大きな問題が発生します．
n = 1000のときなど，nの値が大きくなると，計算が終わりません．どうやら計算量が大きすぎるようです．nを大きくするために問題を変えたのに，これではいけませんね．
計算量を減らすために，さまざまな工夫の余地がありますが，ここでは数学的に解決してみましょう．
原始ピタゴラス数を見つける仕組み
原始ピタグラス数について以下の命題が成り立ちます．
命題．
正の整数の組(a, b, c)が原始ピタゴラス数であることと，
互いに素な正の整数m，nが m > nを満たし，m，nの一方が奇数でもう一方が偶数であり，
(a, b, c) = (2mn, m^2 - n^2, m^2 + n^2)
または，
(a, b, c) = (m^2 - n^2, 2mn, m^2 + n^2)
であることは同値である．
証明はしませんが，少々奇妙な条件である「一方が奇数でもう一方が偶数」の条件が無かった場合の事を少し考えてみましょう．
m，nが両方とも偶数の場合は，(m^2 - n^2, 2mn, m^2 + n^2)が全て偶数となり，(a, b, c)の最大公約数が2以上の数になってしまい，原始ピタゴラス数にはなりません．同様に，m，nが両方とも奇数の場合も，(m^2 - n^2, 2mn, m^2 + n^2)が全て偶数となり，原始ピタグラス数では無くなります．
さて，このような条件のm，nを探していくことで，原始ピタゴラス数を漏れなく見つけることができます．
改めてピタゴラス数を数えよう！
しかし上記の命題を使っても，原始ピタゴラス数しか見つかりません．ここからもう一工夫必要となります．ここで注目するのが，ピタゴラス数に正の整数を掛けたものもピタゴラス数になるということです．
では問題のおさらいと，計算量が小さくなるように改善したコードを見てみましょう．
問題．
正の整数の入力nに対して，a < b < c \leqq nを満たすピタゴラス数がいくつあるか出力するプログラムを作成せよ．
nを先ほどの命題に合わせて使いたいので，入力の変数を”num”にします．
# 最大公約数を求める関数
def gcd(a,b):
if b < a:
c = b
b = a
a = c
# n=0なら終了
if a == 0:
return b
# n=0ではない場合，m=nとn=(余り)で手順2に戻る
else:
rem = b % a
b = a
a = rem
return gcd(b, a)
# cの探索範囲を入力
num = int(input('num:'))
# ピタゴラス数の個数を記録するための変数count
count = 0
# m < √num
m_max = int(num**0.5)
# mを1からm_maxまで1ずつ増やす．
for m in range(1,m_max+1):
# nを1からn-1まで1ずつ増やす．
for n in range(1,m):
# mとn どちらかが偶数でもう一方が奇数である
if m%2 != n%2:
# m,n互いに素
if gcd(m,n) == 1:
# 定数倍したものもピタゴラス数なので一気に個数を追加
#（cが探索範囲を超えていた場合+0になる）
count += (num//(m**2+n**2))
print(count)まず最初に最大公約数を出力する関数を定義しています．
1行目から15行目までのコードです．ここでは，ユークリッドの互除法を使っています．
ユークリッドの互除法やそのコードについては，以下の記事に詳細があります．ユークリッドの互除法が分からない方はご参照ください．
22行目が計算量改善のポイントその1です．
m^2 + n^2がcの値となりますので，探索範囲の方のn，つまりプログラム中で”num”と表される値をxとすると，以下の式が成り立ちます．
m^2 + n^2 < x
m < \sqrt{x}
これによって探索範囲が一気に狭くなりました．
24行目以降が改善点その2です．上述の命題を活用しています．
そして最後の工夫にして，最重要なポイントが33行目です．
count += (num//(m**2+n**2))とても重要なことですが，今回の問題では，原始ピタゴラス数以外のピタゴラス数も個数に加算する必要があります．また，m^2 + n^2が探索範囲より大きくなってしまう可能性があります．その両方を一度に解決しているのがこの行です．Pythonの”//”は割り算の結果の整数部分のみを取り出し，小数点以下を切り捨てる演算子です．探索範囲内にその原始ピタゴラス数の倍数がいくつあるかを表すことができます．3倍したものまで探索範囲内であれば3になりますし，原始ピタゴラス数のみが範囲内にあれば1にまります．そして，m^2 + n^2が探索範囲より大きくなってしまった場合は，値が0になります．
さて最後にプログラムを実行した結果を見てみましょう．
num:1000
881num:10000
12471今回はどちらも一瞬で結果が出力されました．
これにて，プログラム完成です！
