題意：對移動的凸包積分， 給你點數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;
}
