Problem:


ハッピーエンド問題

「ハッピーエンド問題」と呼ばれる数学の未解決問題に関連したプログラムを書いてみましょう。平面上に与えられたN個の点から、ちょうどk個の点を結んでできる凸多角形のうち、最も面積の小さいものを見つけるプログラムを作成してください。ただし、N個の点の座標を与えられた後、質問として凸多角形の角の個数kがいくつか与えられます。

（補足：ハッピーエンド問題について）

平面上にどの３点も同じ直線上に乗らないようにN個の点を置きます。そのとき、どのように点を置いても、k個の点をうまく選ぶとk個の角をもつ凸多角形が必ず作れると予想されています。 今のところ、三角形ならばN=3、四角形ならばN=5、五角形ならばN=9、六角形ならばN=17であればよいことが、２００６年までに証明されています。また、三角形以上のすべてのk角形に対して、N=1+2k-2という予想がありますが、いまだ証明されていません。これは１００年にわたり研究が進められている難問です。

この問題には、「ハッピーエンド問題」という素敵な名前がつけられています。ある男女がこの問題を研究しているうちに仲良くなって、ついに結婚したことにちなんで、友人の数学者が名付けたそうです。ロマンチックですね。

入力

入力は1つのデータセットからなる。入力データは以下の形式で与えられる。

N
x1 y1
:
xN yN
Q
k1
:
kQ

１行目に平面上の点の個数N(3 ≤ N ≤ 40)が与えられる。続くN行に各点の座標が与えられる。各点には1からNまでの番号が、入力される順番に付けられている。xiとyi(-10000 ≤ xi, yi ≤ 10000)はi番目の点のそれぞれx座標とy座標を表す整数である。x軸の正方向は右向き, y軸の正方向は上向きに取るものとする。

続く１行に質問の数Q(1 ≤ Q ≤ N)が与えられる。続くQ行に質問が与えられる。ki(3 ≤ ki ≤ N)はi番目の質問である凸多角形の角の個数を表す。

なお、入力は以下の条件を満たすものとする。

 入力される点の座標はすべて異なる。 

 どの３点も同じ直線上には乗らない。

 各質問に対して面積最小の凸多角形は１つであり、２番目に小さい凸多角形との面積の差は 0.0001以上。

出力

質問ごとに、面積が最小の凸多角形の全頂点を1行に出力する。凸多角形の頂点で最も下にあるものの中で最も左にある頂点から順に、反時計周りで頂点の番号を出力する。頂点の番号の間は空白1つで区切る。行の終わりには空白文字を出力しない。ただし、凸多角形が作れない場合はNAと出力する。

入力例 1

5
0 0
3 0
5 2
1 4
0 3
3
3
4
5

出力例 1

1 4 5
1 2 4 5
1 2 3 4 5

入力例 2

6
0 0
3 0
5 2
1 4
0 3
3 2
4
3
4
5
6

出力例 2

6 3 5
6 3 4 5
1 2 6 4 5
NA



Sample input:
5
0 0
3 0
5 2
1 4
0 3
3
3
4
5


Sample output:
1 4 5
1 2 4 5
1 2 3 4 5

Submission status: Accepted

Solution code:
def cross(z1, z2):
    return z1.real * z2.imag - z1.imag * z2.real

def ccw(p1, p2, p3):
    return cross(p2 - p1, p3 - p1) > 0

def triangle_area(p1, p2, p3):
    # returns signed trangle area
    return cross(p2 - p1, p3 - p1) / 2

from sys import stdin
file_input = stdin

N = int(file_input.readline())
P = []
M = {}

for i in range(1, N + 1):
    x, y = map(int, file_input.readline().split())
    z = x + y * 1j
    P.append(z)
    M[z] = i # point to id

Q = int(file_input.readline())
query = [int(file_input.readline()) for i in range(Q)]
ql = max(query) - 2

P.sort(key=lambda x: (x.imag, x.real))

max_area = 20000 * 20000
CP_area = [max_area for i in range(ql)]
CP_points = [None for i in range(ql)]

from cmath import phase
from itertools import combinations

