BZOJ 1502: [NOI2005]月下柠檬树 [辛普森积分 解析几何 圆]

1502: [NOI2005]月下柠檬树

Time Limit: 5 Sec  Memory Limit: 64 MB
Submit: 1070  Solved: 596
[Submit][Status][Discuss]

Description

Input

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

Output

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

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

HINT

1≤n≤500，0.3

%%%Claude画图真好看 http://blog.csdn.net/wzq_qwq/article/details/48310417

把每条线段和每个点的投影找出来，然后计算F函数时遍历所有线段和圆找最大值行了

这里的线保存了k和b，解析几何的感觉

找点和圆的切线用了射影定理，小新讲过然而并不会(因为我只会三角形相似)，又去学了一下

概述图中，在Rt△ABC中，∠ABC=90°，BD是斜边AC上的高，则有射影定理如下：

BD²=AD·CD

AB²=AC·AD

BC²=CD·AC

找圆的公切线直接用三角函数....

//
//  main.cpp
//  bzoj1502
//
//  Created by Candy on 2017/2/1.
//  Copyright © 2017年 Candy. All rights reserved.
//

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
typedef long long ll;
const int N=2005;
const double INF=1e9;
const double eps=1e-8;
const double pi=acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1; c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘; c=getchar();}
    return x*f;
}
inline int sgn(double x){
    if(abs(x)<eps) return 0;
    else return x<0?-1:1;
}
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
    void print(){printf("%lf %lf\n",x,y);}
};
typedef Vector Point;

struct Line{
    Point s,t;
    double k,b;
    Line(){}
    Line(Point a,Point c):s(a),t(c){
        k=(t.y-s.y)/(t.x-s.x);
        b=s.y-k*s.x;
    }
    double f(double x){return k*x+b;}
}L[N];
int nl;
struct Circle{
    double x,r;
    Circle(){}
    Circle(double x,double r):x(x),r(r){}
}C[N];
void addCommonTangent(Circle a,Circle b){
    nl++;
    double sina=(a.r-b.r)/(b.x-a.x);
    double cosa=sqrt(1-sina*sina);
    double tana=sina/cosa;
    L[nl].s=Point(a.x+a.r*sina,a.r*cosa);
    L[nl].t=Point(b.x+b.r*sina,b.r*cosa);
    L[nl].k=-tana;
    L[nl].b=L[nl].s.y-L[nl].k*L[nl].s.x;
}
int n;
double alpha,h[N],lb=INF,rb;
inline double F(double x){
    double re=0;
    for(int i=1;i<=nl;i++) if(x>=L[i].s.x&&x<=L[i].t.x) re=max(re,L[i].f(x));
    for(int i=1;i<=n;i++) if(x>=C[i].x-C[i].r&&x<=C[i].x+C[i].r)
        re=max(re,sqrt(C[i].r*C[i].r-(x-C[i].x)*(x-C[i].x)));
    return re;
}
inline double cal(double l,double r){
    return (F(l)+F(r)+4*F((l+r)/2))*(r-l)/6;
}
double simpson(double l,double r,double now){
    double mid=(l+r)/2,p=cal(l,mid),q=cal(mid,r);
    if(abs(now-p-q)<eps) return now;
    else return simpson(l,mid,p)+simpson(mid,r,q);
}

Point p;
int main(int argc, const char * argv[]){
    scanf("%d%lf",&n,&alpha);
    for(int i=1;i<=n+1;i++) scanf("%lf",&h[i]),h[i]+=h[i-1];
    for(int i=1;i<=n;i++) scanf("%lf",&C[i].r);
    double ta=tan(alpha);
    p=Point(h[n+1]/ta,0);
    rb=max(rb,p.x);
    {
        C[n].x=h[n]/ta;
        double x=C[n].x,r=C[n].r;
        lb=min(lb,x-r);
        rb=max(rb,x+r);
        if(x+r<p.x){
            double l=r*r/(p.x-x);//she ying ding li
            double h=sqrt(r*r-l*l);
            L[++nl]=Line(Point(x+l,h),p);
        }
    }
    for(int i=n-1;i>=1;i--){
        C[i].x=h[i]/ta;
        double x=C[i].x,r=C[i].r;
        lb=min(lb,x-r);
        rb=max(rb,x+r);
        if(sgn(C[i+1].x-C[i].x-abs(C[i+1].r-C[i].r))>0)//d-abs(R-r)<=0 nei han
            addCommonTangent(C[i],C[i+1]);
    }
    printf("%.2f\n",2*simpson(lb,rb,cal(lb,rb)));
    return 0;
}