SRM 679 div1 med: RedAndBluePoints

問題

TopCoder Statistics - Problem Statement

解法

pekempey さんの解法を参考にしました。
pekempey.hatenablog.com

凸多角形は上側と下側に分けることができ, 上側の頂点では, 座標を時計回りに見ていった時, 今見ている頂点 p と次の頂点 q について, (p.x < q.x または (p.x == q.x かつ p.y < q.y)) が成り立ちます。
下側の頂点では (p.x > q.x または (p.x == q.x かつ p.y > q.y)) になります。

ということで, 凸多角形を考えるときは, 頂点を x 座標, y 座標 についてソートしておくと, 頂点を順番に見ていくだけで大丈夫になります。

この順番を使って, dp を組みます。
dpup[last][second][first] = (上側の頂点を調べる時, 最後の頂点が last, 最後から 2 番目の頂点が second, 最初の頂点が first である時の, 含まれる青い頂点の最大数)
dpdown[last][second][first] = (下側の頂点を調べる時, 最後の頂点が last, 最後から 2 番目の頂点が second, 最初の頂点が first である時の, 含まれる青い頂点の最大数)
を求めます(dpup と dpdown は上側を調べる -> 下側を調べる という順番でやれば区別する必要がなくなるので下のコードでは dp という配列で同時にやっています)。

次に含める頂点を next とすると, 多角形の中に含まれる頂点は三角形(first, last, next) の中にあるものを考えれば十分です。この中に赤い頂点があれば ng です。
また, 凸多角形になるかを調べるために, last から next へのベクトルと second から last へのベクトルの行列式を調べる必要があります。

double eps = 1e-6;

double add(double a, double b) {
    if (abs(a+b) < eps * (abs(a)+abs(b))) return 0;
    return a+b;
}

bool equal(double a, double b) {
    return add(a, -b) == 0;
}

struct P {
    double x, y;
    P() {}
    P(double x, double y) : x(x), y(y) {}
    P operator+(P p) const {return P(add(x, p.x), add(y, p.y));}
    P operator-(P p) const {return P(add(x, -p.x), add(y, -p.y));}
    P operator*(double d) const {return P(x*d, y*d);}
    double dot(P p) const {return add(x*p.x, y*p.y);} // 内積
    double det(P p) const {return add(x*p.y, -y*p.x);} // 外積
    double dist(P p) const {return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y));} // 距離
    void normalize() {double d = sqrt(x*x+y*y); x /= d; y /= d;} // 正規化
    bool operator<(const P& rhs) const {
        if (x != rhs.x) return x < rhs.x;
        return y < rhs.y;
    }
    bool operator==(const P& rhs) const {
        return equal(x, rhs.x) && equal(y, rhs.y);
    }
};

// 点 P が三角形の内部に存在するか
bool inTri(P p, P a, P b, P c) {
    P ab = b-a;
    P bc = c-b;
    P ca = a-c;
    P ap = p-a;
    P bp = p-b;
    P cp = p-c;
    double c1 = ab.det(bp);
    double c2 = bc.det(cp);
    double c3 = ca.det(ap);
    if (c1 > 0 && c2 > 0 && c3 > 0) return true;
    if (c1 < 0 && c2 < 0 && c3 < 0) return true;
    return false;
}

const int INF = 1000;
// last, second, first
int dp[55][55][55];
bool ok[55][55][55];
int cnt[55][55][55];

class RedAndBluePoints {
public:
    int find(vector <int> blueX, vector <int> blueY, vector <int> redX, vector <int> redY) {
        int B = blueX.size(), R = redX.size();
        vector<P> bp(B), rp(R);
        for (int i = 0; i < B; i++) bp[i] = P(blueX[i], blueY[i]);
        for (int i = 0; i < R; i++) rp[i] = P(redX[i], redY[i]);
        sort(bp.begin(), bp.end());
        memset(cnt, 0, sizeof(cnt));
        for (int i = 0; i < B; i++) for (int j = 0; j < B; j++) {
            for (int k = 0; k < B; k++) dp[i][j][k] = -INF;
        }
        for (int i = 0; i < B; i++) for (int j = 0; j < B; j++) for (int k = 0; k < B; k++) {
            ok[i][j][k] = true;
            if (i == j || j == k || k == i) continue;
            for (int l = 0; l < R; l++) ok[i][j][k] &= !inTri(rp[l], bp[i], bp[j], bp[k]);
            for (int l = 0; l < B; l++) cnt[i][j][k] += inTri(bp[l], bp[i], bp[j], bp[k]);
        }
        // left to right
        for (int last = 0; last < B; last++) {
            dp[last][last][last] = 0;
            for (int first = 0; first < B; first++) {
                for (int second = 0; second < B; second++) if (dp[last][second][first] > -INF) {
                    for (int next = last+1; next < B; next++) {
                        if (!ok[first][last][next]) continue;
                        P e1 = bp[last]-bp[second];
                        P e2 = bp[next]-bp[last];
                        if (e1.det(e2) < 0) continue;
                        dp[next][last][first] = max(dp[next][last][first], dp[last][second][first]+cnt[first][last][next]+1);
                    }
                }
            }
        }
        // left to right
        for (int last = B-1; last >= 0; last--) {
            for (int first = B-1; first >= 0; first--) {
                for (int second = B-1; second >= 0; second--) if (dp[last][second][first] > -INF) {
                    for (int next = last-1; next >= 0; next--) {
                        if (!ok[first][last][next]) continue;
                        P e1 = bp[last]-bp[second];
                        P e2 = bp[next]-bp[last];
                        if (e1.det(e2) < 0) continue;
                        dp[next][last][first] = max(dp[next][last][first], dp[last][second][first]+cnt[first][last][next]+1);
                    }
                }
            }
        }
        int ret = min(2, B);
        for (int i = 0; i < B; i++) for (int j = 0; j < B; j++) {
            ret = max(ret, dp[i][j][i] + (i==j));
        }
        return ret;
    }
};