for i, bp in enumerate(P[:-2], start=1):
    tP = P[i:]
    tP.sort(key=lambda x: phase(x - bp))
    pn = N - i
    
    tM = dict(zip(range(pn), tP)) # id to point
    
    # making triangles
    t_rec = [[None for k in range(pn)] for j in range(pn)]
    for j, k in combinations(range(pn), 2):
        pj, pk = tM[j], tM[k]
        ta = triangle_area(bp, pj, pk)
        t_rec[j][k] = [ta, [bp, pj, pk]]
        if ta < CP_area[0]:
            CP_area[0] = ta
            CP_points[0] = [bp, pj, pk]
    
    # making convex polygons
    prev = t_rec.copy()
    cur = [[None for k in range(pn)] for j in range(pn)]
    for i in range(1, ql):
        for j, k, s in combinations(range(pn), 3):
            pre_cp = prev[j][k]
            if pre_cp:
                pj, pk, ps = tM[j], tM[k], tM[s]
                if ccw(pj, pk, ps):
                    ta = pre_cp[0] + t_rec[k][s][0]
                    if not cur[k][s] or cur[k][s][0] > ta:
                        cur[k][s] = [ta, pre_cp[1] + [tM[s]]]
                        if ta < CP_area[i]:
                            CP_area[i] = ta
                            CP_points[i] = cur[k][s][1]
        
        prev = cur.copy()
        cur = [[None for k in range(pn)] for j in range(pn)]

# output
for q in query:
    if CP_points[q-3]:
        print(' '.join(map(lambda x: str(M[x]), CP_points[q-3])))
    else:
        print("NA")

Submission status: Accepted

Solution code:
def triangle(a,b,c):
    return cross(b - a, c - a) / 2
def cos(a, b):
    return dot(a, b) / (abs(a) * abs(b))
def dot(a, b):
    return a.real * b.real + a.imag * b.imag
def cross(a, b):
    return a.real * b.imag - a.imag * b.real

def solve(xyi, xy):
    # anti_clockwise_sort:sort(key=cos((1,0)(0,0)(xyk-xyi)),reverse=True)
    xy = sorted(xy, key=lambda a:(a.real - xyi.real) / abs(a - xyi), reverse=True)

    area = [[[(float('inf'),) for i in range(len(xy))] for j in range(len(xy))] for k in range(max(q) + 1)]
    # triangle
    for j in range(len(xy) - 1):
        for k in range(j + 1, len(xy)):
            area[3][j][k] = triangle(xyi, xy[j], xy[k]), (xyi, xy[j], xy[k])

    # square,pentagon,,,, 
    for i in range(4, max(q) + 1):
        for j in range(len(xy)):
            for k in range(j + 1, len(xy)):
                for s in range(k + 1, len(xy)):                                
                    if math.isinf(area[i - 1][j][k][0]):
                        continue                
                    if triangle(xy[j], xy[k], xy[s]) <= 0:
                        continue #skip not convex polygon
                    tmp = area[3][k][s][0] + area[i - 1][j][k][0]
                    if not math.isnan(area[i][k][s][0]) and area[i][k][s][0] < tmp:
                        continue
                    area[i][k][s] = tmp, area[i - 1][j][k][1] + (xy[s],)
                
    min_space = [(float('inf'),)] * (max(q) + 1)
    for i in range(3, max(q) + 1):
        min_space[i] = min([y for x in area[i] for y in x], key=operator.itemgetter(0))
    
    return min_space

import sys
import operator
import math

f = sys.stdin

xy =[[int(i) for i in f.readline().split()] for _ in range(int(f.readline()))]
q = [int(f.readline()) for _ in range(int(f.readline()))]

xy = [x + y * 1j for x, y in xy]
id = {xyi:i+1 for i, xyi in enumerate(xy)}

xy.sort(key=operator.attrgetter('imag', 'real'))

min_space = [(float('inf'),) for i in range(max(q) + 1)]
while 2 < len(xy):
    xyi = xy.pop(0)
    min_space = [min(i, j, key=operator.itemgetter(0)) for i,j in zip(min_space,solve(xyi, xy))]

for qi in q:
    if math.isinf(min_space[qi][0]):
        print('NA')
    else:
        print(*[id[i] for i in min_space[qi][1]])