BZOJ--1502（Simpson积分）

2015-07-20 20:01:26

【传送门】

题意：给出多个圆台构成的柠檬树，求其在底面上的投影的面积。给出圆台的参数和月光与底面的夹角。（月光考虑为平行光）

思路：考虑圆台 / 圆锥在底面上的平行投影是若干个圆的并，他的轮廓就是若干个圆以及相邻圆之间的两条外公切线构成，所以列出所有圆以及公切线，考虑对X轴微分，f(x)为对应的y值，然后Simpson积分，由于只算了上半部分，最后答案乘2即可。

#include <cstdio>
#include <ctime>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <string>
#include <iostream>
#include <algorithm>
using namespace std;

#define getmid(l,r) ((l) + ((r) - (l)) / 2)
#define MP(a,b) make_pair(a,b)
#define PB push_back

typedef long long ll;
typedef pair<int,int> pii;
const double eps = 1e-5;
const int INF = (1 << 30) - 1;
const int MAXN = 510;
const double PI = acos(-1.0);

int n,L_cnt,C_cnt;
double alpha,r[MAXN],h[MAXN];
map<double,double> mp;

struct Line{
    double l,r,k,b;
}L[MAXN];

struct Circle{
    double a,rr;
    double left,right;
}C[MAXN];

inline double F(double x){
    if(mp.find(x) != mp.end()) return mp[x];
    double tmax = 0.0;
    for(int i = 1; i <= L_cnt; ++i){
        if(x < L[i].l || x > L[i].r) continue;
        tmax = max(tmax,x * L[i].k + L[i].b);
    }
    for(int i = 1; i <= C_cnt; ++i){
        if(x < C[i].left || x > C[i].right) continue;
        double y = sqrt(C[i].rr - (x - C[i].a) * (x - C[i].a));
        tmax = max(tmax,y);
    }
    mp[x] = tmax;
    return tmax;
}

inline double Simpson(double tl,double tr){
    return (F(tl) + 4.0 * F(getmid(tl,tr)) + F(tr)) * (tr - tl) / 6.0;
}

double Integral(double tl,double tr){
    double sum = Simpson(tl,tr);
    double mid = getmid(tl,tr);
    double suml = Simpson(tl,mid);
    double sumr = Simpson(mid,tr);
    if(fabs(sum - suml - sumr) < eps) return sum;
    else return Integral(tl,mid) + Integral(mid,tr);
}

int main(){
    scanf("%d%lf",&n,&alpha);
    double tmp = 1.0 / tan(alpha);
    for(int i = 0; i <= n; ++i){
        scanf("%lf",&h[i]);
        h[i] *= tmp;
    }
    for(int i = 1; i <= n; ++i){
        scanf("%lf",&r[i]);
    }
    r[n + 1] = 0.0;
    double sum_h = h[0];
    double tmin = 0,tmax = 0;
    for(int i = 1; i <= n; ++i){
        //Circle
        ++C_cnt;
        C[C_cnt].left = sum_h - r[i];
        tmin = min(tmin,C[C_cnt].left);
        C[C_cnt].right = sum_h + r[i];
        tmax = max(tmax,C[C_cnt].right);
        C[C_cnt].a = sum_h;
        C[C_cnt].rr = r[i] * r[i];
        //Line
        if(h[i] > 0){
            ++L_cnt;
            double sina = (r[i + 1] - r[i]) / h[i];
            double cosa = sqrt(h[i]*h[i] - (r[i+1]-r[i])*(r[i+1]-r[i])) / h[i];
            double Y1 = r[i] * cosa;
            L[L_cnt].l = sum_h - sina * r[i];
            L[L_cnt].r = sum_h + h[i] - sina * r[i + 1];
            L[L_cnt].k = sina / cosa;
            L[L_cnt].b = Y1 - L[L_cnt].k * L[L_cnt].l;
        }
        sum_h += h[i];
    }
    tmax = max(tmax,sum_h);
    printf("%.2f\n",2.0 * Integral(tmin,tmax));
    return 0;
}