题解：P13456 [GCJ 2008 Qualification] Fly Swatter

· · 题解

题意

三维转二维了，这个转化不难。

你有一个圆环，大圆半径为 R，小圆半径为 R-t，圆内有若干有厚度的弦，弦宽度的一半为 r，两根弦间的距离为 g。

现在有一个半径为 f 的圆，中心从大圆中均匀随机选取。求这个圆与圆环或弦相交的概率。

题解

首先把 f 去掉，即圆转点。那么就让 r \leftarrow r+f,t \leftarrow t+f,g \leftarrow g-2f，反正就是判圆相交只和圆心距和圆半径之和有关，那么不如把这个 +f 的部分转到另外一个圆上。

问题转为求圆环和弦的组合图形的面积与大圆的面积的比值。而圆环的面积是好求的，于是问题转为求小圆内弦的面积。

我们可以用容斥：算出横向的弦总面积和纵向的弦总面积，再减去横向弦与纵向弦交叉面积。

易得半径为 r 的圆上一条距圆心 d 的弦会把圆分成两部分，其中一部分是 d\sqrt{r^2-d^2}-\arcsin (\frac{d}{r})r^2。记 f_0(r,d) 用于计算此部分的面积。那么半径为 r 的圆上一条中心线距圆心分别为 a,b 的弦间的距离就是 f_1(r,a,b)=f_0(r,b)-f_0(r,a)。注意边界情况即可。那么横向弦和纵向弦总面积就很好计算，为 2\sum_{i=-k}^{k}f_1(R-t,is-r,is+r)。

现在的问题是：已知圆内一条横向弦和一条纵向弦，求相交部分的面积。这部分我用几何法没有求出来，只好上数值积分了。

将横向弦和纵向弦视作两组平行线，相交部分是一个矩形。再考虑到弦在园内的限制，相交部分就是一个矩形和圆的重叠面积。那么我们把圆心当做 (0,0)，矩形左下和右上坐标分别是 (x_1,y_1),(x_2,y_2)，这个积分是

\int_{x_1}^{x_2}\max(0,\min(y_2,\sqrt{(R-t)^2-x^2})-\max(y_1,-\sqrt{(R-t)^2-x^2}))dx

用自适应辛普森法求解。

由于题目中保证横向弦和纵向弦各不超过 1000 根，因此暴力统计重叠部分的复杂度是正确的。

实际写的时候注意精度问题，在开根和求反正弦时 double 的精度可能会不够，还有自适应辛普森的 \epsilon 也要开小一些，但是全 long double 会超时，整个卡常还是挺难受的。

代码如下：

#include<iostream>
#include<vector>
#include<cmath>
#include<cassert>
#include<algorithm>
using namespace std;
template<typename T,typename D=double>
D _simpson(T f,D l,D r){
    auto mid=(l+r)/2;
    return (f(l)+4*f(mid)+f(r))*(r-l)/6;
}
template<typename T,typename D=double>
D simpson(T f,D l,D r,D eps,D ans){
    auto mid=(l+r)/2,l_=_simpson(f,l,mid),r_=_simpson(f,mid,r);
    if(fabs(l_+r_-ans)<=15*eps) return l_+r_+(l_+r_-ans)/15;
    else    return simpson(f,l,mid,eps/2,l_)+simpson(f,mid,r,eps/2,r_);
}
template<typename T,typename D=double>
D simpson(T f,D l,D r,D eps=1e-8){
    return simpson(f,l,r,eps,_simpson(f,l,r));
}
long double gr,gb,gd;
double f(long double x){
    long double y=sqrtl((0.0l+gr-x)*(gr+x));
    return max(0.0,(double)(min(gd,y)-max(gb,-y)));
}
constexpr long double pie=3.1415926535897932384626433832795l;
double f0(double r0,double d){
    if(d>r0)    d=r0;
    if(d<-r0)   d=-r0;
    return d*sqrtl((0.0l+r0-d)*(r0+d))+asinl((long double)(d)/r0)*r0*r0;
}
double f1(double r0,double a,double b){
    return max(0.0,f0(r0,b)-f0(r0,a));
}
double cs[4001];
double calc(double r0,double r,double g){
    double s=2*r+g,ans=0;
    int cnt=int((r0+r)/s);
    double eps=1e-6/(1ll*cnt*cnt);
    for(int i=-cnt;i<=cnt;i++){
        cs[i+cnt]=i*s;
        ans+=f1(r0,cs[i+cnt]-r,cs[i+cnt]+r);
    }
    ans*=2;
    gr=r0;
    for(int i=-cnt;i<=cnt;i++){
        gb=cs[i+cnt]-r;
        gd=cs[i+cnt]+r;
        for(int j=-cnt;j<=cnt;j++){
            ans-=simpson(f,max(cs[j+cnt]-r,-r0),min(cs[j+cnt]+r,r0),eps);
        }
    }
    return ans;
}
int main(){
    cout.flags(ios::fixed);
    cout.precision(8);
    ios::sync_with_stdio(false);
    cin.tie(0);
    int T;
    cin >> T;
    for(int _=1;_<=T;_++){
        double f,R,t,r,g;
        cin >> f >> R >> t >> r >> g;
        r+=f;
        t+=f;
        g-=2*f;
        if(g<0 || abs(g)<1e-6){
            cout << "Case #" << _ << ": ";
            cout << "1.0\n";
            continue;
        }
        double s3=calc(R-t,r,g);
        cout << "Case #" << _ << ": ";
        cout << 1.0-(pie*(R-t)*(R-t)-s3)/(pie*R*R) << '\n';
//      cout << "s1=" << s1 << ", s2=" << s2 << ", s3=" << s3 << ".\n";
    }
    return 0;
}