bzoj 1502: [NOI2005]月下柠檬树

.

题意：一组与地面呈$\alpha$角的平行光线射向一颗由圆台组成的树，问地面投影面积。
题解：自适应辛普森积分法。
圆投影下来还是圆，母线投影下来就是公切线（如果上底和下底的投影重合，那么母线的投影就不用考虑了）。求面积可以用辛普森积分法：$\int_a^b f(x)dx\approx\frac{b-a}{6}\Bigl[f(a)+4f\bigl(\frac{a+b}{2}\bigr)+f(b)\Bigr]$这个不太精确，所以可以先用这个算出左右两边的面积，看看和是否等于整个一起算的，不等于就递归下去搞。至于怎么算$f(x)$,看代码吧。
代码：

#include<bits/stdc++.h>
using namespace std;
int n;
double alpha,h[510],r[510];
const double eps=1e-8;
struct circle
{
    double x,r;
}c[510];
struct line
{
    double k,b,l,r;
};
vector<line>l;
inline bool in(circle &x,circle &y)
{
    return abs(x.x-y.x)-abs(x.r-y.r)<eps;
}
double f(double x)
{
    double ans=0;
    for(int i=1;i<=n;i++)
    if(abs(c[i].x-x)-c[i].r<eps)
    ans=max(ans,sqrt(c[i].r*c[i].r-(c[i].x-x)*(c[i].x-x)));
    static int sz=l.size();
    for(int i=0;i<sz;i++)
    if(x-l[i].l>-eps&&x-l[i].r<eps)
    ans=max(ans,l[i].k*x+l[i].b);
    return ans*2;
}
inline double cal(double l,double r)
{
    return (r-l)/6*(f(l)+4*f((l+r)/2)+f(r));
}
double simpson(double l,double r)
{
    double md=(l+r)/2;
    double s=cal(l,r),ls=cal(l,md),rs=cal(md,r);
    if(abs(ls+rs-s)<eps)
    return s;
    else
    return simpson(l,md)+simpson(md,r);
}
int main()
{
    scanf("%d%lf",&n,&alpha);
    n++;
    for(int i=1;i<=n;i++)
    scanf("%lf",&h[i]);
    for(int i=1;i<n;i++)
    scanf("%lf",&r[i]);
    r[n]=0;
    for(int i=1;i<=n;i++)
    {
        h[i]+=h[i-1];
        c[i].r=r[i];
        c[i].x=h[i]/tan(alpha);
    }
    for(int i=1;i<n;i++)
    {
        if(in(c[i],c[i+1]))
        continue;
        double Sin=(c[i+1].r-c[i].r)/(c[i+1].x-c[i].x),Cos=sqrt(1-Sin*Sin);
        double x1=c[i].x-c[i].r*Sin,y1=c[i].r*Cos,x2=c[i+1].x-c[i+1].r*Sin,y2=c[i+1].r*Cos;
        line hh;
        hh.k=(y1-y2)/(x1-x2);
        hh.b=y1-hh.k*x1;
        hh.l=x1;
        hh.r=x2;
        if(hh.l-hh.r>eps)
        swap(hh.l,hh.r);
        l.push_back(hh);
    }/*
    puts("circles:");
    for(int i=1;i<=n;i++)
    printf("%lf %lf\n",c[i].x,c[i].r);
    puts("lines:");
    for(auto i:l)
    printf("%lf %lf %lf %lf\n",i.k,i.b,i.l,i.r);*/
    double ll=1e10,rr=-1e10;
    for(int i=1;i<=n;i++)
    {
        ll=min(ll,c[i].x-c[i].r);
        rr=max(rr,c[i].x+c[i].r);
    }
    printf("%.2lf\n",simpson(ll,rr));
}