Hướng dẫn cho LQDOJ CUP 2022 - Round 7 - TRICOVER

Chỉ sử dụng khi thực sự cần thiết như một cách tôn trọng tác giả và người viết hướng dẫn này.

Chép code từ bài hướng dẫn để nộp bài là hành vi có thể dẫn đến khóa tài khoản.

Authors: letangphuquy

Subtask \(1\) (\(5\%\) số điểm): \(n = 1\).

Tutorial

Chỉ cần tính diện tích của tam giác. Giải trong \(O(1)\).
Không nên dùng công thức Heron vì sai số lớn, nên sử dụng tích chéo.

Solution

C++

#include<bits/stdc++.h>

using namespace std;

const long double EPS = 1e-9;

bool equals(const long double &x, const long double &y) { 
    return abs(x - y) <= EPS; 
}

bool lesser(const long double &x, const long double &y) { 
    return !equals(x, y) && x < y; 
}

bool leq(const long double &x, const long double &y) { 
    return equals(x, y) or x < y; 
}

const long long INF = 1e18;
const long double PI = acos(-1);

//Lời giải có tham khảo KACTL Notebook
template<typename T> struct Point {
    T x,y;

    Point(T x = 0, T y = 0): x(x), y(y) {}

    bool operator < (const Point &p) const { 
        return equals(x, p.x) ? lesser(y, p.y) : lesser(x, p.x); 
    }

    T cross(const Point &p) const { 
        return x * p.y - y * p.x; 
    }

    T cross(const Point &a, const Point &b) const { 
        return (a - *this).cross(b - *this); 
    } 
};

long double areaOf(const vector<Point<long double>> &H) {
    long double area = 0;
    int n = H.size();
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        area += H[i].cross(H[j]);
    }
    return abs(area) / 2;
}

void subtask1() {
    int numAreas, width, height; 
    cin >> numAreas >> width >> height;
    long double answer = 0;
    for (int i = 1; i <= numAreas; i++) {
        vector<Point<long double>> tri(3);
        for (Point<long double>& p : tri) {
            cin >> p.x >> p.y;
        }
        answer += areaOf(tri);
    }
    cout << fixed << setprecision(15) << answer;
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0); 
    cout.tie(0);

    freopen("TRICOVER.inp", "r", stdin);
    freopen("TRICOVER.out", "w", stdout);

    subtask1();

    return 0;
}

Subtask \(2\) (\(30\%\) số điểm): \(n \leq 14\).

Tutorial

Với \(n\) nhỏ, ta nghĩ đến bao hàm loại trừ.
Để ý rằng:

Mỗi tam giác là một đa giác lồi
Khi ta "cắt" ("cắt" = lấy phần chung) một đa giác lồi bởi một đa giác lồi khác, ta thu được phần giao cũng là một đa giác lồi (hoặc rỗng).

Như vậy bài toán cần giải quyết ở đây: Tìm diện tích phần giao của một tập các tam giác \(\Rightarrow\) "Cắt" một đa giác bằng một tam giác.
Giả sử cần lấy phần chung của đa giác và tam giác. Xét các trường hợp như sau:

Một trong hai hình bao chứa hình còn lại \(\Rightarrow\) phần giao là hình bé hơn
Hai hình không có điểm chung (cạnh của hình này không cắt cạnh của hình kia) \(\Rightarrow\) phần giao là rỗng
Hai hình có điểm chung: Xét tập \(S\) gồm các điểm: Giao điểm của các cạnh, các đỉnh của hình này mà nằm trọn trong hình kia. Phần giao chính là đa giác có các đỉnh trong tập \(S\). Nhận thấy cách tính này bao quát luôn cả 2 trường hợp đặc biệt trên.

Độ phức tạp: \(\mathcal{O}(2^n \times n^2)\)

Solution

C++

#include<bits/stdc++.h>

using namespace std;

const long double EPS = 1e-9;

bool equals(const long double &x, const long double &y) { 
    return abs(x - y) <= EPS; 
}

bool lesser(const long double &x, const long double &y) { 
    return !equals(x, y) && x < y; 
}

bool leq(const long double &x, const long double &y) { 
    return equals(x, y) || x < y; 
}

