Hướng dẫn cho SGAME5

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: SPyofgame

Subtask 1: \(n = 1\)

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.

Code mẫu

C++

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

#define all(v) (v).begin(), (v).end()
#define cst(T) const T&

typedef long long Int;
typedef long double Real;

const Real EPS = 1e-9;
bool equals(cst(Real) x, cst(Real) y) { return abs(x-y) <= EPS; }
bool lesser(cst(Real) x, cst(Real) y) { return !equals(x,y) && x<y; }
bool leq(cst(Real) x, cst(Real) y) { return equals(x,y) or x<y; }

const Int INF = 1e18;
const Real 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< (cst(Point) p) const { return equals(x,p.x) ? lesser(y,p.y) : lesser(x,p.x); }
    T cross(cst(Point) p) const { return x*p.y - y*p.x; }
    T cross(cst(Point) a, cst(Point) b) const { return (a-*this).cross(b-*this); } 
};

using Pnt = Point<Real>;
typedef vector<Pnt> Polygon;

Real areaOf(cst(Polygon) H) {
    Real 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;
    Real answer = 0;
    for (int i = 1; i <= numAreas; i++) {
        Polygon tri(3);
        for (auto& p : tri) cin >> p.x >> p.y;
        answer += areaOf(tri);
    }
    cout << fixed << setprecision(15) << answer;
}

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

    #define task "tricover"
    if (fopen(task".inp", "r")) {
        freopen(task".inp", "r", stdin);
        freopen(task".out", "w", stdout);
    }
    subtask1();
}

Subtask 2: \(7 \le n \le 14\)

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 \(O(2^n \times n^2)\)

Code tham khảo

C++

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

#define all(v) (v).begin(), (v).end()
#define cst(T) const T&

typedef long long Int;
typedef long double Real;

const Real EPS = 1e-9;
bool equals(cst(Real) x, cst(Real) y) { return abs(x-y) <= EPS; }
bool lesser(cst(Real) x, cst(Real) y) { return !equals(x,y) && x<y; }
bool leq(cst(Real) x, cst(Real) y) { return equals(x,y) or x<y; }

const Int INF = 1e18;
const Real 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< (cst(Point) p) const { return equals(x,p.x) ? lesser(y,p.y) : lesser(x,p.x); }
    Point operator- (cst(Point) p) const { return Point(x-p.x, y-p.y); }
    Point operator+ (cst(Point) p) const { return Point(x+p.x, y+p.y); }
    Point operator* (cst(T) k) const { return Point(x*k, y*k); }
    Point operator/ (cst(T) k) const { return Point(x/k, y/k); }
    T dot(cst(Point) p) const { return x*p.x + y*p.y; }
    T cross(cst(Point) p) const { return x*p.y - y*p.x; }
    T cross(cst(Point) a, cst(Point) b) const { return (a-*this).cross(b-*this); } 
};

int sign(Real 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),
        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 {all(s)};
}

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

using Pnt = Point<Real>;
typedef vector<Pnt> Polygon;

