【计算几何】【simpson自适应公式】【NOI2004】降雨量

【问题描述】
M国是个多雨的国家，尤其是P城，频繁的降雨给人们的出行带来了不少麻烦。为了方便行人雨天过马路，有关部门在每处人行横道的上空都安装了一种名为“自动伞”的装置。（如图1所示）

每把自动伞都可以近似地看作一块长方形的板，厚度不计。这种伞有相当出色的吸水能力，落到上面的雨水都会完全被伞顶的小孔吸入，并通过管道运走。不下雨时，这些伞闲置着。一旦下雨，它们便立刻开始作匀速率直线往返运动：从马路的一边以固定的速率移动到另一边，再从另一边以相同的速率返回，如此往复，直到雨停为止。任何时刻自动伞都不会越过马路的边界。有了自动伞，下雨天没带伞的行人只要躲在伞下行走，就不会被雨淋着了。
由于自动伞的大小有限，当需要使用自动伞过马路的人比较多时，一把自动伞显然是不够的，所以有关部门在几处主要的人行横道上空安装了多把自动伞。每把自动伞的宽度都等于人行横道的宽度，高度各不相同，长度不一定相同，移动速率也不一定相同。
现在已知一处人行横道的详细情况，请你计算从开始下雨到T秒钟后的这段时间内，一共有多少体积的雨水降落到人行横道上。

【输入文件】
第一行有四个整数N，W，T，V。N表示自动伞的数目，W表示马路的宽度，T表示需要统计从开始下雨到多长时间后的降雨情况，V表示单位面积单位时间内的降雨体积。
为了描述方便，我们画出了一个如图2所示的天空中五把伞的剖面图，取马路左边界为原点，取向右为x轴正方向，取向上为y轴正方向，建立平面直角坐标系。这样，每把自动伞都可以看作平面上的一条线段。

接下来的N行，每行用三个整数描述一把自动伞。第一个数x是伞的初始位置，用它左端点的横坐标表示。第二个数l是伞的长度，即x方向上的尺寸。第三个数v是伞的速度，v的大小表示移动的速率。如果v>0，表示开始时伞向右移动；如果v<0，表示开始时伞向左移动；如果v=0，表示伞不动。

【输出文件】
输出文件只包含一个实数，表示从开始下雨到T秒钟后，一共有多少体积的水降落到人行横道上。输出结果精确到小数点后第二位。

【约定】
雨点均匀地匀速竖直下落
自动伞和马路两者都是水平的
自动伞的宽度和人行横道的宽度相等，都等于1

N <= 10
W <= 100
T <= 100
V <= 50

 所有自动伞的往返次数之和不超过250，一来一回算一个往返。

【样例输入】
2 4 3 10
0 1 1
3 1 -1

【样例输出】
65.00

以位置为横坐标，时间为纵坐标，可以得到一些相交的区域，如图所示。


那么可以很容易用simpson自适应公式算出有色区域的面积（若干个平行四边形的面积并），而题目所求即为总面积减去有色区域的面积（无色区域的面积）乘以V即可。

考试的时候把斜率弄反了，结果全错……
但是这样计算精度不够，只能过90分。
代码：