const long long INF = 1e18;
const long double PI = acos(-1);

//Lời giải có tham khảo KACTL Notebook
template<typename T> struct Point {
    T x,y;

    Point(T x = 0, T y = 0): x(x), y(y) {}

    bool operator < (const Point &p) const { 
        return equals(x, p.x) ? lesser(y, p.y) : lesser(x, p.x); 
    }

    Point operator - (const Point &p) const { 
        return Point(x - p.x, y - p.y); 
    }

    Point operator + (const Point &p) const { 
        return Point(x + p.x, y + p.y); 
    }

    Point operator * (const T &k) const { 
        return Point(x * k, y * k); 
    }

    Point operator / (const T &k) const { 
        return Point(x / k, y / k); 
    }

    T dot(const Point &p) const { 
        return x * p.x + y * p.y; 
    }

    T cross(const Point &p) const { 
        return x * p.y - y * p.x; 
    }

    T cross(const Point &a, const Point &b) const { 
        return (a - *this).cross(b - *this); 
    } 
};

int sign(long double x) { 
    return (x > 0) - (x < 0); 
}

template<typename P> bool onSegment(P s, P e, P p) {
    return equals(p.cross(s,e), 0) && leq((s - p).dot(e - p), 0);
}

template<typename P> vector<P> segInter(P a, P b, P c, P d) {
    auto oa = c.cross(d,a), ob = c.cross(d,b);
    auto oc = a.cross(b,c), od = a.cross(b,d);
    if (sign(oa) * sign(ob) < 0 && sign(oc) * sign(od) < 0) {
        return {(a * ob - b * oa) / (ob - oa)};
    }
    set<P> s;
    if (onSegment(c, d, a)) {
        s.insert(a);
    }
    if (onSegment(c, d, b)) {
        s.insert(b);
    }
    if (onSegment(a, b, c)) {
        s.insert(c);
    }
    if (onSegment(a, b, d)) {
        s.insert(d);
    }
    return {s.begin(), s.end()};
}

const int N = 333;
const int V = (3 * N) * (3 * N);

long double areaOf(const vector<Point<long double>> &H) {
    long double area = 0;
    int n = H.size();
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        area += H[i].cross(H[j]);
    }
    return abs(area) / 2;
}

bool isInside(vector<Point<long double>> H, Point<long double> P) {
    if (H.empty()) {
        return false;
    }
    long double sumArea = 0;
    int n = H.size();
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        sumArea += areaOf(vector<Point<long double>>({P,H[i],H[j]}));
    }
    return equals(abs(sumArea), areaOf(H));
};

void regulate(vector<Point<long double>>& A) {
    if (A.empty()) {
        return;
    }
    sort(A.begin(), A.end());
    sort(A.begin()+1, A.end(), [&](Point<long double> u, Point<long double> v) {
        return A[0].cross(u, v) > 0;
    });
}

vector<Point<long double>> intersect(vector<Point<long double>> A, vector<Point<long double>> B) {
    vector<Point<long double>> res;
    for (Point<long double> p : A) {
        if (isInside(B, p)) {
            res.push_back(p);
        }
    }
    for (Point<long double> p : B) {
        if (isInside(A, p)) {
            res.push_back(p);
        }
    }
    int sza = A.size(), szb = B.size();
    for (int ia = 0; ia < sza; ia++) {
        int ja = (ia + 1) % sza;
        for (int ib = 0; ib < szb; ib++) {
            int jb = (ib + 1) % szb;
            vector<Point<long double>> get = segInter(A[ia], A[ja], B[ib], B[jb]);
            if (get.size() != 1) {
                continue;
            }
            res.push_back(get[0]);
        }
    }
    regulate(res);
    return res;
}