Real areaOf(cst(Polygon) H) {
    Real 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(Polygon H, Pnt P) {
    if (H.empty()) return false;
    Real sumArea = 0;
    int n = H.size();
    for (int i = 0; i < n; i++) {
        int j = (i+1)%n;
        sumArea += areaOf(Polygon({P,H[i],H[j]}));
    }
    return equals(abs(sumArea), areaOf(H));
};
void regulate(Polygon& A) {
    if (A.empty()) return ;
    sort(all(A));
    sort(A.begin()+1, A.end(), [&] (Pnt u, Pnt v) {
        return A[0].cross(u,v) > 0;
    });
}
Polygon intersect(Polygon A, Polygon B) {
    Polygon res;
    for (auto p : A) if (isInside(B,p)) res.push_back(p);
    for (auto 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;
            auto 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<Polygon> areas;
    for (int i = 1; i <= numAreas; i++) {
        Polygon tri(3);
        for (auto& p : tri) cin >> p.x >> p.y;
        regulate(tri);
        areas.push_back(tri);
    }
    Real answer = 0;
    for (int msk = 1; msk < (1<<numAreas); msk++) {
        int first = __builtin_ctz(msk);
        Polygon 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_base::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);

    #define task "tricover"
    if (fopen(task".inp", "r")) {
        freopen(task".inp", "r", stdin);
        freopen(task".out", "w", stdout);
    }
    subtask2();
}

Subtask 3:

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 \(O((W+H+n)^2)\) hoặc \(O(n.\log), \dots\) tùy cách cài đặt.

Code mẫu

C++

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

#define all(v) (v).begin(), (v).end()
#define cst(T) const T&

using Pnt = pair<int,int>;
#define x first
#define y second
typedef vector<Pnt> Polygon;

const int MAX = 1e4 + 50;
vector<pair<int,int>> adds[MAX];

void subtask3() {
    int numAreas, width, height; 
    cin >> numAreas >> width >> height;
    vector<Polygon> triangles;
    for (int i = 1; i <= numAreas; i++) {
        Polygon tri(3);
        for (auto& 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 (auto p : triangles) {
        sort(all(p));
        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 (auto [l,r] : adds[x]) {
            psum[l]++; psum[r+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_base::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);

    #define task "tricover"
    if (fopen(task".inp", "r")) {
        freopen(task".inp", "r", stdin);
        freopen(task".out", "w", stdout);
    }
    subtask3();
}

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: \(n \le 1000\)

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 \(O((3n)^2.\log + n^3)\)

Code tham khảo

C++

#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")

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

#define all(v) (v).begin(), (v).end()
#define cst(T) const T&

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

typedef long long Int;
typedef long double Real;

const Real EPS = 1e-9;
bool equals(cst(Real) x, cst(Real) y) { return abs(x-y) <= EPS; }
bool lesser(cst(Real) x, cst(Real) y) { return !equals(x,y) && x<y; }
bool leq(cst(Real) x, cst(Real) y) { return equals(x,y) or x<y; }

#define DBG(x) cerr << #x << " = " << x << ' '
#define DBGn(x) cerr << #x << " = " << x << endl

/*CONSTANTS*/
const Real PI = acos(-1);

template<typename T> struct Point {
    T x,y;
    Point(T x = 0, T y = 0): x(x), y(y) {}
    bool operator< (cst(Point) p) const { return equals(x,p.x) ? lesser(y,p.y) : lesser(x,p.x); }
    Point operator- (cst(Point) p) const { return Point(x-p.x, y-p.y); }
    Point operator+ (cst(Point) p) const { return Point(x+p.x, y+p.y); }
    Point operator* (cst(T) k) const { return Point(x*k, y*k); }
    Point operator/ (cst(T) k) const { return Point(x/k, y/k); }
    T dot(cst(Point) p) const { return x*p.x + y*p.y; }
    T cross(cst(Point) p) const { return x*p.y - y*p.x; }
    T cross(cst(Point) a, cst(Point) b) const { return (a-*this).cross(b-*this); } 
    Real angle(void) const { Real al = atan2(y,x); if (lesser(al,0)) al += 2*PI; return al; }
    Point rotate(cst(Real) a) const {
        return Point(x*cos(a) - y*sin(a), x*sin(a) + y*cos(a));
    }
    friend ostream& operator<< (ostream& os, cst(Point) p) {
        return os << p.x << ' ' << p.y; 
    }
};

int sign(Real 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),
        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 {all(s)};
}

using Pnt = Point<Real>;
vector<Pnt> points;
struct Comparator {
    bool operator() (cst(int) i, cst(int) j) const {
        return points[i] < points[j];
    }
};

typedef array<Pnt,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));
}

Real angleBetween(cst(Pnt) u, cst(Pnt) v) {
    Real al = v.rotate(-u.angle()).angle();
    al = PI-al;
    if (lesser(al,0)) al += 2*PI;
    if (leq(2*PI,al)) al -= 2*PI;
    return al;
}
Real refine(Real a) {
    a += PI/2;
    if (leq(2*PI, a)) a -= 2*PI;
    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 (auto& p : tri) cin >> p.x >> p.y;
        areas.push_back(tri);
    }

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

    map<Pnt, int> getIndex;
    points.clear();
    points.push_back(Pnt()); //1-based index
    int numNodes = 0;
    auto getNode = [&] (Pnt P) {
        auto& 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 Real 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<Pnt>(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);
                    Pnt mid = (points[idx] + points[last]) / 2;
                    Pnt 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(all(ids), 1);
    sort(all(ids), Comparator());

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

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

        do {
            if (adj[v].empty()) break;
            int w = -1; Real best = 1e9;
            for (auto ww : adj[v]) {
                if (w == s) break;
                if (ww == u or check[{v,ww}]) continue;
                Real angle = angleBetween(points[v] - points[u], points[ww] - points[v]);
                if (!equals(angle,0) and (umin(best, angle) or 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;
    cerr << "answer = " << answer << " vs. area = " << answer*numAreas;
    for (int u = 0; u <= numNodes; u++) adj[u].clear();
}

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

    #define task "tricover"
    if (fopen(task".inp", "r")) {
        freopen(task".inp", "r", stdin);
        freopen(task".out", "w", stdout);
    }
    solve();
}

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

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;

#ifdef LOCAL
#include <cp/debugger.hpp>
#else
#define debug(...) 42
#endif

#define sz(v) ((int)(v).size())
#define all(v) (v).begin(), (v).end()

typedef int64_t int64;
typedef pair<int, int> ii;

#define double long double

const double PI = acos(-1);

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

#define EPS 1e-6
inline int cmp(double a, double b) { return (a < b - EPS) ? -1 : ((a > b + EPS) ? 1 : 0); }
struct Point {
  double x, y;
  Point() { x = y = 0.0; }
  Point(double x, 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*(double k) const { return Point(x * k, y * k); }
  Point operator/(double k) const { return Point(x / k, y / k); }

  double operator*(const Point &a) const { return x * a.x + y * a.y; }  // dot product
  double operator%(const Point &a) const { return x * a.y - y * a.x; }  // cross product
  double norm() { return x * x + y * y; }
  double len() { return sqrt(norm()); }  // hypot(x, y);
  Point rotate(double alpha) {
    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 {
  double a, b, c;
  Point A, B;  // Added for polygon intersect line. Do not rely on assumption that these are valid

  Line(double a, double b, 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, double m) {
    a = -m;
    b = 1;
    c = -((a * P.x) + (b * P.y));
  }
  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;
  double dx = l1.b * l2.c - l2.b * l1.c;
  double dy = l1.c * l2.a - l2.c * l1.a;
  double d = l1.a * l2.b - l2.a * l1.b;
  p = Point(dx / d, dy / d);
  return true;
}

typedef vector< Point > Polygon;

Polygon polygon_cut(const Polygon &P, Line l) {
  Polygon 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(false);
  cin.tie(nullptr);
  #define task "tricover"
  if (fopen(task".inp", "r")) {
    freopen(task".inp", "r", stdin);
    freopen(task".out", "w", stdout);
  }
  int n, w, h;
  cin >> n >> w >> h;
  vector<Polygon> 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]);
  }

  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 (sz(cuts) == 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 (auto [p, is_start]: cut_points) {
        if (cur == 0) {
          res += (is_positive ? 1 : -1) * (p.x - last.x) * (p.y + last.y) / 2;
        }
        cur += (is_start ? 1 : -1);
        last = p;
      }
    }
  }  

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

  return 0;
}

Bình luận

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