hdu 4697 Convex hull 对移动的凸包积分 利用叉积的分配率
题意:对移动的凸包积分, 给你点数n和时间T,0时刻的点为p,每个点的速度(用向量给出)为v,问你每单位时间的平均面积。
思路:对于这类题目,应该很容易想到把时间分段分别求凸包。
1. 所以我们先求出所有时间点,即点集中有三点共线的时候,就是一个时间点。
对于i,j,k三个点共线,那么用叉积可以列出一个关于时间t的一元二次方程(可能a==0),
解出所有时间点即可。这里要注意技巧,很容易证明叉积是满足分配率的,
所以你在结构体point中定义一个叉乘运算,对于求方程的系数是非常方便而且不容易出错。
当然对于时间0和T也算一个时间点
2. 枚举所有相邻时间段,取时间段中的中间时间,对这一时刻求一次凸包,
算出当前时间在凸包上的点。然后我们用这些点求面积,面积是一个关于时间t的函数
为a*x^2+b*x+c,求面积的方法:以某个凸包上的点为中心,枚举所有相邻的两点做叉积,
跟求时间点的时候的方程类似,累加系数a,b,c。积分得a/3*x^3+b/2*x^2+c*x,
把两端时间代入相减累加到答案中。注意这里的面积是有向面积,这点坑了我好久。
灵活运用叉积的分配率可以很快地解决此题。
#include <cstdio> #include <cstring> #include <cmath> #include <algorithm> using namespace std; const double eps = 1e-9; inline int dcmp(double x) { if (fabs(x) < eps) return 0; return x > eps ? 1 : -1; } struct point { double x, y; int id; inline void in(int i) { scanf("%lf%lf", &x, &y); id = i; } point(double x = 0, double y = 0) : x(x), y(y) { } inline point operator -(const point &t) const { return point(x - t.x, y - t.y); } inline point operator +(const point &t) const { return point(x + t.x, y + t.y); } point operator /(const double &t) const { return point(x / t, y / t); } point operator *(const double &t) const { return point(x * t, y * t); } inline double operator *(const point &t) const { //叉乘 return x*t.y-y*t.x; } bool operator <(const point &t) const { if(y == t.y) return x < t.x; return y < t.y; } }p[55], v[55], tp[55], st[55]; double t[55*55*55], T, cur, ans; int sz, n, m; double a, b, c; inline double cross(const point &o, const point &a, const point &b) { return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x); } int graham(point *p, int n, point *st) { //凸包 sort(p, p + n); int top = 0; int i; for (i = 0; i < n; i++) { while (top >= 2 && cross(st[top - 2], st[top - 1], p[i]) < eps) top--; st[top++] = p[i]; } int t = top + 1; for (i = n - 2; i >= 0; i--) { while (top >= t && cross(st[top - 2], st[top - 1], p[i]) < eps) top--; st[top++] = p[i]; } return top; } inline void gao(int i, int j, int k) { //计算系数的函数 point a1 = p[i]-p[j], a2 = p[i]-p[k]; point b1 = v[i]-v[j], b2 = v[i]-v[k]; a += b1*b2; b += a1*b2+b1*a2; c += a1*a2; } inline void solve(double a, double b, double c) { //解方程,注意a==0的情况 double x; if(a == 0) { if(b!= 0) { x = -c/b; if(x >= 0 && x <= T) t[sz++] = x; } return; } double dlt = b*b-4*a*c; if(dlt < 0) return; if(dlt == 0) { x = -b*0.5/a; if(x >= 0 && x <= T) t[sz++] = x; return; } dlt = sqrt(dlt); x = 0.5*(-b+dlt)/a; if(x >= 0 && x <= T) t[sz++] = x; x = 0.5*(-b-dlt)/a; if(x >= 0 && x <= T) t[sz++] = x; } inline double F(double x) { //积分求值函数 return a*x*x*x/3.0+b*x*x/2.0+c*x; } int main() { int i, j, k; while( ~scanf("%d%lf", &n, &T)) { for(i = 0; i < n; i++) p[i].in(i), v[i].in(i); if(n <= 2) { printf("%.10f\n", 0.0); continue; } t[0] = 0; t[1] = T; //处理出所有时间点 sz = 2; for(i = 0; i < n; i++) for(j = i+1; j < n; j++) for(k = j+1; k < n; k++) { a = b = c = 0; gao(i, j, k); solve(a, b, c); } sort(t, t+sz); ans = 0; for(i = 0; i < sz-1; i++) { //枚举相邻时间段 cur = 0.5*(t[i]+t[i+1]); a = b = c = 0; for(j = 0; j < n; j++) { tp[j] = p[j]+v[j]*cur; tp[j].id = p[j].id; } m = graham(tp, n, st);//凸包求出该时间段在凸包上的点 for(j = 2; j < m; j++) gao(st[0].id, st[j-1].id, st[j].id); //求系数 ans += F(t[i+1])- F(t[i]);//手动积分算面积 } printf("%.10f\n", fabs(ans*0.5/T)); } return 0; }
补充:软件开发 , C++ ,