Project Euler 91（1）

Problem 91
点P(x₁, y₁)とQ(x₁, y₂)が格子点上にあり、これらに原点O(0, 0)を加えて△OPQを形成する。
どの座標も0と2の間にある場合、14個の直角三角形がある（中略）。
0 ≤ x_1, y_1, x_2, y_2 ≤ 50のとき、直角三角形はいくつあるか。
http://projecteuler.net/index.php?section=problems&id=91

L=50として、大きさLの正方形内のPとQを単純に生成し（ただし、PはOQの右下）、直角三角形か判定するだけで答えは出ます（Code1）。しかし、これではあまりに雑ですね。計算量もO(L⁴)です。もうちょっとなんとかしましょう。

Oが直角の三角形はL²個あります。ここから、Pが直角だとして考えます。Pが軸上にあるとき、三角形は2L²個あります。そうでないとき、例えばP(4, 2)とすると、OPに垂直な方向は(1, -2)なので、この方向の格子点にQを取れます。左上には4/1 = 4個、右下には2/2 = 1個取れます。このようにすれば計算量がO(L²logL)（？）になって劇的に速くなります。

もうちょっと速くなると思うのですが。

# Code1
from itertools import ifilter, product

def is_right_triangle(P, Q):
    if P[1] == 0 and Q[0] == 0:
        return True
    elif P[0] * (Q[0] - P[0]) == P[1] * (P[1] - Q[1]):
        return True
    elif Q[0] * (P[0] - Q[0]) == Q[1] * (Q[1] - P[1]):
        return True
    return False

def is_right_lower(P, Q):
    return P[0] * Q[1] - P[1] * Q[0] > 0

L = 50
counter = 0
for P in product(range(L + 1), repeat = 2):
    for Q in product(range(L + 1), repeat = 2):
        if is_right_lower(P, Q) and is_right_triangle(P, Q):
            counter += 1
print counter

# Code2
from itertools import ifilter, product
from fractions import gcd

L = 1000
counter = L * L * 3
for x, y in product(range(1, L + 1), repeat = 2):
    d = gcd(x, y)
    dx = y / d
    dy = x / d
    counter += min((L - x) / dx, y / dy) + min(x / dx, (L - y) / dy)
print counter