jzoj6678 苏菲的世界 (辛普森积分，圆与圆并)

积分，圆与圆并)"/>

jzoj6678 苏菲的世界 (辛普森积分，圆与圆并)

题意

20个球，求面积并

做法

辛普森积分+圆并。
辛普森积分：知道一个函数在 l , r , 0.5 ( l + r ) l,r,0.5(l+r) l,r,0.5(l+r)处的值，可以估计其区间积分为 ( r − l ) ( f ( l ) + f ( r ) + 4 f ( 0.5 ( l + r ) ) ) 6 (r-l)\frac {(f(l)+f(r)+4f(0.5(l+r)))} 6 (r−l)6(f(l)+f(r)+4f(0.5(l+r)))​。二分处理，直到 r − l < = g a p r-l<=gap r−l<=gap且当前区间的估计值与子区间的估计值之差小于eps为止。
圆并做法：先去重合，去包含，然后考虑上轮廓减下轮廓，枚举每一个圆，用一些几何求出所有被覆盖的区间，然后再对没有被覆盖的区间做积分。

#include <bits/stdc++.h>
using namespace std;
typedef double db;
const db gap = 0.1, pi = acos(-1), eps = 1e-6;
const int N = 30;
int n, x[N], y[N], z[N], r[N];
struct pot{db x,y;pot(db _x=0, db _y=0) {x=_x,y=_y;}pot operator-(pot b){return pot(x-b.x,y-b.y);}pot operator+(pot b){return pot(x+b.x,y+b.y);}db operator*(pot b){return x*b.y-y*b.x;}pot operator*(db r){return pot(x*r,y*r);}
} o[N];
db R[N];
int m, del[N];
#define sqr(x) ((x)*(x))
#define zero(x) (abs(x)<=eps)
db len2(pot a) {return sqr(a.x)+sqr(a.y);
}pair<db,int> tp[N * 2];
int cc = 0;
void addseg(db l, db r) {tp[++cc]=make_pair(l,1);tp[++cc]=make_pair(r,-1);
}db ciruni() {db ret = 0;for(int i = 1; i <= m; i++) {cc = 0;tp[++cc]=make_pair(0,0);tp[++cc]=make_pair(-pi,0);tp[++cc]=make_pair(pi,0);for(int j = 1; j <= m; j++) if (i != j) {db d2 = len2(o[i] - o[j]);db d = sqrt(d2);if (d >= R[i] + R[j]) continue;db cs = (sqr(R[i]) + d2 - sqr(R[j])) * 0.5 / (R[i] * d);db ang = acos(cs);db b = atan2(o[j].y - o[i].y, o[j].x - o[i].x);db l = b - ang, r = b + ang;if (-pi <= l && r <= pi)addseg(l,r);else if (r > pi) {addseg(l,pi), addseg(-pi,r-2*pi);} else if (l < -pi) {addseg(-pi,r),addseg(l+2*pi,pi);}}sort(tp+1,tp+1+cc);int js=0;for(int j=1;j<cc;j++){js += tp[j].second;if(!zero(tp[j].first - tp[j + 1].first) && js == 0) {db jl = tp[j].first, jr = tp[j + 1].first;pot l = o[i] + pot(cos(jl), sin(jl)) * R[i];pot r = o[i] + pot(cos(jr), sin(jr)) * R[i];db w = jr - jl;ret -= (l.y + r.y) * 0.5 * (r.x - l.x);ret += w * R[i] * R[i] * 0.5 - R[i] * R[i] * sin(w) * 0.5;}}}return ret;
}db calc_z(db tz) {m=0;db xmx=-1000,xmi=1000;for(int i = 1; i <= n; i++) {if (z[i] - r[i] <= tz && tz <= z[i] + r[i]) {++m;o[m]=pot(x[i],y[i]);R[m]=sqrt(r[i]*r[i]-sqr(tz-z[i]));xmx=max(xmx,x[i]+R[m]);xmi=min(xmi,x[i]-R[m]);}}if (m==0)return 0;if (m==1)return pi*R[1]*R[1];memset(del,0,sizeof del);for(int i = 1; i <= m; i++) {for(int j = i + 1; j <= m; j++) {if (zero(o[i].x-o[j].x) && zero(o[i].y-o[j].y) && zero(R[i]-R[j])) {del[j] = 1;}}}for(int i=1;i<=m;i++){for(int j=1;j<=m;j++)if(j!=i&&R[i]>R[j]){if (len2(o[i]-o[j])<=sqr(R[i]-R[j]))del[j]=1;}}int tm=0;for(int i=1;i<=m;i++)if(!del[i]){o[++tm]=o[i];R[tm]=R[i];}m=tm;db v = ciruni();
//	printf("%.5lf %.5lf\n",tz,v);return v;
}db simp_z(db l, db r, db fl, db fr, db &fmid) {return (r-l)/6*(fl+fr+4*(fmid=calc_z((l+r)*0.5)));
}db sum_z(db lz, db rz, db s, db fl, db fr, db fmid) {db mid = 0.5 * (lz + rz), lmid = 0, rmid = 0;db ansl = simp_z(lz, mid, fl, fmid, lmid);db ansr = simp_z(mid, rz, fmid, fr, rmid);if (rz - lz > gap || abs(ansl + ansr - s) >= eps) {return sum_z(lz, mid, ansl, fl, fmid, lmid) + sum_z(mid, rz, ansr, fmid, fr, rmid);} else return ansl + ansr;
}db start_z(db lz, db rz) {db fl = calc_z(lz), fr = calc_z(rz), fmid = 0, z = simp_z(lz, rz, fl, fr, fmid);return sum_z(lz, rz, z, fl, fr, fmid);
}int main() {freopen("sphere.in","r",stdin);freopen("sphere.out","w",stdout);cin >> n;int zmx=-1000,zmi=1000;for(int i=1;i<=n;i++){scanf("%d %d %d %d",&x[i],&y[i],&z[i],&r[i]);zmx=max(zmx,z[i]+r[i]);zmi=min(zmi,z[i]-r[i]);}db ans = start_z(zmi, zmx);printf("%.3lf\n",ans);
}