void subtask2() {
    int numAreas, width, height; 
    cin >> numAreas >> width >> height;
    vector<vector<Point<long double>>> areas;
    for (int i = 1; i <= numAreas; i++) {
        vector<Point<long double>> tri(3);
        for (Point<long double>& p : tri) {
            cin >> p.x >> p.y;
        }
        regulate(tri);
        areas.push_back(tri);
    }
    long double answer = 0;
    for (int msk = 1; msk < (1<<numAreas); msk++) {
        int first = __builtin_ctz(msk);
        vector<Point<long double>> common = areas[first];
        for (int i = first+1; i < numAreas; i++) {
            if (msk>>i & 1) {
                common = intersect(common, areas[i]);
            }
        }
        int sgn = (__builtin_popcount(msk) & 1) ? +1 : -1;
        answer += areaOf(common) * sgn;
    }
    answer /= numAreas;
    cout << fixed << setprecision(15) << answer;
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0); 
    cout.tie(0);

    freopen("TRICOVER.inp", "r", stdin);
    freopen("TRICOVER.out", "w", stdout);

    subtask2();

    return 0;
}

Subtask \(3\) (\(25\%\) số điểm): Các tam giác là tam giác vuông cân có hai cạnh song song với trục tọa độ, và chỉ thuộc vào một trong hai dạng sau:

Tutorial

Nếu bạn biết bài https://oj.vnoi.info/problem/triangle thì việc giải subtask này không khó.
Ý tưởng: Dùng sweep line: Quét qua từng tọa độ \(x\), với mỗi tọa độ \(x_0\), ta tính tổng diện tích được giới hạn trong hai đường thẳng \(x = x_0\) và \(x = x_0+1\). Tạm gọi phần mặt phẳng bị giới hạn bởi 2 đường dọc này là một "lát cắt".
Vì các tam giác vuông cân chia đôi các ô vuông đơn vị nên phần diện tích tại mỗi lát cắt không phải là số nguyên.
Để dễ xử lý, ta "gấp đôi" trục Oy. Lúc này, toàn bộ các diện tích đều là số nguyên. Quan sát hình sau:

nửa dưới của ô vuông đơn vị tương ứng với ô vuông lẻ, và nửa trên tương ứng với ô có tọa độ chẵn.
Ví dụ:

Tư tưởng của thuật toán sweep line: Tại tọa độ \(x_0\), tôi dùng CTDL để lưu trữ trạng thái ở lát cắt đó. Khi đi qua tọa độ \(x_0+1\), tôi chỉ cập nhật những gì cần thiết để có được trạng thái tại lát cắt \(x_0+1\). Trong trường hợp này, ta có thể dùng Segment tree Lazy update, hoặc đơn giản hơn cả là prefix sum để quản lý. (Tăng một đoạn liên tiếp lên 1 đơn vị, đếm số vị trí dương).
Vì tọa độ rất bé nên không cần nén tọa độ hoặc sử dụng thêm CTDL khác.

Độ phức tạp: \(\mathcal{O}((W+H+n)^2)\) hoặc \(O(n.\log), \dots\) tùy cách cài đặt.

Solution

C++

#include<bits/stdc++.h>

using namespace std;

const int MAX = 1e4 + 50;

struct Point {
    int x, y;

    Point() {}

    Point(int _x, int _y) {
        x = _x; y = _y;
    }

    bool operator < (const Point &p) const {
        return x < p.x || (x == p.x && y < p.y);
    }
};

vector<pair<int,int>> adds[MAX];

void subtask3() {
    int numAreas, width, height; 
    cin >> numAreas >> width >> height;
    vector<vector<Point>> triangles;
    for (int i = 1; i <= numAreas; i++) {
        vector<Point> tri(3);
        for (Point& p : tri) {
            cin >> p.x >> p.y;
        }
        triangles.push_back(tri);
    }
    //tọa độ của ô = tọa độ của góc trái trên 
    for (vector<Point> p : triangles) {
        sort(p.begin(), p.end());
        if (p[0].x == p[1].x) {
            assert(p[2].y == p[0].y);
            int r = 2 * p[1].y - 1;
            for (int x = p[1].x; x < p[2].x; x++, r -= 2) {
                adds[x].push_back({2*p[0].y + 1, r});
            }
        } else {
            assert(p[2].y == p[0].y);
            int l = 2 * p[0].y, r = l;
            for (int x = p[0].x; x < p[1].x; x++, l -= 2) {
                adds[x].push_back({l,r});
            }
        }
    }
    long double answer = 0;
    for (int x = 0; x <= width; x++) {
        vector<int> psum(2 * height + 10, 0);
        for (pair<int,int> a : adds[x]) {
            psum[a.first]++; 
            psum[a.second+1]--;
        }
        for (int y = 1; y <= 2 * height + 4; y++) {
            psum[y] += psum[y-1];
            if (psum[y]) {
                answer++;
            }
        }
    }
    cout << fixed << setprecision(15) << answer / (2*numAreas);
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0); 
    cout.tie(0);

    freopen("TRICOVER.inp", "r", stdin);
    freopen("TRICOVER.out", "w", stdout);

    subtask3();

    return 0;
}

