[BZOJ 1502][NOI 2005]月下柠檬树(自适应Simpson积分)

题目链接

http://www.lydsy.com/JudgeOnline/problem.php?id=1502

思路

首先不会自适应Simpson积分求圆并(SPOJ CIRU)的请看我的上篇文章(http://blog.csdn.net/qpswwww/article/details/44201333)

这个题同样是可以用杀人越货之利器自适应Simpson积分做出来的，首先我们很容易可以发现不管光线是以何种角度投到地面上，光线永远是平行光线，投下来的圆的大小不会变，形状也不会变，变的只是圆与圆的圆心的距离，而且由于柠檬树是直的，所以很容易可以发现影子上所有圆心都在同一直线上，为了方便起见就定这条直线为$x$轴。

只要你有想象力的话，就很容易可以发现不管树是什么样的，影子永远是下面这样的形式：

就是一堆圆心在$x$轴上的大小不一的圆，圆心相邻的圆之间由两条公切线链接，公切线之间和所有的圆都是影子的一部分。显然这样的图形，扫描线被覆盖部分的长度的函数$f(x)$在这个图形的区间中是连续的，因此不必考虑将整个图形拆成若干个一坨一坨的图形再求积分，少了不少细节。

更令人愉悦的是，无论扫描线在何处，它被覆盖的部分也是永远是连续的，因此完全可以暴力找每个圆是否和扫描线有交，每条公切线段是否和扫描线有交，然后取扫描线被覆盖长度的最大值即可，这个的正确性也是很显然的。

确定每个圆的位置、求公切线都是很简单的事情，学过初中数学的都会做，渣渣我就不细讲了。

因此这个题比SPOJ的CIRU要简单很多，要考虑的情况也少很多，实际上只需要特判下某个圆是否被另一个圆完全盖住了，如果被完全盖住的话，就忽略掉这个圆。

代码

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 550
#define EPS 1e-6

using namespace std;

int n; //n个圆

struct Point
{
    double x,y;
    Point(){}
    Point(double _x,double _y):x(_x),y(_y){}
}points[MAXN];

struct Circle
{
    Point o;
    double r;
    Circle(){}
    Circle(Point _o,double _r):o(_o),r(_r){}
}circles[MAXN];

struct Line
{
    Point st,ed;
    double k,b;
    Line(){}
    Line(Point _st,Point _ed):st(_st),ed(_ed)
    {
        if(st.x>ed.x) swap(st,ed);
        k=(st.y-ed.y)/(st.x-ed.x);
        b=st.y-st.x*k;
    }
    double F(double x) //求直线上横坐标为x的点的纵坐标
    {
        return k*x+b;
    }
}lines[MAXN];
int tot=0; //直线总数

double f(double x)
{
    double ans=0; //扫描线被覆盖部分的长度
    for(int i=1;i<=n;i++) //枚举圆i与扫描线是否有交，若有交，更新被覆盖的长度
    {
        double dist=fabs(x-circles[i].o.x);
        if(dist-circles[i].r>-EPS)
            continue;
        double len=2*sqrt(circles[i].r*circles[i].r-dist*dist);
        ans=max(ans,len);
    }
    for(int i=1;i<=tot;i++) //枚举直线i是否与扫描线有交
        if(x>=lines[i].st.x&&x<=lines[i].ed.x)
            ans=max(ans,2*lines[i].F(x));
    return ans;
}

double calc(double len,double fL,double fM,double fR)
{
    return (fL+4*fM+fR)*len/6;
}

double Simpson(double L,double M,double R,double fL,double fM,double fR,double sqr)
{
    double M1=(L+M)/2,M2=(M+R)/2;
    double fM1=f(M1),fM2=f(M2);
    double g1=calc(M-L,fL,fM1,fM),g2=calc(R-M,fM,fM2,fR); //!!!!!!靠居然挂在了这里!!!!
    if(fabs(g1+g2-sqr)<EPS) return g1+g2;
    return Simpson(L,M1,M,fL,fM1,fM,g1)+Simpson(M,M2,R,fM,fM2,fR,g2);
}

int main()
{
    double lowerBound=1e12,upperBound=-1e12;
    double alpha; //光线的角度
    scanf("%d%lf",&n,&alpha);
    n++;
    for(int i=1;i<=n;i++) //注意第n个圆是个点！！！
    {
        double h;
        scanf("%lf",&h); //注意是每层的高度，不要被坑了
        circles[i].o.y=0;
        circles[i].o.x=h/tan(alpha);
        circles[i].o.x+=circles[i-1].o.x;
    }
    for(int i=1;i<n;i++)
        scanf("%lf",&circles[i].r);
    for(int i=1;i<=n;i++)
    {
        lowerBound=min(lowerBound,circles[i].o.x-circles[i].r);
        upperBound=max(upperBound,circles[i].o.x+circles[i].r);
    }
    for(int i=1;i<n;i++)
    {
        double dist=circles[i+1].o.x-circles[i].o.x;
        if(dist-fabs(circles[i+1].r-circles[i].r)<-EPS) continue; //特判：一个圆把另一个圆完全盖住了
        double sinAlpha=(circles[i].r-circles[i+1].r)/dist;
        double cosAlpha=sqrt(1-sinAlpha*sinAlpha);
        lines[++tot]=Line(Point(circles[i].o.x+circles[i].r*sinAlpha,circles[i].r*cosAlpha),Point(circles[i+1].o.x+circles[i+1].r*sinAlpha,circles[i+1].r*cosAlpha));
    }
    double L=lowerBound,R=upperBound;
    double M=(L+R)/2;
    double fL=f(L),fM=f(M),fR=f(R);
    printf("%.2lf\n",Simpson(L,M,R,fL,fM,fR,calc(R-L,fL,fM,fR)));
    return 0;
}