Subtask 1: ~M, N \le 80~

Duyệt trâu, thử loại bỏ từng cặp tranh, thực hiện tổng dồn trên bảng ~N \cdot N~ với ~M - 2~ ô còn lại, kết quả là các ô có giá trị bằng ~0~, cập nhật đáp án.

Độ phức tạp ~O(M^2 \cdot N^2)~.

Subtask 2: ~M \le 2000~

• Sử dụng tổng dồn để tính ~a[i][j]~ là số hình chữ nhật phủ ô ~i, j~.

• Duyệt tất cả các cặp tranh hình chữ nhật

• Với mỗi cặp hình chữ nhật ~(A, B)~

  • Ta chỉ quan tâm những ô vuông được phủ bởi 1 hoặc 2 bức tranh.

  • Kết quả là (số ô có ~a[i][j]=1~ của A) + (số ô có ~a[i][j]=1~ của B) + (số ô có ~a[i][j]=2~ của (~A \cap B~)) + (số ô có ~a[i][j] = 0~).

Độ phức tạp: ~O(M^2)~.

Subtask 3: ~M \le 20000, N \le 100~

Nhận xét: Ta nhận thấy để có kết quả tối ưu thì cần xóa ~2~ hình chữ nhật sao cho tổng số lượng ô bằng ~1~ và bằng ~2~ nằm trong phần giao của cả ~2~ hình là nhiều nhất.

• Duyệt tất cả các ô ~a[i][j] = 2~, với mỗi ô tìm ~2~ hình chữ nhật bao phủ ô đó và cập nhật kết quả.

Độ phức tạp: ~O(N ^ 2 \cdot M)~.

Subtask 4: Ràng buộc gốc

• Trường hợp 1: 2 hình chữ nhật không giao nhau.

  • Kết quả là 2 hình chữ nhật có số ô bằng 1 nhiều nhất.

• Trường hợp 2: 2 hình chữ nhật giao nhau

  • Gọi ~sum[i][j]~ là tổng chỉ số các hcn phủ ô ~i, j~, ~sqr[i][j]~ là tổng bình phương chỉ số các hcn phủ ô ~i, j~.

  • Với mỗi ô ~a[i][j] = 2~ tìm 2 hình chữ nhật ~u, v~ phủ ô ~i, j~.

  • Khi đó ~sum[i][j] = u + v~ và ~sqr[i][j] = u^2 + v^2~. Từ đây có thể tìm ~u, v~ dễ dàng.

Độ phức tạp: ~O(N^2)~.

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 2e5 + 5, M = 5002;

int n, m, ps[3][M][M], cnt1[N], cnt2[M][M];

array<int, 4> rect[N];  

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

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

    cin >> n >> m;
    for(int i = 1; i <= n; i++){
        cnt1[i] = 0;

        int x1, x2, y1, y2; cin >> x1 >> x2 >> y1 >> y2;
        rect[i] = {x1, y1, x2, y2};
        int v[3] = {1, i, i * i};
        for(int j = 0; j < 3; j++){
            ps[j][x1][y1] += v[j], ps[j][x2 + 1][y2 + 1] += v[j];
            ps[j][x1][y2 + 1] -= v[j], ps[j][x2 + 1][y1] -= v[j];
        }
    }

    int cnt0 = 0, mx1 = 0, mx2 = 0;
    for(int i = 1; i <= m; i++){
        for(int j = 1; j <= m; j++){
            for(int k = 0; k < 3; k++){
                ps[k][i][j] += ps[k][i - 1][j] + ps[k][i][j - 1] - ps[k][i - 1][j - 1];
            }
            cnt2[i][j] += cnt2[i - 1][j] + cnt2[i][j - 1] - cnt2[i - 1][j - 1];
            if(ps[0][i][j] == 1) cnt1[ps[1][i][j]]++;
            else if(ps[0][i][j] == 0) cnt0++;
            else if(ps[0][i][j] == 2) cnt2[i][j]++;
        }
    }

    for(int i = 1; i <= n; i++){
        if(cnt1[i] > mx1) mx2 = mx1, mx1 = cnt1[i];
        else mx2 = max(mx2, cnt1[i]);
    }
    int ans = cnt0 + mx1 + mx2;

    for(int i = 1; i <= m; i++){
        for(int j = 1; j <= m; j++){
            if(ps[0][i][j] == 2){
                int x = ps[1][i][j];
                int y = ps[2][i][j];

                int a = 2;
                int b = -2 * x;
                int c = x * x - y;
                int delta = b * b - 4 * a * c;
                int u = (-b + sqrt(delta)) / (2 * a);
                int v = x - u;

                int x1 = max(rect[u][0], rect[v][0]), y1 = max(rect[u][1], rect[v][1]);
                int x2 = min(rect[u][2], rect[v][2]), y2 = min(rect[u][3], rect[v][3]);
                int ct2 = cnt2[x2][y2] - cnt2[x1 - 1][y2] - cnt2[x2][y1 - 1] + cnt2[x1 - 1][y1 - 1];
                ans = max(ans, cnt0 + cnt1[u] + cnt1[v] + ct2);
            }
        }
    }
    cout << ans << '\n';
}