BZOJ 1502 计算几何+自适应Simpson积分 解题报告

1502: [NOI2005]月下柠檬树

Time Limit: 5 Sec Memory Limit: 64 MB
Submit: 1205 Solved: 636
[Submit][Status][Discuss]
Description

Input

文件的第1行包含一个整数n和一个实数alpha，表示柠檬树的层数和月亮的光线与地面夹角(单位为弧度)。第2行包含n+1个实数h0,h1,h2,…,hn，表示树离地的高度和每层的高度。第3行包含n个实数r1,r2,…,rn，表示柠檬树每层下底面的圆的半径。上述输入文件中的数据，同一行相邻的两个数之间用一个空格分隔。输入的所有实数的小数点后可能包含1至10位有效数字。

Output

输出1个实数，表示树影的面积。四舍五入保留两位小数。

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

HINT

1≤n≤500，0.3

[解题报告]
首先，关于那个牛B轰轰的自适应Simpson积分
http://blog.csdn.net/greatwall1995/article/details/8639135
然后因为我其实也不理解一些细节，所以贴几个大神的题解
http://www.cnblogs.com/DaD3zZ-Beyonder/p/5676841.html
http://blog.csdn.net/qpswwww/article/details/44204425

代码如下：

/**************************************************************
    Problem: 1502
    User: onepointo
    Language: C++
    Result: Accepted
    Time:164 ms
    Memory:892 kb
****************************************************************/

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define max(a,b) (a>b)?a:b
#define min(a,b) (a<b)?a:b
#define N 505
#define EPS 1e-5  

struct Point
{  
    long double x,y;         
    Point(long double _ = .0,long double __ = .0):x(_),y(__) {}  
    Point operator +(const Point &a)const
    {return Point(x+a.x,y+a.y);}  
    Point operator -(const Point &a)const
    {return Point(x-a.x,y-a.y);}  
    Point operator *(double a)const    
    {return Point(x*a,y*a);}  
}points[N];
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[N];
struct Circle
{
    Point o;
    double r;
    Circle(){}
    Circle(Point _o,double _r):o(_o),r(_r){}
}circles[N];

int tot=0;//直线总数 
int n;
double alpha;

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;
}
inline 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()
{
    scanf("%d%lf",&n,&alpha);
    ++n;
    for(int i=1;i<=n;++i)
    {
        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);
    }
    double lowerBound=1e12,upperBound=-1e12;
    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;
}