- 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。
問題 91.頂点が整数座標の直角三角形
原文 Problem 91: Right triangles with integer coordinates
問題の要約: 2点$P(x_1,y_1),Q(x_1,y_1),$と原点$O(0,0)$で作られる三角形OPQが直角三角形となるのは $0 \le x_1, x_2, y_1,y_2 \le 50$のとき何個あるか求めよ
$0 \le x_1, x_2, y_1,y_2 \le 2$の時は以下の14個になる。
全探索による解 (Ver 1.0)
まず全探索ですべてのPQの組み合わせで直角三角形になっているかチェックしたのが以下のプログラム。Problem 100番以下であればこのやり方でもOKなのかもしれませんが、今後のこともあるのでもう少し掘り下げてみたいと思います。
# Ver 1.0
from scipy.spatial import distance
from itertools import product,combinations
def lineLen(p):
return distance.euclidean(p[0],p[1])
def isRightTr(p3): # return True if p3 is right triangle
e3 = sorted([lineLen(p2) for p2 in combinations(p3,2) ])
return abs(e3[0]**2+e3[1]**2-e3[2]**2)<0.000001
r = 50
pq = [[i,j] for (i,j) in product(range(r+1),repeat=2) if i+j != 0 ] # all point except (0,0)
print(f"Answer : {sum([isRightTr([p,q,[0,0]]) for p, q in combinations(pq,2)])}")
直角三角形を探索する別解（Ver 2.0）
まず直角三角形を、以下の2つに分類して、各々の数を数えることを考えます。
- 2辺がそれぞれx軸、y軸と平行なもの
- それ以外
1. 2辺がそれぞれx軸、y軸と平行なもの
例題の14個の直角三角形の内赤丸をつけた12個です。
これは原点を頂点とする長方形１つに対して3個あることが分かります。上の例では左下の正方形の中に１～３の数字を書いた3個です。したがって
原点を頂点とする長方形の数はは \\
0 \le x_1, x_2, y_1,y_2 \le r のときr^2個なので\\
直角三角形の数は3r^2となります。
2. それ以外の直角三角形
それ以外は例題の$r=2$では2個しかないので簡単かなと思ったら全然違いました。かなり複雑ですが以下のようなステップで考えて行きます。
- 原点から領域内のすべての格子点へ線を引く（赤線）
- そこから左右に垂線を引いて格子点にあたる点を求める（青の点）
- その格子点から原点もどる線を引くと直角三角形が出来る。（青線）
これを$r=5$で傾きベクトル$(2,1)$の線を引いた場合で考えると以下のように4つの直角三角形が見つかります。
これをもれなく探すアルゴリズムを求めて行きます。
領域内のすべての格子点へ原点から線を引く
まともに引くと$r^2$本必要ですが以下のように工夫すると減らすことが出来ます。
- x軸から45度以内を考えれば45度～90度は対称なので2倍すればよい。（ただし45度線上は一つだけ）
- 上の赤線の傾きベクトル$(2,1)$の例で分かるように$(2,1)$から$(4,2)$、$(6,3)$と伸ばしていける。
- この時 $+90,-90度$回転させたベクトルは、それぞれに対して$(-1,2)$、$(1,-2)$となりこの回転ベクトルを加えると必ず格子点になる。
従って見つける必要があるのは原点から直線を引いたときに最初に当たる点。これはすなわちユークリッドの果樹園(Wikipedia)ということになるので$gcd(x,y)=1$となる点を探すことと同値になります。
ここでまたひと工夫。$gcd(x,y)=1$の点を探すのはやはり$O(n^2)$となってしまうので。見方を変えて$0 < \frac{y}{x} < 1$の既約分数を探す問題に置き換えることが出来ます。
するとProblem 73: 範囲内の分数の数の別解として紹介したStern-Brocot treeを使えば効率よく生成することが出来ます。以下のコードでは分母が4以下の分数を生成してx,yを入れ替えると最初に当たる点がリストされます。
### Stern–Brocot tree
def sternBrocot(nl, dl, nr, dr, Dmax):
nm, dm = nl+nr, dl+dr
if dm > Dmax: return []
return sternBrocot(nl, dl, nm, dm, Dmax)+sternBrocot(nm, dm, nr, dr, Dmax)+[(nm, dm)]
rpf = sternBrocot(0,1,1,1, 4) # 0/1 < n/d < 1/1, Dmax = 4
print(f"{len(rpf)}, {[pf[::-1] for pf in rpf]}")
# 5, [(4, 1), (3, 1), (4, 3), (3, 2), (2, 1)]
アルゴリズムのまとめ
以上をまとめると格子のサイズが$r$の時のアルゴリズムは以下のようなステップになります。
- Stern-Brocot treeを使って分母が$r-1$以下の$0 < \frac{y}{x} < 1$の既約真分数をすべて列挙する。
- その分数の$(x,y)$をベクトルとしてを原点から次々と伸ばした点の各々で、左右90度にベクトルをまた次々に伸ばして格子の範囲ならカウントする。
この動きを$r=12$でベクトル$v=(2,1),(3,1), (4,1),(3,2)$に関して図示したのが以下の図です。countがそこで見つかった直角三角形の数を表しています。
このステップをコードにしたのが以下になります。注意点としては、
- Stern-Brocot treeの出力には１は含まれないのでベクトル(1,1)を追加する
- ベクトル(1,1)は左右対称なので右90だけ見る
- x軸から45度以内を求めた数を2倍する
- 最後に2辺がそれぞれx軸、y軸と平行なもの数$3r^2$を加える
# Ver 2.0
import itertools
import numpy as np
def inR(r,pt): #check within 0-r range
return (np.all(pt<=r) and np.all(pt>=0))
def c90Tri(r, pf, pfi, rot): # rotate left/right rotate
pfrot = np.dot(rot,pf)
for j in itertools.count(1):
if not inR(r, pfi + pfrot*j): break #check within 0-r range
return j-1
l90,r90 = [[0,-1],[1,0]], [[0,1],[-1,0]]
def cTri(r,pf):
cnt = 0
for pfi in itertools.count(pf, pf): #pfi: pf x i
if not inR(r,pfi): break # if out of range
cnt += c90Tri(r, pf, pfi, l90) # left 90 degree rotate
if np.all(pf == (1,1)): continue # for vector (1,1) only check left90 (symmetric)
cnt += c90Tri(r, pf, pfi, r90) # right 90 degree rotate
return cnt
def countTri(r, rpf):
return sum([cTri(r,np.array(pf[::-1])) for pf in rpf])*2
r = 50
# list up all proper fractions: max denominator = r-1
rpf = sternBrocot(0,1,1,1, r-1) # 0/1 < n/d < 1/1, Dmax = r-1
print(f"Answer: {3*(r**2)+countTri(r, np.array(rpf+[(1,1)]))}")
パフォーマンス測定
Ver 1.0 と 2.0で実行時間を測定した結果です。かなりの改善が行われていることが分かります。
| 50 | 100 | 500 | 1000 | |
|---|---|---|---|---|
| Ver 1.0 | 155s | - | - | - |
| Ver 2.0 | 0s | 0s | 17s | 72s |
（開発環境：Google Colab）