[BZOJ1502][NOI2005]月下柠檬树（辛普森积分）

=== ===

这里放传送门

=== ===

题解

据说这题正解好像是各种分类讨论一块块求面积？

然而一开始知道这题辛普森积分就可以过但是差一点还是写成了那种分类讨论的方法。。就一开始非常无脑求了公切线然后想它不就是会把曲线分成一段直线一段圆弧交替那样吗然后就一段段求不就行了吗但是忽略了可能有很多复杂的情况比如说有圆相交或内含的情况然后公切线就把圆割的乱七八糟一坨坨的。。。然后差点想各种分类讨论一下然后突然反应过来那我还用什么辛普森积分直接求面积不就可以了吗。。。

好吧其实正确的打开方式是直接把所有公切线都求出来然后算某个横坐标处的覆盖长度的时候就枚举所有公切线和圆求个交线然后做个线段覆盖就出来了。。注意公切线是线段所以要好好判断一下边界。

这里求公切线的方法纠结了好久。。一开始做数学做惯了上来就想用解析几何的方法联立圆方程然后想这不是有毒么然后就上网各种查怎么求圆的公切线然后终于发现了一种比较好用的方法。。。

因为这个题有一个方便的地方就是坐标啥的可以自己定，那我们就把圆心都定在x轴上。如下图

一个大点的圆和一个小点的圆，大圆的半径是$r_1$，小圆的半径是$r_2$。我们可以做一个辅助圆（右图蓝色），让它的半径是$r_1-r_2$。就相当于把大圆和小圆都往里“缩”了$r_2$那么多的一块，那么小圆就缩成一个点了，也就是小圆的圆心。这个时候从小圆的圆心往辅助圆做切线（右图浅蓝色），然后再把这个切线沿着和它垂直的方向（右图粉色）平移$r_2$那么多，就相当于还原回去，就求出了原来两个圆的公切线。

这里求过一个点的圆切线的方法用了一点解析几何。如果我们可以求出切点，那么就有了直线上的两个点——小圆的圆心是已知的，那么直线就有了。关键就是怎么求出切点。这里就用到了圆心在x轴上的便利性——可以发现切点弦一定是和x轴垂直的，那么也就是用一条垂直于x轴的直线去截那个辅助圆，如果知道了切点弦的方程，用勾股定理就可以方便地求出切点的纵坐标。而过圆$(x+a)^2+(y+b)^2=r^2$外一点$(x_0,y_0)$的两条切线形成的切点弦方程式$(x+a)(x_0+a)+(y+b)(y_0+b)=r^2$。这里$b$和$y_0$又是0，那也就是$(x+a)(x_0+a)=r^2$，化简一下就可以了。

这里平移直线的时候是先求出与当前直线垂直的单位法向量再乘以平移长度那样的，注意法向量的方向必须是“朝外”的。

还有一个问题就是内含的圆不会有贡献，直接判掉就可以了。不过因为是做线段覆盖所以不判也无所谓。。。

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const double eps=1e-7;
const double inf=1e9;
int n,cnt,tot;
double alpha,h[510],r[510],L,R,f1,f2,f3;
bool del[510];
struct Vector{
    double x,y;
    Vector (double X=0,double Y=0){x=X;y=Y;}
    Vector operator + (const Vector &a){return Vector(x+a.x,y+a.y);}
    Vector operator - (const Vector &a){return Vector(x-a.x,y-a.y);}
    double operator * (const Vector &a){return x*a.y-y*a.x;}
    Vector mul(double t){return Vector(x*t,y*t);}
    double len(){return sqrt(x*x+y*y);}
};
struct Line{
    Vector P,v;
    Line(){P=Vector();v=Vector();}
    void get(Vector A,Vector B){P=A;v=B-A;}
}t[510];
double getHeight(double A,double C){return sqrt(C*C-A*A);}
Vector GLI(Line A,Line B){
    Vector u=A.P-B.P;
    double t=(B.v*u)/(A.v*B.v);
    return A.P+A.v.mul(t);
}
void Trans(Vector &A,Vector &B,double dis){
    Vector v=B-A,w;
    double t;//平移直线
    w=Vector(-v.y,v.x);
    t=1/w.len();w=w.mul(t*dis);
    A=A+w;B=B+w;
}
void GetTan(int pos){
    double c1,c2,r1,r2,X,hgt,tr;
    Vector A,B,ip;
    c1=h[pos];c2=h[pos+1];
    r1=r[pos];r2=r[pos+1];
    if (r1-r2>eps){swap(c1,c2);swap(r1,r2);}
    tr=(r2-r1);//固定r1为半径较小的那个圆
    X=(tr*tr+c2*(c1-c2))/(c1-c2);//算出切点弦方程
    hgt=getHeight(fabs(X-c2),tr);//求出交点的纵坐标
    A=Vector(c1,0);B=Vector(X,hgt);
    if (A.x-B.x>eps) swap(A,B);
    Trans(A,B,r1);//平移直线
    t[++cnt].get(A,B);
}
bool In(int c1,int c2){
    double dis=fabs(h[c1]-h[c2]);
    return r[c2]-dis-r[c1]>-eps;
}
double F(double x){
    double nl,nr,Max=0;
    Line tmp;
    Vector ip;
    tmp.P=Vector(x,0);tmp.v=Vector(0,1);
    for (int i=1;i<=n;i++)
      if (del[i]==false){
          nl=h[i]-r[i];nr=h[i]+r[i];
          if (x-nl>eps&&nr-x>eps)
            Max=max(Max,getHeight(fabs(h[i]-x),r[i]));
      }
    for (int i=1;i<=cnt;i++){
        nl=t[i].P.x;nr=nl+t[i].v.x;
        if (x-nl>eps&&nr-x>eps){
            ip=GLI(tmp,t[i]);
            Max=max(Max,ip.y);
        }
    }
    return Max;
}
double calc(double fl,double fmid,double fr,double len){
    return len*(fl+4*fmid+fr)/6;
}
double Simpson(double l,double r,double fl,double fmid,double fr,double now){
    double mid=(l+r)/2,lmid=F((l+mid)/2),rmid=F((mid+r)/2),p,q;
    p=calc(fl,lmid,fmid,mid-l);
    q=calc(fmid,rmid,fr,r-mid);
    if (fabs(now-p-q)<eps) return now;
    return Simpson(l,mid,fl,lmid,fmid,p)+Simpson(mid,r,fmid,rmid,fr,q);
}
int main()
{
    scanf("%d%lf",&n,&alpha);
    for (int i=1;i<=n+1;i++){
        double p;scanf("%lf",&p);
        h[i]=h[i-1]+(p/tan(alpha));
    }
    for (int i=1;i<=n;i++) scanf("%lf",&r[i]);
    L=h[1]-r[1];R=h[n+1];
    for (int i=1;i<=n;i++){
        L=min(L,h[i]-r[i]);
        R=max(R,h[i]+r[i]);
    }
    for (int i=1;i<=n;i++)
      if (del[i]==false)
        for (int j=1;j<=n;j++)
          if (j!=i&&del[i]==false&&In(i,j))
            del[i]=true;//判掉内含的圆
    for (int i=1;i<=n;i++) GetTan(i);
    f1=F(L);f2=F((L+R)/2);f3=F(R);
    printf("%.2lf\n",Simpson(L,R,f1,f2,f3,calc(f1,f2,f3,R-L))*2);
    return 0;
}

偏偏在最后出现的补充说明

以前做辛普森积分的时候老是想把它当个函数求个函数值什么的。。其实如果做线段覆盖的话就非常的方便并且实际上$O(n)$的求也非常快。