Bonus

Làm sao để từ input, xác định được đó là test thuộc subtask 3
Đoạn code mẫu nói trên chưa hoàn chỉnh lắm vì chỉ đọc vào và giải mỗi subtask này thôi chứ không kiểm tra. Trong thực tế, có thể bạn sẽ muốn chạy nhiều lời giải khác nhau cho các subtask khác nhau, nên cần bước này.

C++

//lưu ý: đọc input trước khi kiểm tra
bool ok = true;
for (auto p : triangles) {
    ok &= p.size() == 3;
    sort(all(p));
    ok &= ((p[0].x == p[1].x && p[1].x < p[2].x)
        or (p[0].x < p[1].x && p[1].x == p[2].x))
        and (p[0].y == p[2].y);
}
if (ok) {
    subtask3(); 
}

Subtask \(4\) (\(40\%\) số điểm): Không có ràng buộc gì thêm.

Tutorial

Quan sát: \(n\) tam giác cắt mặt phẳng thành nhiều vùng khác nhau. Ta muốn lấy ra được toàn bộ các vùng thuộc phần diện tích bị bao phủ và từ đó tính được tổng diện tích bao phủ (bằng tổng diện tích các vùng).
Tuy nhiên, điểm khó ở đây là các vùng không phải lúc nào cũng là đa giác lồi, mà có thể là đa giác lõm, hoặc thậm chí là hình vành khuyên!!. Quan sát:

Trong trường hợp lý tưởng, nếu toàn bộ các vùng tạo thành đều là đa giác, có thể làm thuật toán như sau:
Duyệt qua \(O(n^2)\) cặp cạnh để lấy ra các giao điểm. Coi mỗi giao điểm như một đỉnh của đồ thị. Ta sẽ dựng một đồ thị phẳng mô phỏng các vùng.
Tiến hành duyệt qua từng cạnh của tam giác, với mỗi đoạn thẳng (cạnh), nối 2 đỉnh liên tiếp cùng nằm trên đoạn thẳng đó.

Với mỗi cạnh \((u,v)\) thuộc đồ thị, ta thêm cạnh (cạnh vô hướng) này vào đa tập \(S\) một hoặc hai lần:

Nếu cạnh này là cạnh biên: một bên của nó chứa phần diện tích được bao phủ, bên còn lại là mặt ngoài (không được phủ bởi các tam giác): thêm vào một lần
Ngược lại, thêm vào hai lần (vì cạnh này nằm trọn trong phần bao phủ nên sẽ thuộc về 2 đa giác khác nhau).

Để lấy ra được các vùng, ta làm như sau: Ưu tiên chọn các đỉnh bắt đầu theo thứ tự tăng dần về \(x\), rồi tới \(y\). Bắt đầu từ một đỉnh \(u\), chọn một đỉnh \(v\) kề \(u\), mà \((u,v) \in S\) và hệ số góc của \(v-u\) nhỏ nhất có thể. Đây sẽ là cạnh đầu tiên của đa giác ta muốn lấy ra. Loại \((u,v)\) khỏi \(S\). Sau đó, lặp lại quá trình sau:

Chọn đỉnh \(w\) kề \(v\) mà \((w,v) \in S, w \neq u\) sao cho góc tạo bởi cạnh này và cạnh liền trước là nhỏ nhất có thể. Tức là: góc \(\alpha\) tạo bởi vector \(v-u\) và vector \(w-v\) lớn nhất (vì \(\pi-\alpha\) chính là góc giữa 2 cạnh này). Nếu không tìm thấy thì kết thúc.
Xóa cạnh \((w,v)\) khỏi \(S\).
Nếu \(w\) trùng với đỉnh ta chọn lúc đầu thì kết thúc. Lúc này ta đã lấy ra được một đa giác.
Gán \(u \leftarrow v, v \leftarrow w\), và quay lại các bước trên.

