[Simpson积分] BZOJ 1502 [NOI2005]月下柠檬树

%%% PoPoQQQ ：http://blog.youkuaiyun.com/popoqqq/article/details/39252719

题目大意：给定一棵由圆台和圆锥构成的柠檬树，月光以α的夹角平行射向地面，求阴影部分面积

求公切线

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define PI acos(-1.0)
#define EPS (1e-6)
using namespace std;

inline int dcmp(double a,double b){
	if (fabs(a-b)<1e-6) return 0;
	if (a<b) return -1;
	return 1;
}

struct Point{
	double x,y;
	Point(double x=0,double y=0):x(x),y(y) { }
};

struct Line{
	Point A,B;
	double k,b;
	Line(){
	}
	Line(Point iA,Point iB){
		A=iA; B=iB;
		calc();
	}
	void calc(){
		k=(A.y-B.y)/(A.x-B.x);  
        b=B.y-k*B.x; 
	}
	double f(double x){
		if (dcmp(x,A.x)<0 || dcmp(x,B.x)>0) return 0;
		return k*x+b;
	}
}L[505];
int tot;

struct Circle{
	double c,r;
	double f(double x){
		if (dcmp(c-x,r)>0) return 0;
		return sqrt(r*r-(c-x)*(c-x)); 
	}
}C[505];
int n;

inline double F(double x){
	double ret=0;
	for (int i=1;i<=n;i++)
		ret=max(ret,C[i].f(x));
	for (int i=1;i<=tot;i++)
		ret=max(ret,L[i].f(x));
	return ret;
}

inline double Simpson(double l,double mid,double r,double fl,double fm,double fr){
	double flm=F((l+mid)/2),fmr=F((mid+r)/2);
	double lret=(fl+4*flm+fm)*(mid-l)/6,rret=(fm+4*fmr+fr)*(r-mid)/6,ret=(fl+4*fm+fr)*(r-l)/6;
	if (fabs(lret+rret-ret)<EPS)
		return ret;
	else
		return Simpson(l,(l+mid)/2,mid,fl,flm,fm)+Simpson(mid,(mid+r)/2,r,fm,fmr,fr);
}

int main()
{
	double l,r,alpha;
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	scanf("%d%lf",&n,&alpha); 
	n++; alpha=1/tan(alpha);
	for (int i=1;i<=n;i++)
        scanf("%lf",&C[i].c),(C[i].c*=alpha)+=C[i-1].c;
    for (int i=1;i<n;i++)
    	scanf("%lf",&C[i].r);
    for(int i=1;i<=n;i++)  
    {  
    	l=min(l,C[i].c-C[i].r);  
        r=max(r,C[i].c+C[i].r);  
    }
    for (int i=2;i<=n;i++)
    {
    	double Len=C[i].c-C[i-1].c;  
        if(dcmp(Len,fabs(C[i-1].r-C[i].r))<=0)
            continue;
        double sin_alpha=(C[i-1].r-C[i].r)/Len;
        double cos_alpha=sqrt(1-sin_alpha*sin_alpha);  
        L[++tot]=Line(Point(C[i-1].c+C[i-1].r*sin_alpha,C[i-1].r*cos_alpha),Point(C[i].c+C[i].r*sin_alpha,C[i].r*cos_alpha));
	}
    printf("%.2lf\n",2*Simpson(l,(l+r)/2,r,F(l),F((l+r)/2),F(r)));
    return 0;
}