Project Euler 246 - inamori’s diary

Q246.
与えられた楕円の外の格子点で、そこから楕円に向かって引いた2本の接線が成す角度が45度より大きいものはいくつあるか。

これも前問と並んで最難問題となっているが、実際には易しい。1点ずつ調べても解ける。45度になる点の軌跡は解析的にはわからなかったが、それがわかる必要もない。対称性から、楕円の中心を原点とする座標系で第1象限だけ調べればよいが、軸上の点はすぐにわかる。x軸上から一つ上がって境界を探す、というのを繰り返せば十分に速い。速度を気にするコーディングをしなくてもいいレベル。



from math import *

from itertools import count
def calc_tangential_lines_angle(x0, y0):

    A = (a * y0) ** 2 + (b * x0) ** 2

    B = (a * b) ** 2 * x0

    C = a ** 4 * (b * b - y0 * y0)

    D = max(0.0, B * B - A * C)

    det_root = sqrt(D)

    x1 = (B + det_root) / A

    x2 = (B - det_root) / A

    y1 = b * b / y0 * (1 - x1 * x0 / (a * a))

    y2 = b * b / y0 * (1 - x2 * x0 / (a * a))

    

    dx1 = x1 - x0

    dy1 = y1 - y0

    dx2 = x2 - x0

    dy2 = y2 - y0

    return acos( (dx1 * dx2 + dy1 * dy2) / sqrt(

            (dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2)))
def is_in_region(x, y):

    return calc_tangential_lines_angle(x, y) > limit_angle
def count_points():

    counter = 0

    x_bound = int(x_max)

    for y0 in xrange(1, int(y_max) + 1):

        if y0 <= b:

            x_start = int(a * sqrt(1 - y0 * y0 / (b * b))) + 1

        else:

            x_start = 1

        

        if is_in_region(x_bound, y0):

            x_bound += 1

            while is_in_region(x_bound, y0):

                x_bound += 1

            x_bound -= 1

        else:

            x_bound -= 1

            while not is_in_region(x_bound, y0):

                x_bound -= 1

        

        counter += x_bound - x_start + 1

    return counter
M = (-2000.0, 1500.0)

G = ( 8000.0, 1500.0)

R = 15000.0

limit_angle = pi / 4
L = G[0] - M[0]

a = R / 2

b = sqrt(R ** 2 - L ** 2) / 2
th1 = atan(b / a * tan(pi / 8))

th2 = atan(b / a * tan(3 * pi / 8))

y_max = b / sin(th1)

x_max = a / cos(th2)
print count_points() * 4 + (int(y_max) - int(b) + int(x_max) - int(a)) * 2