Lưu ý: với một cạnh được thêm vào 2 lần thì 2 đa giác chứa nó sẽ đi qua cạnh đó theo 2 hướng khác nhau (\(u \rightarrow v\) và \(v \rightarrow u\)).
Với mỗi vùng được phát hiện, ta cộng diện tích của nó vào tổng.

Làm sao để đưa về trường hợp lý tưởng như trên?
Cách 1) Diện tích có dấu
Việc ưu tiên đỉnh có tọa độ bé hơn trước làm đỉnh bắt đầu, cũng như chọn cạnh có hệ số góc giúp ta đảm bảo sẽ duyệt qua các vùng theo chiều ccw (ngược kim đồng hồ, rẽ trái) trước. Diện tích duyệt theo chiều này sẽ là dương. Nếu sau đó vẫn còn dư các cạnh, và duyệt được một đa giác theo chiều ngược lại (rẽ phải) thì diện tích sẽ âm. Ta cộng vào như bình thường, coi như là bù trừ của các hình.

Cách 2) Kéo dài các cạnh thành các đường thẳng.
Bây giờ ta sẽ làm thuật toán như trên, nhưng xét với toàn bộ \(3n \choose 2\) giao điểm. Nhận xét: Hình được tạo ra từ "sự cắt" của các đường thẳng sẽ chỉ bao gồm các đa giác lồi.

Độ phức tạp: \(\mathcal{O}((3n)^2.\log + n^3)\)

Solution

C++

#include <bits/stdc++.h>

using namespace std;

template<class A, class B> bool umin(A& var, const B& val) {
    return (val < var) ? (var = val, true) : false;
}

const long double EPS = 1e-9;

bool equals(const long double& x, const long double& y) {
    return abs(x - y) <= EPS; 
}

bool lesser(const long double& x, const long double& y) { 
    return !equals(x, y) && x < y; 
}

bool leq(const long double& x, const long double& y) { 
    return equals(x, y) || x < y; 
}

const long double PI = acos(-1);

template<typename T> struct Point {
    T x, y;

    Point(T x = 0, T y = 0): x(x), y(y) {}

    bool operator < (const Point& p) const { 
        return equals(x, p.x) ? lesser(y, p.y) : lesser(x, p.x); 
    }

    Point operator - (const Point& p) const { 
        return Point(x - p.x, y - p.y); 
    }

    Point operator + (const Point& p) const { 
        return Point(x + p.x, y + p.y); 
    }

    Point operator * (const T& k) const { 
        return Point(x * k, y * k); 
    }

    Point operator / (const T& k) const { 
        return Point(x / k, y / k); 
    }

    T dot(const Point& p) const { 
        return x * p.x + y * p.y; 
    }

    T cross(const Point& p) const { 
        return x * p.y - y * p.x; 
    }

    T cross(const Point& a, const Point& b) const { 
        return (a - *this).cross(b - *this); 
    } 

    long double angle(void) const { 
        long double al = atan2(y,x); 
        if (lesser(al,0)) {
            al += 2*PI;
        }
        return al; 
    }

    Point rotate(const long double& a) const {
        return Point(x * cos(a) - y * sin(a), x * sin(a) + y * cos(a));
    }

    friend ostream& operator << (ostream& os, const Point& p) {
        return os << p.x << ' ' << p.y; 
    }
};

int sign(long double x) { 
    return (x > 0) - (x < 0); 
}

vector<Point<long double>> points;

struct Comparator {
    bool operator() (const int& i, const int& j) const {
        return points[i] < points[j];
    }
};

typedef array<Point<long double>,3> Triangle;
const int N = 333;
const int V = (3 * N) * (3 * N);
multiset<int> adj[V];

void connect(int u, int v) {
    adj[u].insert(v);
    adj[v].insert(u);
}

void deconnect(int u, int v) {
    adj[u].erase(adj[u].find(v));
    adj[v].erase(adj[v].find(u));
}

