Project Euler 91

プロジェクトオイラー
http://projecteuler.net/index.php

Q91.
原点と0≤x,y≤50に2点とって、直角三角形を作る。いくつ作れるか。

素直に点を打って、面積が正とある角が直角を条件に数え上げればよい。

from itertools import product
def is_positive_area(x1, y1, x2, y2):

    return x1 * y2 - x2 * y1 > 0
def is_right_angle(x1, y1, x2, y2):

    if x1 * x2 + y1 * y2 == 0:

        return True

    if x1 * (x2 - x1) + y1 * (y2 - y1) == 0:

        return True

    if x2 * (x2 - x1) + y2 * (y2 - y1) == 0:

        return True

    return False
N = 50
counter = 0

for x1, y1, x2, y2 in product(xrange(N + 1), repeat = 4):

    if is_positive_area(x1, y1, x2, y2):

        if is_right_angle(x1, y1, x2, y2):

            counter += 1

print counter

しかし、これだとO(N⁴)になってしまう。
原点をはさむ角が直角である場合はN²通り、原点から左回りに考えて、次の1点を決めてそこから直角に曲がる場合を数えて、それと対称形状が必ずあるから2倍する。これだと、O(N²)で済む（最大公約数の計算が入るから厳密には違うが）。

from itertools import product

import fractions
def calc(x1, y1):

    d = fractions.gcd(x1, y1)

    dx = -y1 / d

    dy =  x1 / d

    

    if dx == 0:

        return N

    else:

        return min(-x1 / dx, (N - y1) / dy)
N = 50
counter = N * N

for x1, y1 in product(xrange(1, N + 1), xrange(N)):

    counter += calc(x1, y1) * 2

print counter