HDU4697 Convex hull

来源:互联网 发布:java语言绝对值 编辑:程序博客网 时间:2024/06/05 02:05

X姐说没那么难写。

但是A完这题,我觉得我可以去死了= =

该死的精度问题。


基本解题思路是这样的,对任意三点共线的情况求出其发生的时间t,将其标记为一个事件,可以证明,只有发生事件的时候才有可能使凸包的形状发生变化。

对所有事件按照时间排序,求出所有相邻事件间的凸包形状(在凸包上的点集),并根据该点集和凸包面积的积分求出答案。


本题的难点在于:

1. 精度问题(1e-9丧心病狂!!!double和int的转换,写template有优势)

2. 解方程问题(学习了套公式求参数)

3. 求积分问题(也是套公式求参数,这个求参数的函数真有用)


#pragma comment(linker,"/STACK:102400000,102400000")#include <cstdio>#include <iostream>#include <cstring>#include <cmath>#include <ctime>#include <vector>#include <set>#include <queue>#include <stack>#include <map>#include <algorithm>using namespace std;typedef long long ll;const int maxn = 60;const double eps = 1e-9;template <class Type>struct Point{Type x, y;int id;Point(Type x=0, Type y=0): x(x), y(y) {}};template <class Type>bool operator < (const Point<Type> &a, const Point<Type> &b) {return a.x < b.x || (a.x == b.x && a.y <b.y);}int dcmp(double x) {if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;}template <class Type>Point<Type> operator + (Point<Type> A, Point<Type> B) {return Point<Type>(A.x+B.x, A.y+B.y);}template <class Type>Point<Type> operator - (Point<Type> A, Point<Type> B) {return Point<Type>(A.x-B.x, A.y-B.y);}template <class Type>Point<Type> operator * (Point<Type> A, Type p) {return Point<Type>(A.x*p, A.y*p);}template <class Type>Type Dot(Point<Type> A, Point<Type> B) {return A.x*B.x + A.y*B.y;}template <class Type>Type Cross(Point<Type> A, Point<Type> B) {return A.x*B.y - A.y*B.x;}template <class Type>int ConvexHull(Point<Type> *A, int n, Point<Type> *ch){sort(A, A+n);int m = 0;for (int i = 0; i < n; i++){while (m > 1 && dcmp(Cross(ch[m-1]-ch[m-2], A[i]-ch[m-1])) <= 0) m--;ch[m++] = A[i];}int k = m;for (int i = n-2; i >= 0; i--){while (m > k && dcmp(Cross(ch[m-1]-ch[m-2], A[i]-ch[m-1])) <= 0) m--;ch[m++] = A[i];}return m;}int n, T, m;Point<int> p[maxn], v[maxn];Point<double> q[maxn], ch[maxn];double t[250000];int tot;double ans;template <class Type>void get_coef(Type &a, Type &b, Type &c, Type d, Type k_1, Type b_1, Type k_2, Type b_2) // d * (k_1 x + b_1) * (k_2 x + b_2){    a += d * (k_1 * k_2);    b += d * (k_1 * b_2 + k_2 * b_1);    c += d * b_1 * b_2;}inline void solve(int a, int b, int c){    if (a == 0)    {        if (b != 0) t[tot++] = -1.0 * c / b;        return;    }    int delta = b*b-4*a*c;    if (delta < 0) return;    if (delta == 0)    {        t[tot++] = 0.5 * (-b) / a;        return;    }    double tmp = sqrt(1.0*delta);    t[tot++] = 0.5 * (-b - tmp) / a;    t[tot++] = 0.5 * (-b + tmp) / a;    return;}inline double integral(int i, int j, double t1, double t2){    int a = 0, b = 0, c = 0;    get_coef(a, b, c, 1, v[i].x, p[i].x, v[j].y, p[j].y);    get_coef(a, b, c, -1, v[j].x, p[j].x, v[i].y, p[i].y);    return a * (t2*t2*t2-t1*t1*t1)/3.0 + b * (t2+t1)*(t2-t1)/2.0 + c * (t2-t1);}int main(){    while (scanf("%d%d", &n, &T) == 2)    {        for (int i=0;i<n;i++)        {            scanf("%d%d%d%d", &p[i].x, &p[i].y, &v[i].x, &v[i].y);            p[i].id = i;        }        if (n <= 2) {printf("%.10f\n", 0); continue;}        t[0] = 0; t[1] = T;        tot = 2;        for (int i=0;i<n;i++)            for (int j=i+1;j<n;j++)                for (int k=j+1;k<n;k++)                {                    int a = 0, b = 0, c = 0;                    get_coef(a, b, c, 1, v[i].x-v[j].x, p[i].x-p[j].x, v[i].y-v[k].y, p[i].y-p[k].y);                    get_coef(a, b, c, -1, v[i].x-v[k].x, p[i].x-p[k].x, v[i].y-v[j].y, p[i].y-p[j].y);                    solve(a, b, c);                }        sort(t, t+tot);        ans = 0;        for (int k=0;k<tot-1;k++)        {            double mt = 0.5 * (t[k] + t[k+1]);            if (mt < 0 || mt > T) continue;            for (int i=0;i<n;i++)            {                q[i] = Point<double>(p[i].x + v[i].x * mt, p[i].y + v[i].y * mt);                q[i].id = p[i].id;            }            m = ConvexHull(q, n, ch);            for (int i=1;i<m;i++) ans += integral(ch[i-1].id, ch[i].id, t[k], t[k+1]);        }        printf("%.10f\n", abs(ans * 0.5 / T));    }    return 0;}