bzoj1094 粒子运动 计算几何

版权声明:转载请注明:http://blog.csdn.net/lych_cys https://blog.csdn.net/lych_cys/article/details/50785713

       看了好多份博客。。现在来这里详细讲一下。(附:博客1博客2;)

       求任意两点任意时刻的最近距离,实际上我们可以转化为枚举两个点,求这两个点的历史最近距离。

       显然, 对于一个点,它的运动轨迹是k条折线,不妨先来看如何求出一个点的运动轨迹。见下图:


       我们作出向量AB到达边界以后反弹得到的向量BA'。

       首先考虑一直点A,和向量AB(即运动方向),如何得到点B的坐标,已知点A(px,py)和向量AB(vx,vy),假设圆心为(0,0)(否则将点A平移),设时间为λ,那么列出方程:

       (px+λvx)^2+(py+λcy)^2=r^2,解出λ,那么B(px+λvx,py+λvy)。

       然后观察向量AB和向量BA’,发现它们相当于半径OB的法向量BC(即BC⊥OB)轴对称,然后我们考虑两个向量u,v,关于向量t轴对称,已知u,t,如何求出v,见下图:

       

       我们将向量AB平移得到BA'',可以发现BA'显然和BA''关于向量BC轴对称。

       考虑两个向量的点积dot(p,q)=p.x*q.x+p.y*q.y,根据点积的定义为len(p)*len(q)*cos(p,q),我们可以得到dot(BA',BD)=dot(BD,BA''),而后者已知。

       然后作A'D⊥BD,那么dot(BA',BC)=len(BA')*len(BC)*cos∠A'BC=len(BA')*len(BD),不妨令len(BA')=len(BA'')(向量所以长度没有关系),那么可以得到len(BD),那么可以用len(BD)和len(BC)的长度关系得到向量BD,倍长得到BB'。

       将BA''平移得到A'B',那么显然BA'+A'B'=BB',因此BA'=BB'-A'B',得到向量BA',结束。


       然后我们就可以用上述方法得到k个碰撞点,以及碰撞之后的向量了,即得到了一个粒子运动的折线。

       考虑两个点x,y的最近距离,我们可以根据碰撞点把两条折线的距离转化为2*k对线段的距离。那么如何求一对线段的距离呢?

       设线段为p->p+v,p为一个端点,p+v为另一个端点,那么v为线段的向量。那么线段上面任一点可以用p+λv表示,其中λ∈[0,1]。那么x和y的距离可以转化为二次函数(px.x+λvx.x-py.x-λvy.x)^2+(px.y+λvx.y-py.y-λvt.y)^2,然后就可以用二次函数求最小值了(注意退化成一次函数的情况;注意对称轴取不到的情况)。

       虽然讲了那么多,代码还是只有2k的辣~\(≧▽≦)/~!!

AC代码如下(也不知道讲得清不清楚):

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#define N 205
using namespace std;

struct point{ double x,y; }o;
struct line{ point p,v; }a[N][N];
int n,m; double r,c[N][N];
point operator -(point u,point v){
	point t; t.x=u.x-v.x; t.y=u.y-v.y; return t;
}
point operator +(point u,point v){
	point t; t.x=u.x+v.x; t.y=u.y+v.y; return t;
}
point operator *(point u,double t){
	point v; v.x=u.x*t; v.y=u.y*t; return v;
}
double dot(point u,point v){ return u.x*v.x+u.y*v.y; }
double crs(point u,point v){ return u.x*v.y-u.y*v.x; }
double solve(int x,int y,int k1,int k2,double t1,double t2){
	point v1=(a[x][k1].p-a[x][k1].v*c[x][k1])-(a[y][k2].p-a[y][k2].v*c[y][k2]),v2=a[x][k1].v-a[y][k2].v;
	double u=dot(v2,v2),v=2*dot(v1,v2),w=dot(v1,v1),t;
	if (!u){
		if (v>0) t=t1; else t=t2;
		return sqrt(v*t+w);
	} else{
		t=-v/(2*u);
		if (t<t1) t=t1; if (t>t2) t=t2;
		return sqrt(u*t*t+v*t+w);
	}
}
int main(){
	scanf("%lf%lf%lf",&o.x,&o.y,&r);
	scanf("%d%d",&n,&m); int i,j;
	for (i=1; i<=n; i++)
		scanf("%lf%lf%lf%lf",&a[i][0].p.x,&a[i][0].p.y,&a[i][0].v.x,&a[i][0].v.y);
	double t,u,v,w,ans=1e100; point nml,p,q;
	for (i=1; i<=n; i++)
		for (j=1; j<=m; j++){
			p=a[i][j-1].p-o; q=a[i][j-1].v;
			u=dot(q,q); v=2*dot(p,q); w=dot(p,p)-r*r;
			t=(sqrt(v*v-4*u*w)-v)/u/2; 
			c[i][j]=c[i][j-1]+t;
			a[i][j].p=a[i][j-1].p+a[i][j-1].v*t;
			nml=a[i][j].p-o; swap(nml.x,nml.y); nml.x=-nml.x;
			a[i][j].v=nml*(dot(nml,a[i][j-1].v)/dot(nml,nml)*2)-a[i][j-1].v;
		}
	for (i=1; i<n; i++)
		for (j=i+1; j<=n; j++){
			int k1=0,k2=0;
			while (k1<m && k2<m){
				ans=min(ans,solve(i,j,k1,k2,max(c[i][k1],c[j][k2]),min(c[i][k1+1],c[j][k2+1])));
				if (c[i][k1+1]<c[j][k2+1]) k1++; else k2++;
			}
		}
	printf("%.3f\n",ans);
	return 0;
}


by lych

2016.3.3

展开阅读全文

没有更多推荐了,返回首页