BZOJ 1502 NOI 2005 月下柠檬树

题目大意

给定一些圆和梯形，求它们面积的并

解答

参考用辛普森公式积分

$$\int_a^bf(x)dx \approx {b-a\over6}[f(a)+f({a+b\over2})+f(b)]$$
利用二分思想，当：
$$\vert\int_a^{mid}f(x)dx + \int_{mid}^bf(x)dx - \int_a^bf(x)dx\vert \le eps$$
我们就认为这个数就是积分的答案

参考代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <algorithm>

using namespace std;

const double eps = 1e-8;

class Line
{
    public:
        double x1, y1, x2, y2;
};

int n;
double alpha;
double h[800];
double r[800];
double ll = 0, rr = 0;
Line ls[1600];
int lt = 0;

inline double sqr(double x)
{
    return x*x;
}

void readin()
{
    scanf("%d %lf", &n, &alpha);
    double tmp = 1/tan(alpha);
    double x;
    for (int i = 0; i <= n; i++) {
        scanf("%lf", &x);
        x *= tmp;
        h[i] = x + (i ? h[i-1] : 0);
    }
    for (int i = 0; i < n; i++) {
        scanf("%lf", &r[i]);
        ll = min(ll, h[i]-r[i]);
        rr = max(rr, h[i]+r[i]);
    }
    r[n] = 0;
    rr = max(rr, h[n]);
}

void preWork()
{
    double d;
    for (int i = 0; i < n; i++) {
        d = h[i+1]-h[i];
        if (d > fabs(r[i+1]-r[i])) {
            double sinbeta = (r[i]-r[i+1])/d;
            double cosbeta = sqrt(1-sqr(sinbeta));
            ls[lt].x1 = h[i]+r[i]*sinbeta;
            ls[lt].y1 = r[i]*cosbeta;
            ls[lt].x2 = h[i+1]+r[i+1]*sinbeta;
            ls[lt++].y2 = r[i+1]*cosbeta;
        }
    }
}

double f(double l)
{
    double ret = 0;
    for (int i = 0; i < n; i++) {
        if (fabs(h[i]-l) < r[i])
            ret = max(ret, sqrt(sqr(r[i])-sqr(h[i]-l)));
    }
    for (int i = 0; i < lt; i++) {
        if (ls[i].x1 < l && ls[i].x2 > l) {
            double k = (ls[i].y2-ls[i].y1) / (ls[i].x2-ls[i].x1);
            ret = max(ret, ls[i].y1+k*(l-ls[i].x1));
        }
    }
    return ret;
}

double simpson(double lp, double rp)
{
    double mid = (lp+rp)/2;
    return (rp-lp)*(f(lp)+4*f(mid)+f(rp))/6;
}

double work(double lp, double rp)
{
    double mid = (lp+rp)/2;
    double cmp = simpson(lp, mid) + simpson(mid, rp);
    double cmp1 = simpson(lp, rp);
    if (fabs(cmp1-cmp) < eps)
        return cmp;
    return work(lp, mid) + work(mid, rp);
}

int main()
{
    readin();
    preWork();
    printf("%.2lf", work(ll, rr)*2);
    return 0;
}