50个点,每个点有一个恒定速度,给定一个时刻T
求,其中 S(t) 为所有点在t时刻组成的凸包面积。
找到所有特殊的时间点: 只有在存在三点共线的时候凸包上的点才会改变。
如果凸包上的点保持不变,其面积随时间变化是一个一元二次函数(考虑相邻两点的位置向量是一个关于时间的一次函数,叉乘后得到一个一元二次的式子)。
所以求出a,b,c后求积分即可。求a,b,c时可以取开始时间点,结束时间点,再取一个中间点。这个过程在网上看的代码,比较巧妙。。
精度问题,除以2 用乘以0.5代替
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <ctime>
#include <algorithm>
#include <iostream>
#include <vector>
#define pb push_back
#define sqr(x) ((x)*(x))
#define cube(x) ((x)*(x)*(x))
using namespace std;
const double eps = 1e-9;
const int maxN = 50;
int n, T;
vector<double> t;
int dcmp(double x){
if (fabs(x)<eps) return 0;
return x>0?1:-1;
}
struct Point{
double x, y;
int id;
Point(){}
Point(double xx, double yy):x(xx), y(yy){}
void input(){ scanf("%lf%lf", &x, &y); }
friend Point operator + (const Point& a, const Point& b){
return Point(a.x+b.x, a.y+b.y);
}
friend Point operator - (const Point& a, const Point& b){
return Point(a.x-b.x, a.y-b.y);
}
friend double operator * (const Point& a, const Point& b){
return a.x*b.y - a.y+b.x;
}
friend Point operator * (const Point& a, const double& k){
return Point(a.x*k, a.y*k);
}
friend bool operator < (const Point& a, const Point& b){
return a.x<b.x || a.x==b.x&&a.y<b.y;
}
}p[maxN+2], q[maxN+2], v[maxN+2], ch[maxN+2];
int det(const Point& a, const Point& b){
return dcmp(a.x*b.y - a.y*b.x);
}
int Hull(int n, Point p[], Point ch[]){
sort(p+1, p+1+n);
int m=0;
for (int i=1; i<=n; i++){
while (m>1 && det(ch[m]-ch[m-1], p[i]-ch[m-1])<=0) m--;
ch[++m] = p[i];
}
int pre=m;
for (int i=n; i; i--){
while (m>pre && det(ch[m]-ch[m-1], p[i]-ch[m-1])<=0) m--;
ch[++m] = p[i];
}
return m;
}
void get_abc(double& a, double& b, double& c, int i, int j, int k){
Point p1=p[j]-p[i], v1=v[j]-v[i];
Point p2=p[k]-p[i], v2=v[k]-v[i];
a = v1*v2;
b = p1*v2 + p2*v1;
c = p1*p2;
}
void solve(double a, double b, double c){
double x;
if (a==0){
if (dcmp(b)){
x = -c/b;
if (x>=0 && x<=T) t.pb(x);
}
return;
}
double delta = b*b - 4*a*c;
if (delta<0) return;
if (!dcmp(delta)){
x = -b*0.5/a;
if (x>=0 && x<=T) t.pb(x);
return;
}
double tmp = sqrt(delta);
x = (-b+tmp)*0.5/a;
if (x>=0 && x<=T) t.pb(x);
x = (-b-tmp)*0.5/a;
if (x>=0 && x<=T) t.pb(x);
}
int main(){
while (scanf("%d%d", &n, &T) != EOF){
for (int i=1; i<=n; i++){
p[i].input();
v[i].input();
}
t.clear();
t.pb(0); t.pb(T);
for (int i=1; i<=n; i++)
for (int j=i+1; j<=n; j++)
for (int k=j+1; k<=n; k++){
double a, b, c;
get_abc(a, b, c, i, j, k);
solve(a, b, c);
}
int tot=t.size();
sort(t.begin(), t.end());
double ans=0;
for (int i=0; i<tot-1; i++){
double mt = (t[i]+t[i+1])*0.5;
for (int i=1; i<=n; i++){
q[i] = p[i] + v[i]*mt;
q[i].id = i;
}
int m = Hull(n, q, ch);
for (int i=1; i<=m; i++)
printf("\t%d\n", ch[i].id);
for (int i=1; i<m; i++){
double a, b, c;
get_abc(a, b, c, ch[1].id, ch[i].id, ch[i+1].id);
ans += a/3*(cube(t[i+1]) - cube(t[i])) + b*0.5*(sqr(t[i+1]) - sqr(t[i])) + c*(t[i+1]-t[i]);
}
}
printf("%.10f\n", ans/T*0.5);
}
return 0;
}