BZOJ 1502([NOI2005]月上柠檬树-Simpson积分)

1502: [NOI2005]月下柠檬树李哲非常非常喜欢柠檬树,特别是在静静的夜晚,当天空中有一弯明月温柔地照亮地面上的景物时,他必会悠闲地坐在他亲手植下的那棵柠檬树旁,独自思索着人生的哲理。
李哲是一个喜爱思考的孩子,当他看到在月光的照射下柠檬树投在地面上的影子是如此的清晰,马上想到了一个问题:树影的面积是多大呢?
李哲知道,直接测量面积是很难的,他想用几何的方法算,因为他对这棵柠檬树的形状了解得非常清楚,而且想好了简化的方法。
李哲将整棵柠檬树分成了 n 层,由下向上依次将层编号为 1,2,...,n。从第 1到 n-1 层,每层都是一个圆台型,第 n 层(最上面一层)是圆锥型。对于圆台型,其上下底面都是水平的圆。对于相邻的两个圆台,上层的下底面和下层的上底面重合。第 n 层(最上面一层)圆锥的底面就是第 n-1 层圆台的上底面。所有的底面的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为h1,h2,...,hn,第 1 层圆台的下底面距地面的高度为 h0,以及每层的下底面的圆的半径 r1,r2,...,rn。李哲用熟知的方法测出了月亮的光线与地面的夹角为 alpha。

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

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

Sample Input2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output171.97

HINT

1≤n≤500，0.3<alpha<π 2，0<hi≤100，0<ri≤100。
 10%的数据中，n=1。
30%的数据中，n≤2。
60%的数据中，n≤20。
100%的数据中，n≤500。

这题是Simpson模版题。有一个关键是圆的公切线的求法。首先延长圆的切线求l，再求θ，最后用相似三角形。

#include<cstdio>#include<cstdlib>#include<cstring>#include<cmath>#include<iostream>#include<algorithm>#include<functional>using namespace std;#define MAXN (500+10)#define eps 1e-8double h[MAXN],s[MAXN],r[MAXN],alpha,l1=0,r1=0;double sqr(double x){return x*x;}int n,size=0;struct S{double x,y,p,q;}c[MAXN];double f(double l){double t=0.0;for (int i=0;i<n;i++){if (fabs(s[i]-l)<r[i]) t=max(t,sqrt(sqr(r[i])-sqr(s[i]-l)));}for (int i=1;i<=size;i++){if (c[i].x<l&&l<c[i].p)t=max(t,c[i].y+(c[i].q-c[i].y)*(l-c[i].x)/(c[i].p-c[i].x));}return t;}double simpson(double l,double r,double fl,double fmid,double fr){return (fl+4*fmid+fr)*(r-l)/6;}double rsimpson(double l,double r,double fl,double fmid,double fr){double m=(l+r)/2;double p=f((l+m)/2),q=f((m+r)/2);double x=simpson(l,r,fl,fmid,fr),y=simpson(l,m,fl,p,fmid),z=simpson(m,r,fmid,q,fr);if (fabs(x-y-z)<eps)return y+z;return rsimpson(l,m,fl,p,fmid)+rsimpson(m,r,fmid,q,fr);}int main(){scanf("%d%lf",&n,&alpha);alpha=1/tan(alpha);for (int i=0;i<=n;i++){scanf("%lf",&h[i]);if (i) h[i]+=h[i-1];s[i]=h[i]*alpha;}for (int i=0;i<n;i++){scanf("%lf",&r[i]);l1=min(l1,s[i]-r[i]);r1=max(r1,s[i]+r[i]);}r[n]=0;for (int i=0;i<n;i++){double d=s[i+1]-s[i];if (d>fabs(r[i]-r[i+1])){c[++size].x=s[i]-r[i]*(r[i+1]-r[i])/d;c[size].y=sqrt(sqr(r[i])-sqr(c[size].x-s[i]));c[size].p=s[i+1]-r[i+1]*(r[i+1]-r[i])/d;c[size].q=sqrt(sqr(r[i+1])-sqr(c[size].p-s[i+1]));}}r1=max(r1,s[n]);printf("%.2lf",2*rsimpson(l1,r1,0,f((l1+r1)/2),0));return 0;}