/*****************************\
 * @prob: NOI2004 rainfall   *
 * @auth: Wang Junji         *
 * @stat: WA: 90             *
 * @date: June. 10th, 2012   *
 * @memo: simpson自适应公式   *
\*****************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>

const int maxN = 20, maxW = 110, maxM = 1010;
const double zero = 1e-10, INF = 1e198;

struct UM
{
    int L0, R0, v, len;
    double T0, T;
} um[maxN];
struct Seg
{
    double L, R; Seg() {}
    Seg(double L, double R): L(L), R(R) {}
};
int n, m, W, T, V;

inline bool work(double& x, double& L, double& R,
                 double L0, double R0,
                 double T0, double T1, double v)
{
    L = (x - L0) / v, R = (x - R0) / v; //
    if (L - R > zero) std::swap(L, R);
    if (R - T0 < -zero || L - T1 > zero) return 0;
    if (R - T1 > zero) R = T1;
    if (L - T0 < -zero) L = T0;
    if (L - R > -zero) return 0;
    return 1;
}

inline bool cmp(const Seg& a, const Seg& b) {return a.L < b.L;}

double f(double x)
{
    int cnt = 0; static Seg seg[maxM];
    for (int i = 0; i < n; ++i)
    {
        if (!um[i].v || um[i].len == W) //
        {
            if (x >= um[i].L0 && x <= um[i].R0) return T;
            else continue;
        }
        double L, R;
        if (um[i].T0 - T > -zero)
        {
            if (work(x, L, R, um[i].L0, um[i].R0, 0., T, um[i].v))
                seg[cnt++] = Seg(L, R);
        }
        else
        {
            if (work(x, L, R, um[i].L0, um[i].R0, 0., um[i].T0, um[i].v))
                seg[cnt++] = Seg(L, R);
            double tmp = um[i].T0, v = -um[i].v, len = um[i].len,
                   M = v < 0 ? W * 2 - len : len,
                   l0 = M - um[i].R0, r0 = M - um[i].L0;
            while (tmp - T < -zero)
            {
                if (work(x, L, R, l0, r0, tmp, std::min(tmp + um[i].T, (double)T), v))
                    seg[cnt++] = Seg(L, R);
                v = -v, M = v < 0 ? W * 2 - len : len;
                double l = M - r0, r = M - l0;
                l0 = l, r0 = r, tmp += um[i].T;
            }
        }
    }
    if (!cnt) return 0;
    std::sort(seg, seg + cnt, cmp);
    double L = seg[0].L, R = seg[0].R, len = 0;
    for (int i = 1; i < cnt; ++i)
    if (seg[i].L - R > zero) len += R - L, L = seg[i].L, R = seg[i].R;
    else if (seg[i].R - R > zero) R = seg[i].R;
    return len += R - L;
}

inline double simpson(double L, double R, double fL, double fM, double fR)
{return (R - L) * (fL + 4 * fM + fR) / 6;}

double calc(double L, double Mid, double R, double fL, double fM, double fR, double pre)
{
    double LM = (L + Mid) / 2, RM = (Mid + R) / 2,
           fLM = f(LM), sL = simpson(L, Mid, fL, fLM, fM),
           fRM = f(RM), sR = simpson(Mid, R, fM, fRM, fR);
    if (fabs(sL + sR - pre) < zero) return pre;
    else return calc(L, LM, Mid, fL, fLM, fM, sL) +
                calc(Mid, RM, R, fM, fRM, fR, sR);
}

inline void get_T0(double& T0, int L, int R, int v)
{
    if (!v) T0 = INF;
    else if (v < 0) T0 = -(double)L / v;
    else if (v > 0) T0 = (double)(W - R) / v;
    return;
}

int main()
{
    freopen("rainfall.in", "r", stdin);
    freopen("rainfall.out", "w", stdout);
    scanf("%d%d%d%d", &n, &W, &T, &V);
    for (int i = 0; i < n; ++i)
        scanf("%d%d%d", &um[i].L0, &um[i].len, &um[i].v),
        um[i].R0 = um[i].L0 + um[i].len,
        um[i].T = fabs((double)(W - um[i].len) / um[i].v),
        get_T0(um[i].T0, um[i].L0, um[i].R0, um[i].v);
    double L = 0, R = W, Mid = (L + R) / 2,
           fL = f(L), fR = f(R), fM = f(Mid),
           pre = simpson(L, R, fL, fM, fR);
    printf("%.2lf\n", (W * T - calc(L, Mid, R, fL, fM, fR, pre)) * V);
    return 0;
}

以时间为横坐标，位置为纵坐标，建立平面直角坐标系。
样例如图所示。

目标是求出无色区域的面积，需要先求出有色区域的面积。
那么可以把所有平行四边形间的交点全部算出来，从左到右排个序，这样就把整个x轴分成了若干段，每一段的面积都可以等效成一个梯形来计算。

当然计算梯形的底边是通过计算该位置的竖线被所有平行四边形的截距而得出的。
代码：

/*****************************\
 * @prob: NOI2004 rainfall   *
 * @auth: Wang Junji         *
 * @stat: Accepted.          *
 * @date: June. 11th, 2012   *
 * @memo: 计算几何            *
\*****************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>

const int maxN = 20, maxW = 110, maxM = 2010;
const double zero = 1e-10, INF = 1e198;

inline int sgn(const double b) {return fabs(b) < zero ? 0 : b > 0 ? 1 : -1;}

struct vec
{
    double x, y; vec() {}
    vec(double x, double y): x(x), y(y) {}
    vec operator+(const double& b) const {return vec(x, y + b);}
    vec operator-(const vec& b) const {return vec(x - b.x, y - b.y);}
    double operator*(const vec& b) const {return x * b.y - y * b.x;}
    bool operator<(const vec& b) const
    {return sgn(x - b.x) < 0 || (sgn(x - b.x) == 0 && sgn(y - b.y) < 0);}
} p[maxN][maxM];
struct UM
{
    int L0, v, len;
    double T0, T;
} um[maxN];
struct Seg
{
    double L, R; Seg() {}
    Seg(double L, double R): L(L), R(R) {}
};
double inter[maxM], height[maxM];
int cnt[maxN], len[maxN], cnt_inter, n, m, Width, Time, V;

inline bool cmp(const Seg& a, const Seg& b) {return a.L < b.L;}

inline double get_T0(int L, int R, int v)
{
    if (!v) return INF;
    else if (v < 0) return -(double)L / v;
    else if (v > 0) return (double)(Width - R) / v;
}

inline bool test_cross(vec& A, vec& B, vec& C, vec& D)
{
    vec AC = C - A, AD = D - A, BC = C - B, BD = D - B;
    return (sgn(AC * AD) ^ sgn(BC * BD)) == -2
        && (sgn(AC * BC) ^ sgn(AD * BD)) == -2;
}

inline void Find_inter(vec A, vec B, vec C, vec D)
{
    if (test_cross(A, B, C, D))
        inter[cnt_inter++] =
            (B.x - A.x) * fabs(((C - A) * (D - A))
            / ((B - A) * (D - C))) + A.x;
    return;
}

inline void calc_height(int ths)
{
    double x = inter[ths];
    static Seg seg[maxM]; int cnt_seg = 0;
    for (int i = 0; i < n; ++i)
    {
        int pos = std::lower_bound(p[i], p[i] + cnt[i], vec(x, INF)) - p[i] - 1;
        double tmp = (p[i][pos + 1].y - p[i][pos].y)
            / (p[i][pos + 1].x - p[i][pos].x)
            * (x - p[i][pos].x) + p[i][pos].y;
        seg[cnt_seg++] = Seg(tmp, tmp + len[i]);
    }
    std::sort(seg, seg + cnt_seg, cmp);
    double L = seg[0].L, R = seg[0].R, len = 0;
    for (int i = 1; i < cnt_seg; ++i)
    if (seg[i].L - R > zero) len += R - L, L = seg[i].L, R = seg[i].R;
    else if (seg[i].R - R > zero) R = seg[i].R;
    height[ths] = len += R - L;
    return;
}

inline bool eq(const double& a, const double& b) {return !sgn(a - b);}

int main()
{
    freopen("rainfall.in", "r", stdin);
    freopen("rainfall.out", "w", stdout);
    scanf("%d%d%d%d", &n, &Width, &Time, &V);
    for (int i = 0, y, v; i < n; ++i)
    {
        scanf("%d%d%d", &y, len + i, &v);
        double T = fabs((double)(Width - len[i]) / v);
        p[i][(cnt[i] = 0)++] = vec(0, y);
        if (len[i] == Width || v == 0) p[i][1] = vec(Time, y);
        else
        {
            double X = get_T0(y, y + len[i], v),
                   Y = v > 0 ? Width - len[i] : 0;
            if (X - Time > -zero)
            {
                p[i][cnt[i]] = vec(Time, y + v * Time);
                continue;
            }
            p[i][cnt[i]++] = vec(X, Y), inter[cnt_inter++] = X;
            //一定要把各个平行四边形的拼接处也算上。
            while (Y = Width - len[i] - Y, v = -v, (X += T) - Time < -zero)
                inter[cnt_inter++] = X, p[i][cnt[i]++] = vec(X, Y);
            p[i][cnt[i]] = vec(Time, Y - v * (X - Time));
        }
        for (int j = 0; j < cnt[i]; ++j)
        for (int i1 = 0; i1 < i; ++i1)
        for (int j1 = 0; j1 < cnt[i1]; ++j1)
            Find_inter(p[i][j], p[i][j + 1], p[i1][j1], p[i1][j1 + 1]),
            Find_inter(p[i][j] + len[i], p[i][j + 1] + len[i],
                    p[i1][j1], p[i1][j1 + 1]),
            Find_inter(p[i][j], p[i][j + 1],
                    p[i1][j1] + len[i1], p[i1][j1 + 1] + len[i1]),
            Find_inter(p[i][j] + len[i], p[i][j + 1] + len[i],
                    p[i1][j1] + len[i1], p[i1][j1 + 1] + len[i1]);
    }
    inter[cnt_inter++] = 0;
    std::sort(inter, inter + cnt_inter);
    cnt_inter = std::unique(inter, inter + cnt_inter, eq) - inter;
    inter[cnt_inter] = Time;
    for (int i = 0; i < cnt_inter + 1; ++i) calc_height(i);
    double area = 0;
    for (int i = 0; i < cnt_inter; ++i)
        area += (inter[i + 1] - inter[i]) * (height[i] + height[i + 1]) / 2;
    printf("%.2lf\n", (Width * Time - area) * V);
    return 0;
}