long double angleBetween(const Point<long double>& u, const Point<long double>& v) {
    long double al = v.rotate(-u.angle()).angle();
    al = PI - al;
    if (lesser(al, 0)) {
        al += PI * 2;
    }
    if (leq(PI * 2, al)) {
        al -= PI * 2;
    }
    return al;
}

long double refine(long double a) {
    a += PI / 2;
    if (leq(PI * 2, a)) {
        a -= PI * 2;
    }
    return a;
}

template<typename P> pair<int,P> lineInter(P s1, P e1, P s2, P e2) {
    auto det = (e1-s1).cross(e2-s2);
    if (!det) {
        return {-(s1.cross(e1,s2) == 0), P(0, 0)};
    }
    auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
    return {1, (s1 * p + e1 * q) / det};
}

void solve() {
    int numAreas, width, height; 
    cin >> numAreas >> width >> height;
    vector<Triangle> areas;
    assert(!areas.size());  
    for (int i = 1; i <= numAreas; i++) {
        Triangle tri;
        for (Point<long double>& p : tri) {
            cin >> p.x >> p.y;
        }
        areas.push_back(tri);
    }

    auto isInside = [&](int id, Point<long double> P) {
        int s0 = sign(P.cross(areas[id][0], areas[id][1])); 
        int s1 = sign(P.cross(areas[id][1], areas[id][2])); 
        if (s0 != s1) {
            return false;
        }
        int s2 = sign(P.cross(areas[id][2], areas[id][0]));
        return s1 == s2;
    };
    auto isRelevant = [&](Point<long double> P) {
        for (int i = 0; i < numAreas; i++) {
            if (isInside(i, P)) {
                return true;
            }
        }
        return false;
    };

    map<Point<long double>, int> getIndex;
    points.clear();
    points.push_back(Point<long double>()); // 1-based index
    int numNodes = 0;
    auto getNode = [&](Point<long double> P) {
        int& idx = getIndex[P];
        if (!idx) {
            idx = ++numNodes;
            points.push_back(P);
        }
        return idx;
    };

    set<int, Comparator> verOnEdge[numAreas][3]; 
    // vì không có 3 điểm nào thẳng hàng nên xét theo cạnh
    for (int i = 0; i < numAreas; i++) {
        for (int j = 0; j < 3; j++) {
            int idx = getNode(areas[i][j]);
            verOnEdge[i][j].insert(idx);
            verOnEdge[i][(j - 1 + 3) % 3].insert(idx);
        }
    }
    cerr << fixed << setprecision(4);
    const long double shift = 1e-6;
    int numEdges = 0;
    for (int i = 0; i < numAreas; i++) {
        for (int j = i+1; j < numAreas; j++) {
            for (int i1 = 0; i1 < 3; i1++) {
                int i2 = (i1 + 1) % 3;
                for (int j1 = 0; j1 < 3; j1++) {
                    int j2 = (j1 + 1) % 3;
                    auto get = lineInter<Point<long double>>(areas[i][i1], areas[i][i2], areas[j][j1], areas[j][j2]);
                    if (get.first != 1) {
                        continue;
                    }
                    int idx = getNode(get.second);
                    verOnEdge[i][i1].insert(idx);
                    verOnEdge[j][j1].insert(idx);
                }
            }
        }
    }

    for (int i = 0; i < numAreas; i++) {
        for (int i1 = 0; i1 < 3; i1++) {
            int last = 0;
            for (auto idx : verOnEdge[i][i1]) {
                if (last) {
                    ++numEdges;
                    connect(last, idx);
                    Point<long double> mid = (points[idx] + points[last]) / 2;
                    Point<long double> dir = points[idx] - points[last];
                    swap(dir.x, dir.y); 
                    dir.x *= -1;
                    bool g1 = isRelevant(mid + dir * shift); 
                    bool g2 = isRelevant(mid - dir * shift); 
                    if (g1 == g2) {
                        if (g1) {
                            connect(last, idx);
                        }
                        else {
                            deconnect(last, idx);
                        }
                    }
                }
                last = idx;
            }
        }
    }

    map<pair<int,int>, bool> check;
    vector<int> ids(numNodes);
    iota(ids.begin(), ids.end(), 1);
    sort(ids.begin(), ids.end(), Comparator());

    long double answer = 0;
    for (int idx = 0; idx < numNodes; ) {
        int s = ids[idx];
        if (adj[s].empty()) { 
            idx++; 
            continue; 
        }
        int u = s, v = -1; 
        long double best = 1e9;
        for (auto vv : adj[u]) {
            if (!check[{u, vv}]) {
                long double cur = refine((points[vv] - points[u]).angle());
                if (v == -1 || lesser(cur, best)) {
                    v = vv;
                    best = cur;
                }
            }
        }
        if (v == -1) { 
            idx++; 
            continue; 
        }

        long double A = 0;
        deconnect(u,v);
        check[{u,v}] = true;
        A += points[u].cross(points[v]);

        do {
            if (adj[v].empty()) {
                break;
            }
            int w = -1; 
            long double best = 1e9;

            for (int ww : adj[v]) {
                if (w == s) {
                    break;
                }
                if (ww == u || check[{v, ww}]) {
                    continue;
                }
                long double angle = angleBetween(points[v] - points[u], points[ww] - points[v]);
                if (!equals(angle, 0) && (umin(best, angle) || ww == s)) {
                    w = ww;
                }
            }
            u = v; v = w;
            deconnect(u,v);
            check[{u, v}] = true;
            A += points[u].cross(points[v]);
        } while (v != s);
        if (v == s) {
            answer += abs(A);
        }
    }
    answer /= 2 * numAreas;
    cout << fixed << setprecision(10) << answer;
    for (int u = 0; u <= numNodes; u++) {
        adj[u].clear();
    }
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0); 
    cout.tie(0);

    freopen("TRICOVER.inp", "r", stdin);
    freopen("TRICOVER.out", "w", stdout);

    solve();

    return 0;
}

