关于自适应Simpson法的介绍可以去看我的另一篇blog

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

题意：空间里圆心在同一直线上且底面与地面平行的若干个圆台和顶层的圆锥以$\alpha$的角度投影到地面，求投影的面积。

---

（其实我是看po姐博客来的x）

首先把圆锥的顶点也看成一个半径为0的圆锥，对于每个高度为$h$的圆投影下去的坐标是$h/tan(\alpha)$，半径不变，而对于圆台的侧面投影下去是上下底两个圆的切线。

关于两个圆的切线可以像图上这样求（参考po姐博客的x）

$\sin(\alpha)=\frac{r_{i-1}-r_i}{x_i-x_{i-1}}$

这个式子同样可以适用于其他情况，然后三角函数搞一下求出其他点的坐标，根据坐标的范围用Simpson法来求面积。

只是把求点的函数值改成扫一遍所有的点和圆，找最大值。

哦然后我的Simpson递归的时候本来是让eps乘上$\frac{1}{2}$的…但是这样会T掉，不乘又会wa…orz

然后改成$\frac{2}{3}$就可以了…跑900+ms真神奇

#include
#include
#include
#include
using namespace std;
const int N=515;
struct point
{
double x,y;
};
struct line
{
point p1,p2;
double k,b;
line(){}
void modify(double x1,double y1,double x2,double y2)
{
p1.x=x1;p1.y=y1;
p2.x=x2;p2.y=y2;
k=(p1.y-p2.y)/(p1.x-p2.x);
b=p1.y-k*p1.x;
}
double f(double x)
{
return k*x+b;
}
}ls[N];
struct circle
{
double x,r;
}cs[N];
int n,tot;double alpha;
inline double f(double x)
{
double res=0;double eps=1e-6;
for(register int i=1;i<=n;i++) { double dis=fabs(x-cs[i].x); if(dis-cs[i].r>-eps)continue;
double y=sqrt(cs[i].r*cs[i].r-dis*dis);
res=max(res,y);
}
for(register int i=1;i<=tot;i++)
{
if(!(ls[i].p1.x<=x&&x<=ls[i].p2.x))continue;
double y=ls[i].f(x);
res=max(res,y);
}
return res;
}
inline double calc(double l,double r)
{
double mid=l+(r-l)/2;
return (f(l)+4.0*f(mid)+f(r))*(r-l)/6.0;
}
inline double asr(double l,double r,double area,double eps)
{
double mid=l+(r-l)/2;
double L=calc(l,mid),R=calc(mid,r);
if(fabs(L+R-area)<=eps*10.0)return L+R+(L+R-area)/10.0;
return asr(l,mid,L,eps*2.0/3.0)+asr(mid,r,R,eps*2.0/3.0);
}
inline double solve(double l,double r,double eps)
{
return asr(l,r,calc(l,r),eps);
}
int main()
{
//freopen("input.in","r",stdin);
double eps=1e-6,l,r;l=r=0;
scanf("%d%lf",&n,&alpha);alpha=1/(tan(alpha));
for(register int i=1;i<=n+1;i++)
{
scanf("%lf",&cs[i].x);
cs[i].x*=alpha;
cs[i].x+=cs[i-1].x;
}
for(register int i=1;i<=n;i++)
{
scanf("%lf",&cs[i].r);
}
for(register int i=1;i<=n+1;i++)
{
l=min(l,cs[i].x-cs[i].r);
r=max(r,cs[i].x+cs[i].r);
}
for(register int i=2;i<=n+1;i++)
{
double L=cs[i].x-cs[i-1].x;
if(L-fabs(cs[i].r-cs[i-1].r)<eps)continue;
//double sin_alpha=fabs(cs[i].r-cs[i-1].r)/L;
double sin_alpha=(cs[i-1].r-cs[i].r)/L;
double cos_alpha=sqrt(1-sin_alpha*sin_alpha);
ls[++tot].modify(cs[i-1].x+cs[i-1].r*sin_alpha,cs[i-1].r*cos_alpha,
cs[i].x+cs[i].r*sin_alpha,cs[i].r*cos_alpha);
}
printf("%.2lf",solve(l,r,eps)*2.0);
return 0;
}