AOJ2256 ケーキ分割問題 - ferinの競プロ帳

問題ページ
直線 $x=0$ 上の点 $p=(0,p_y)$ を決めたときに， $x=w$ 上のどの位置に点を置けば2等分できるか？ $p$ を原点として $2N$ 点を偏角でソートする．このとき $N$ 番目と $N+1$ 番目の点の間を通るような線分を引くと，2等分することができる． $p$ と $N$ 番目の点を結ぶ線分が $x=w$ のときに取る $y$ 座標から $p$ と $N+1$ 番目の点を結ぶ線分が $x=w$ のときに取る $y$ 座標までの区間に点を置けばよい．

点 $p$ を $0 \sim H$ まで動かしたときに $N,N+1$ 番目の点が変化するタイミングがいつか考える．2点 $p_1,p_2$ を通る線分が $N$ 個ずつに2等分する場合，この線分が $x=0$ のときの $y$ 座標で切り替わる．このような点は $O(N)$ 個存在する．

これらの点で切り分けたときの $O(N)$ 個の区間が答えに寄与する分を計算する．点 $p$ に対して $x=w$ 上の点が取る区間の長さは一次関数になっている．(cf. AOJ : 2256 Divide the Cake ケーキ分割問題 - 衛生のゴミ箱(二勝目)) １次関数 $f(x)$ の積分は台形公式より $\int_{a}^{b} f(x) dx = (b-a) \times f((a+b)/2)$ である． $f((a+b)/2)$ を計算するには，区間の中間の点に対して $x=w$ 上の点が取る区間の長さを計算すればよい．これは $O(N^2)$ で求められるので，合計 $O(N^3)$ で答えを計算できる．

1次関数の積分は台形公式 $\int_{a}^{b} f(x) dx = (b-a) \times f((a+b)/2)$
3次以下の関数の積分はシンプソン公式 $\int_{a}^{b} f(x) dx = \frac{b-a}{6} \{f(a)+f(b)+4f((a+b)/2)\}$

$o(N^3)$ の解ないのかな

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG_   
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
const ll INF = 1LL<<60;  
  
using R = long double;  
const R EPS = 1e-10;  
int sgn(const R &r) { return (r>EPS) - (r<-EPS); }  
int sgn(const R &a, const R &b) { return sgn(a-b); }  
  
void solve() {  
    ll w, h, n;  
    cin >> w >> h >> n;  
    if(w == 0) exit(0);  
    using P = pair<R,R>;  
    vector<P> ps;  
    vector<ll> x(2*n+4), y(2*n+4);  
    REP(i, 2*n) {  
        cin >> x[i] >> y[i];  
        ps.push_back({x[i], y[i]});  
    }  
    x[2*n] = 0, y[2*n] = 0, ps.push_back({0, 0});  
    x[2*n+1] = 0, y[2*n+1] = h, ps.push_back({0, h});  
    x[2*n+2] = w, y[2*n+2] = 0, ps.push_back({w, 0});  
    x[2*n+3] = w, y[2*n+3] = h, ps.push_back({w, h});  
  
    vector<R> v;  
    v.push_back(0);  
    v.push_back(h);  
    REP(i, ps.size()) FOR(j, i+1, ps.size()) {  
        if(sgn(x[i],x[j]) == 0) continue;  
        // ((x[i],y[i]), (x[j],y[j])) で分断したときN個ずつになるか  
        ll cnt = 0, tc = 0;  
        R a = ((R)y[i]-y[j])/((R)x[i]-x[j]), b = y[i] - a*x[i];  
        REP(k, 2*n) {  
            if(sgn(a*x[k]+b, y[k]) < 0) cnt++;  
            else if(sgn(a*x[k]+b, y[k]) == 0) tc++;  
        }  
        while(cnt < n && tc > 0) cnt++, tc--;  
        if(cnt == n) v.push_back(b);  
    }  
    sort(ALL(v));  
    v.erase(unique(ALL(v)), v.end());  
  
    R ans = 0;  
    REP(i, (ll)v.size()-1) {  
        // [y1,y2] の分を計算する  
        R y1 = v[i], y2 = v[i+1];  
        if(sgn(y1) < 0 || sgn(y2,h) > 0) continue;  
  
        R midy = (y1+y2)/2.0;  
        R miny = h, maxy = 0;  
        for(auto p: ps) {  
            if(sgn(p.first) <= 0) continue;  
            // ((0,midy), p) で分断してN個ずつになるか  
            ll cnt = 0, tc = 0;  
            R a = (p.second-midy)/p.first, b = midy;  
            REP(k, 2*n) {  
                if(sgn(a*x[k]+b, y[k]) < 0) cnt++;  
                else if(sgn(a*x[k]+b, y[k]) == 0) tc++;  
            }  
            while(cnt < n && tc > 0) cnt++, tc--;  
            if(cnt == n) {  
                R ty = a*w+b;  
                chmin(miny, ty);  
                chmax(maxy, ty);  
            }  
        }  
        chmax(miny, 0.0l);  
        chmin(maxy, (R)h);  
  
        ans += (y2-y1) * max(0.0l, maxy-miny);  
    }  
  
    cout << fixed << setprecision(10) << ans/(h*h) << endl;  
}  
  
int main(void) {  
    while(1) solve();  
  
    return 0;  
}