Hướng tiếp cận \(O(n^2)\)

Cách làm này của 6aren.

Solution

C++

#include <bits/stdc++.h>
using namespace std;

const long double PI = acos(-1);
const long double EPS = 1e-6;

long double sqr(long double x) {
  return x * x;
}

inline int cmp(long double a, long double b) { 
    return (a < b - EPS) ? -1 : ((a > b + EPS) ? 1 : 0); 
}

struct Point {
    long double x, y;

    Point() { 
        x = y = 0.0; 
    }

    Point(long double x, long double y): x(x), y(y) {}

    Point operator + (const Point &a) const { 
        return Point(x + a.x, y + a.y); 
    }

    Point operator - (const Point &a) const { 
        return Point(x - a.x, y - a.y); 
    }

    Point operator * (long double k) const { 
        return Point(x * k, y * k); 
    }
    Point operator / (long double k) const { 
        return Point(x / k, y / k); 
    }

    long double operator * (const Point &a) const { // dot product
        return x * a.x + y * a.y; 
    }

    long double operator % (const Point &a) const { // cross product
        return x * a.y - y * a.x; 
    }  

    long double norm() { 
        return x * x + y * y; 
    }

    long double len() { // hypot(x, y);
        return sqrt(norm()); 
    }  

    Point rotate(long double alpha) {
        long double cosa = cos(alpha), sina = sin(alpha);
        return Point(x * cosa - y * sina, x * sina + y * cosa);
    }
};

int ccw(Point a, Point b, Point c) {
    return cmp((b - a) % (c - a), 0);
}

int cmp_point(Point A, Point B) {
    if (cmp(A.x, B.x) == 0) return cmp(A.y, B.y);
    return cmp(A.x, B.x);
}

int TURN_LEFT = ccw(Point(0, 0), Point(0, 1), Point(-1, 1));
int TURN_RIGHT = ccw(Point(0, 0), Point(0, 1), Point(1, 1));

// NOTE: WILL NOT WORK WHEN a = b = 0.
struct Line {
    long double a, b, c;
    Point A, B;  // Added for polygon intersect line. Do not rely on assumption that these are valid

    Line(long double a, long double b, long double c): a(a), b(b), c(c) {}

