您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> C++ >> 關於C++

BZOJ 2178 圓的面積並 - Simpson 積分法求解

編輯：關於C++

問題描述

給定 N (≤ 1000) 個圓，每個圓由圓心坐標 (x, y) 和半徑 r 定義，座標和半徑均為絕對值 < 1000 的整數。求所有圓覆蓋區域的面積並，結果保留三位小數。

核心思路：Simpson 積分 + 掃描線

將二維面積問題轉化為一維積分：對 x 軸進行積分，被積函數 f(x) 為掃描線 x = 常數與所有圓相交形成的線段長度總和（即該 x 位置的「掃描線長度」）。

面積 = ∫ f(x) dx，其中 f(x) = 所有圓在 x 位置的投影線段長度的並集長度。

關鍵子問題：如何計算 f(x)

對於固定的 x，每個圓 (cx, cy, r) 若 |x - cx| ≤ r，則與掃描線相交，形成線段：

[cy - √(r² - (x-cx)²),  cy + √(r² - (x-cx)²)]

將所有相交圓產生的線段按左端點排序，掃描合併重疊區間，累加總長度即得 f(x)。

Simpson 積分公式

自適應 Simpson 積分遞歸公式：

S(l, r) = (f(l) + 4f(mid) + f(r)) * (r-l) / 6
若 |S(l, r) - S(l, mid) - S(mid, r)| < eps，接受結果
否則遞歸計算左右兩半區間

代碼框架（關鍵部分）

const double EPS = 1e-7;
struct Circle { double x, y, r; } c[1005];

// 計算掃描線 x 處的覆蓋長度
double calc_len(double x) {
    vector<pair<double,double>> segs;
    for (int i = 0; i < n; ++i) {
        double dx = x - c[i].x;
        if (fabs(dx) <= c[i].r) {
            double h = sqrt(c[i].r * c[i].r - dx * dx);
            segs.emplace_back(c[i].y - h, c[i].y + h);
        }
    }
    if (segs.empty()) return 0;
    sort(segs.begin(), segs.end());
    double len = 0, curL = segs[0].first, curR = segs[0].second;
    for (auto &s : segs) {
        if (s.first > curR) {
            len += curR - curL;
            curL = s.first; curR = s.second;
        } else curR = max(curR, s.second);
    }
    len += curR - curL;
    return len;
}

// Simpson 遞歸
double simpson(double l, double r, double fl, double fm, double fr) {
    return (fl + 4*fm + fr) * (r-l) / 6;
}

void asr(double l, double r, double fl, double fm, double fr, double whole) {
    double mid = (l+r)/2;
    double flm = calc_len((l+mid)/2), frm = calc_len((mid+r)/2);
    double left = simpson(l, mid, fl, flm, fm);
    double right = simpson(mid, r, fm, frm, fr);
    if (fabs(left + right - whole) <= 15*EPS) {
        ans += left + right;
        return;
    }
    asr(l, mid, fl, flm, fm, left);
    asr(mid, r, fm, frm, fr, right);
}

int main() {
    // 讀入 n 和圓參數
    // 確定積分區間 [xmin, xmax] 為所有圓的 x-r 到 x+r
    double fl = calc_len(xmin), fm = calc_len((xmin+xmax)/2), fr = calc_len(xmax);
    asr(xmin, xmax, fl, fm, fr, simpson(xmin, xmax, fl, fm, fr));
    printf("%.3f\n", ans);
}