Project Euler 228

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

Q228.
S_nを、頂点が
x_k = cos(π(2k - 1) / n)
y_k = sin(π(2k - 1) / n)
の正n角形とする。
ミンコフスキー和、S + Tを、それぞれの多角形の内部の点の位置ベクトル同士の和からなる多角形とする。S₁₈₆₄ + S₁₈₆₅ + ... + S₁₉₀₉は何角形か。

普通に組んだらやっぱり遅かったので、プログラミング上の工夫と数学的な工夫を一つずつ入れたら、なんとか答えが出た。非常に遅いし、まだ工夫する余地はあるが、もう気力がない。

from math import sin, cos, atan2, acos, pi, sqrt

from itertools import imap
def add(pt1, pt2):

    x = pt1[0] + pt2[0]

    y = pt1[1] + pt2[1]

    r = sqrt(x * x + y * y)

    theta = atan2(y, x)

    if theta < 0:

        theta += 2 * pi

    

    return (x, y, r, theta)
def dist_theta(pt1, pt2):

    d = abs(pt1[3] - pt2[3])

    if d > pi:

        d = 2 * pi - d

    return d
def get_nearest_point_index(pt0, P):

    i = int(pt0[3] * len(P) / (2 * pi))

    pt = P[i]

    min_d = dist_theta(pt, pt0)

    if pt[3] < pt0[3]:

        while True:

            i = next_index(i, P)

            d = dist_theta(P[i], pt0)

            if d < min_d:

                min_d = d

            else:

                return prev_index(i, P)

    else:

        while True:

            i = prev_index(i, P)

            d = dist_theta(P[i], pt0)

            if d < min_d:

                min_d = d

            else:

                return next_index(i, P)

    return i
def get_nearest_point(pt0, P):

    return P[get_nearest_point_index(pt0, P)]
def prev_index(i, P):

    if i == 0:

        return len(P) - 1

    else:

        return i - 1
def next_index(i, P):

    n = len(P)

    if i == n - 1:

        return 0

    else:

        return i + 1
def is_left(pt, pt1, pt2):

    eps = 1e-12

    x, y, x1, y1, x2, y2 = pt[0], pt[1], pt1[0], pt1[1], pt2[0], pt2[1]

    return (y2 - y1) * (x - x1) < (x2 - x1) * (y - y1) + eps
def add_point_core(pt, i, P):

    if i == 0 and pt[3] < 0:

        P.append(pt)

        return len(P) - 1

    else:

        P.insert(i, pt)

        return i
def add_point(pt, P):

    i1 = get_nearest_point_index(pt, P)

    if pt[3] < P[i1][3]:

        i2 = i1

        i1 = prev_index(i2, P)

    else:

        i2 = next_index(i1, P)

    

    pt1 = P[i1]

    pt2 = P[i2]

    

    if is_left(pt, pt1, pt2):

        return

    

    i2 = add_point_core(pt, i2, P)

    

    while True:

        i3 = next_index(i2, P)

        i4 = next_index(i3, P)

        if not is_left(P[i3], pt, P[i4]):

            break

        P.pop(i3)

        if i3 < i2:

            i2 -= 1

    

    while True:

        i3 = prev_index(i2, P)

        i4 = prev_index(i3, P)

        if not is_left(P[i3], P[i4], pt):

            break

        P.pop(i3)

        if i3 < i2:

            i2 -= 1
def print_point(pt):

    print "(%8.3f, %8.3f, %8.3f, %8.3f)" % (pt[:3] + (pt[3] * 180 / pi, ))
def print_polygon(P):

    print "%d points :" % len(P)

    for pt in P:

        print_point(pt)
def calc_segment_length(pt1, pt2):

    dx = pt2[0] - pt1[0]

    dy = pt2[1] - pt1[1]

    a = dx * dx + dy * dy

    b = pt1[0] * dx + pt1[1] * dy

    c = pt1[2] * pt1[2]

    t = -b / a

    if t < 0:

        return pt1[2]

    elif t > 1:

        return pt2[2]

    else:

        return (pt1[0] * dy - pt1[1] * dx) / sqrt(a)
def sum_polygon(P1, P2):

    P = [ ]

    for pt1 in P1:

        pt2 = get_nearest_point(pt1, P2)

        P.append(add(pt1, pt2))

    

    R = calc_segment_length(P[-1], P[0])

    R = reduce(min, imap(lambda k: calc_segment_length(P[k], P[k+1]),

                                    xrange(len(P) - 1)), R)

    n = len(P2)

    for pt1 in P1:

        r0 = pt1[2]

        th0 = pt1[3]

        phi = acos( (R * R - r0 * r0 - 1) / (2 * r0) )

        k0 = (n * (th0 - phi) / pi - 1) / 2 - 1e-6

        k1 = (n * (th0 + phi) / pi - 1) / 2 + 1e-6

        k0 = int(k0 + n + 1) % n

        k1 = int(k1) % n

        k = k0

        while True:

            add_point(add(pt1, P2[k]), P)

            if k == k1:

                break

            k = next_index(k, P2)

    

    return P
def create_point(k, n):

    theta = pi * (2 * k - 1) / n;

    return ( cos(theta), sin(theta), 1.0, theta )
def create_polygon(n):

    return [ create_point(k, n) for k in range(1, n + 1) ]
M = 1864

N = 1909

print len(reduce(sum_polygon, map(create_polygon, range(M, N + 1))))