    Line(Point A, Point B): A(A), B(B) {
        a = B.y - A.y;
        b = A.x - B.x;
        c = -(a * A.x + b * A.y);
    }
    Line(Point P, long double m) {
        a = -m;
        b = 1;
        c = -((a * P.x) + (b * P.y));
    }
    long double f(Point A) { 
        return a * A.x + b * A.y + c; 
    }
};

bool areParallel(Line l1, Line l2) { 
    return cmp(l1.a * l2.b, l1.b * l2.a) == 0; 
}

bool areSame(Line l1, Line l2) {
  return areParallel(l1, l2) && cmp(l1.c * l2.a, l2.c * l1.a) == 0 && cmp(l1.c * l2.b, l1.b * l2.c) == 0;
}

bool areIntersect(Line l1, Line l2, Point &p) {
  if (areParallel(l1, l2)) {
    return false;
  }
  long double dx = l1.b * l2.c - l2.b * l1.c;
  long double dy = l1.c * l2.a - l2.c * l1.a;
  long double d = l1.a * l2.b - l2.a * l1.b;
  p = Point(dx / d, dy / d);
  return true;
}

vector<Point> polygon_cut(const vector<Point> &P, Line l) {
    vector<Point> Q;
    vector<Point> res;
    for (int i = 0; i < P.size(); ++i) {
        Point A = P[i], B = (i == P.size() - 1) ? P[0] : P[i + 1];
        if (ccw(l.A, l.B, A) != -1) {
            Q.push_back(A);
        }
        if (ccw(l.A, l.B, A) * ccw(l.A, l.B, B) < 0) {
            Point p;
            areIntersect(Line(A, B), l, p);
            Q.push_back(p);
            res.push_back(p);
        }
    }
    if (res.size() == 0) {
        return vector<Point>();
    }
    else {
        return vector<Point>{res.front(), res.back()};
    }
}


int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);

    freopen("TRICOVER.inp", "r", stdin);
    freopen("TRICOVER.out", "w", stdout);

    int n, w, h;
    cin >> n >> w >> h;

    vector<vector<Point>> ps(n, vector<Point>(3));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < 3; j++) {
            cin >> ps[i][j].x >> ps[i][j].y;
        }
        if (ccw(ps[i][0], ps[i][1], ps[i][2]) != TURN_LEFT) {
            swap(ps[i][1], ps[i][2]);
        }
    }

    long double res = 0;

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < 3; j++) {
            Point A = ps[i][j], B = ps[i][(j + 1) % 3];
            Line line(A, B);
            bool is_positive = A.x > B.x;
            if (cmp_point(A, B) > 0) {
                swap(A, B);
            }
            // find intersects
            vector<pair<Point, bool>> cut_points;
            for (int k = 0; k < n; k++) {
                if (k == i) {
                    continue;
                }
                auto cuts = polygon_cut(ps[k], line);
                if (cuts.size() == 0) {
                    continue;
                }
                auto u = cuts[0], v = cuts[1];
                if (cmp_point(u, v) > 0) {
                    swap(u, v);
                }
                if (cmp_point(u, A) < 0) {
                    u = A;
                }
                if (cmp_point(v, B) > 0) {
                    v = B;
                }
                if (cmp_point(u, v) > 0) {
                    continue;
                }
                cut_points.emplace_back(u, true);
                cut_points.emplace_back(v, false);
            }
            cut_points.emplace_back(A, true);
            cut_points.emplace_back(A, false);
            cut_points.emplace_back(B, true);
            cut_points.emplace_back(B, false);

            sort(cut_points.begin(), cut_points.end(), [&](pair<Point, bool> lhs, pair<Point, bool> rhs) {
                return cmp_point(lhs.first, rhs.first) < 0;
            });
            int cur = 0;
            Point last = A;
            for (pair<Point, bool> p : cut_points) {
                if (cur == 0) {
                    res += (is_positive ? 1 : -1) * (p.first.x - last.x) * (p.first.y + last.y) / 2;
                }
                cur += (p.second ? 1 : -1);
                last = p.first;
            }
        }
    }  

    cout << setprecision(10) << fixed << res / n << endl;

    return 0;
}

Bạn có thể tham khảo thêm:

Bình luận

Không có bình luận nào.