【bzoj 1502】月下柠檬树

月下柠檬树

题意

求n个圆与他们的公切线的定积分。

解法

求出圆的公切线就可以了。
特别坑的一点 ：
最两端的圆，有可能会被其他的圆所包含，所以要重新求一下最左端与最右端。
比较坑的一点 ：
精度要设小一点，不然会TLE。

代码如下：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <cctype>
#define INF 2139062143
#define MAX 0x7ffffffffffffff
#define del(a,b) memset(a,b,sizeof(a))
using namespace std;
typedef long long ll;
template<typename T>
inline void read(T&x)
{
    x=0;T k=1;char c=getchar();
    while(!isdigit(c)){if(c=='-')k=-1;c=getchar();}
    while(isdigit(c)){x=x*10+c-'0';c=getchar();}x*=k;
}
const int maxn=(500+5)*2;
struct cir{
    double x,r;
    double h(double k){
        double d=fabs(k-x);
        return sqrt(r*r - d*d);
    }
}c[maxn];
struct lin{
    // y=kx+b
    double k,b,lx,ly;
    double h(double x){
        return x*k+b;
    }
}p[maxn];
int n;
double F(double x){
    double ans=0.0;
    for(int i=1;i<=n;i++){
        if(x>(c[i].x-c[i].r)&&x<(c[i].x+c[i].r))
            ans=max(ans,c[i].h(x));
    }
    for(int i=1;i<n;i++){
        if(x>=p[i].lx&&x<=p[i].ly)
            ans=max(ans,p[i].h(x));
    }
    return ans;
}
double simpson(double a,double b){
    double mid=(a+b)/2;//*0.5???
    return ( F(a) + 4*F(mid) + F(b) ) * (b-a)/6;
}
double asr(double a,double b,double eps,double A){
    double mid = (a+b)/2;
    double L = simpson(a,mid) , R = simpson(mid,b);
    if(fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15;
    return asr(a,mid,eps/2,L) + asr(mid,b,eps/2,R);
}
double asr(double l,double r,double eps){
    return asr(l,r,eps,simpson(l,r));
}
double alpha;
double h[maxn];
int main()
{
    scanf("%d%lf",&n,&alpha);n++;
    alpha= 1.0 / tan(alpha);
    for(int i=1;i<=n;i++) scanf("%lf",&h[i]);
    for(int i=1;i<n;i++) scanf("%lf",&c[i].r);c[n].r=0.0;
    for(int i=1;i<=n;i++) c[i].x=c[i-1].x+h[i];
    for(int i=1;i<=n;i++) c[i].x*=alpha;
    
    for(int i=1;i<n;i++){
        if((c[i+1].x+c[i+1].r)<=(c[i].x+c[i].r)) continue;
        if(fabs(c[i].r-c[i+1].r)<=1e-6){
            p[i].lx=c[i].x,p[i].ly=c[i+1].x;
            p[i].k=0.0;p[i].b=c[i].r;
            continue;
        }
        if(c[i].r<c[i+1].r){
            double _cos=(c[i+1].r-c[i].r)/(c[i+1].x-c[i].x);
            double d=_cos*c[i+1].r;
            double y2 = sqrt(c[i+1].r*c[i+1].r - d*d);
            p[i].ly=c[i+1].x-d;
            d=_cos*c[i].r;
            double y1 = sqrt(c[i].r*c[i].r - d*d);
            p[i].lx=c[i].x-d;
            p[i].k = (y1-y2)/(p[i].lx-p[i].ly);
            p[i].b = y1 - p[i].k*p[i].lx;
        }
        else{
            double _cos=(c[i].r-c[i+1].r)/(c[i+1].x-c[i].x);
            double d=_cos*c[i].r;
            p[i].lx=c[i].x+d;
            double y1 = sqrt(c[i].r*c[i].r - d*d);
            d=_cos*c[i+1].r;
            double y2 = sqrt(c[i+1].r*c[i+1].r - d*d);
            p[i].ly=c[i+1].x+d;
            p[i].k = (y1-y2)/(p[i].lx-p[i].ly);
            p[i].b = y1 - p[i].k*p[i].lx;
        }
    }
    //*****
    double ll=c[1].x-c[1].r,rr=c[n].x;
    for(int i=1;i<=n;i++){
        double l=c[i].x-c[i].r,r=c[i].x+c[i].r;
        ll=min(l,ll);rr=max(rr,r);
    }
    printf("%.2lf",2*asr(ll,rr,1e-6));
    return 0;
}