Shadow of Survival

Shadow of Survival

题目链接：http://dutacm.club:7217/codesheaven/problem.php?id=1083

题目大意：给出一个由圆台构成的塔和一个点光源，问塔的阴影面积.

Simpson积分

注意到圆台底面由点光源投影后仍为圆（对应弦成正比），于是就与这道题类似http://www.cnblogs.com/barrier/p/6512826.html.不再赘述，经测试eps取1e-5时最优.

/*需要注意的是，#define Abs(x) (x>0?x:-x) 这种写法可能会因为x为复合表达式而导致运算错误*/

代码如下：

#include <cstdio>
#include <cmath>
#include <iostream>
#define N 1005
using namespace std;
const double eps=1e-5;
const double inf=1e12;
double d,z;
int n,k;
int dcmp(double x){
    if(fabs(x)<eps)return 0;
    else return x<0?-1:1;
};
struct point{
    double x,y;
    point(double X=0,double Y=0){x=X;y=Y;}
};
struct circle{
    point o;
    double r;
    circle(point O=point(0,0),double R=0){o=O;r=R;}
    bool within(circle tmp){
        return fabs(o.x-tmp.o.x)-fabs(r-tmp.r)<=0;
    }
    bool xin(double x){
        return o.x-r<=x&&x<=o.x+r;
    }
}c[N];
struct line{
    point s,t;
    double k,b;
    line(point S=point(-1,0),point T=point(1,0)){
        s=S;t=T;
        if(s.x>t.x)swap(s,t);
        k=(s.y-t.y)/(s.x-t.x);
        b=s.y-k*s.x;
    }
    bool xin(double x){
        return s.x<=x&&x<=t.x;
    }
    double f(double x){
        return k*x+b;
    }
}l[N];
void init(){
    scanf("%lf",&c[0].r);
    c[0].o=point(0,0);
    k=0;
    for(int i=1;i<n;++i){
        double h,r;
        scanf("%lf%lf",&h,&r);
        c[i]=circle(point(h*d/(z-h),0),z*r/(z-h));
    }
    for(int i=1;i<n;++i){
        if(c[i].within(c[i-1]))continue;
        double sin_theta=(c[i].r-c[i-1].r)/(c[i].o.x-c[i-1].o.x);
        double cos_theta=sqrt(1-sin_theta*sin_theta);
        l[k++]=line(point(c[i-1].o.x-c[i-1].r*sin_theta,c[i-1].r*cos_theta),
                    point(c[i].o.x-c[i].r*sin_theta,c[i].r*cos_theta));
    }
}
double F(double x){
    double len=0;
    for(int i=0;i<n;++i){
        if(c[i].xin(x)){
            double d=c[i].o.x-x;
            len=max(len,2*sqrt(c[i].r*c[i].r-d*d));
        }
    }
    for(int i=0;i<k;++i)
        if(l[i].xin(x))
            len=max(len,2*l[i].f(x));
    return len;
}
double cal(double l,double r){
    return (F(l)+4*F((l+r)/2)+F(r))/6*(r-l);
}
double simpson(double l,double r,double s){
    double mid=(l+r)/2;
    double x=cal(l,mid),y=cal(mid,r);
    if(!dcmp(x+y-s))return s;
    else return simpson(l,mid,x)+simpson(mid,r,y);
}
void solve(){
    double L=inf,R=-inf;
    for(int i=0;i<n;++i){
        L=min(L,c[i].o.x-c[i].r);
        R=max(R,c[i].o.x+c[i].r);
    }
    printf("%.2lf
",simpson(L,R,cal(L,R)));
}
int main(void){
    while(~scanf("%lf%lf%d",&d,&z,&n)){
        init();
        solve();
    }
}