题解 P4682 【[ZJOI2007]粒子运动】
liaohaoping · · 题解
这题数据范围并不大，显然不用什么复杂的优化，求出每个粒子每段（这里简称某粒子两次碰撞之间的部分为一段，包括时刻0到第一次碰撞）的信息就行了，通过观察可发现，只需要处理好以下三个问题（为了方便，我们把所有坐标分别减去圆心的坐标，再把圆心的坐标改为(0,0)）：
1、某粒子在某段与圆下次碰撞的位置
2、某粒子在某次碰撞后vx,vy的变化
3、粒子a的某一段与粒子b的某一段（保证某一时刻这两段同时存在）的历史最近距离
先看问题1，若此段起始点为(stx,sty)，速度为(vx,vy)，设碰撞时间为t，则(stx+vx·t)²+(sty+vy·t)²=r² → stx²+vx²·t²+2·stx·vx·t+sty²+vy²·t²+2·sty·vy·t-r²=0 → (vx²+vy²)·t²+(2·stx·vx+2·sty·vy)·t+(stx²+sty²-r²)=0，我们设a=vx²+vy²,b=2·stx·vx+2·sty·vy,c=stx²+sty²-r²，则t=-b+sqrt(b²-4ac)/2a（舍去小解，小解一定不成立），再通过t算出位置即可。
问题2也不难想，若与圆的碰撞点为(x,y)，原速为(x1,x2)，反弹后为(x2,y2)，不难得出atan2(x2,y2)-atan2(x,y)=atan2(x,y)-atan2(x1,y1)（角相等），即atan2(x2,y2)=2·atan2(x,y)-atan2(x1,y1)，设k=x2/y2，则k=tan(2·atan2(x,y)-atan2(x1,y1))，因为速度相同，所以x1²+y1²=x2²+y2²，设z=x1²+y1²，则y2²+k²·y2²=z → (k²+1)·y2²=z → y2=sqrt(x1²+y1²/k²+1)，然后通过y2求出x2，需要注意的是，(-x2,-y2)在这里同样成立，我们可以通过带入一个点判断是否在圆内来判断取哪个解。
问题3稍微难一点，为了方便，我们设tt=max(t1,t2)（t1,t2分别为两段开始的时刻），若(x1,y1)，(x2,y2)分别为tt时刻两粒子的坐标，(vx1,vy1)，(vx2,vy2)分别为两段的速度，y为两粒子间距离的平方，y=((x1+vx1·t)-(x2+vx2·t))²+((y1+vy1·t)-(y2+vy2·t))² → y=((x1-x2)+(vx1-vx2)·t)²+((y1-y2)+(vy1-vy2)·t)²我们设xx=x1-x2,yy=y1-y2,vx=vx1-vx2,vy=vy1-vy2，那么y=(xx+vx·t)²+(yy+vy·t)² → y=xx²+vx²·t²+2·xx·vx·t+yy²+vy²·t²+2·yy·vy·t → y=(vx²+vy²)·t²+(2·xx·vx+2·yy·vy)·t+(xx²+yy²)，再设a=vx²+vy²,b=2·xx·vx+2·yy·vy,c=xx²+yy²，然后求出最小值y3=4ac-b²/4a即为历史最近距离，特殊情况特判一下就行了。
我们一开始先用1、2点处理出每个粒子的每一段，再两个两个粒子求历史最近距离，找最小值即为答案。
代码
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read()
{
int x=0,w=0;char ch=getchar();
while(ch<'0'||ch>'9')w|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
return w?-x:x;
}
struct line{double stx,sty,vx,vy,edt;};
struct particle{line f[110];}p[110];
int n,k;double r,z=0.000000001;
inline void get_p(int i,int j,double stx,double sty,double vx,double vy)
{
p[i].f[j].stx=stx,p[i].f[j].sty=sty,p[i].f[j].vx=vx,p[i].f[j].vy=vy;
double a=vx*vx+vy*vy,b=2*stx*vx+2*sty*vy,c=stx*stx+sty*sty-r*r;
double t=(-b+sqrt(b*b-4*a*c))/(2*a);
p[i].f[j].edt=p[i].f[j-1].edt+t;//求碰撞点
if(j<k)
{
double edx=stx+vx*t,edy=sty+vy*t;
double kk=tan(atan2(edx,edy)*2-atan2(vx,vy));
double vy2=sqrt((vx*vx+vy*vy)/(kk*kk+1)),vx2=kk*vy2;
double tt=sqrt(r)/100000/((vx2+vy2)/2);tt=tt<0?-tt:tt;
double xx=(edx+vx2*tt),yy=(edy+vy2*tt);
if(xx*xx+yy*yy>=r*r)vx2=-vx2,vy2=-vy2;
get_p(i,j+1,edx,edy,vx2,vy2);//反弹
}
}
inline double dis(double x1,double y1,double x2,double y2)
{
if(x1==x2&&y1==y2)return 0;
double x=x1-x2,y=y1-y2;
return sqrt(x*x+y*y);
}
int main()
{
register double xx,yy;register int i,j,ii,jj;
scanf("%lf%lf%lf",&xx,&yy,&r);n=read(),k=read();
for(i=1;i<=n;++i)
{
double x,y,vx,vy;scanf("%lf%lf%lf%lf",&x,&y,&vx,&vy);
x-=xx;y-=yy;p[i].f[0].edt=0;get_p(i,1,x,y,vx,vy);
}
double ans=1e50;
for(i=2;i<=n;++i)
{
for(j=1;j<i;++j)
{
jj=1;
for(ii=1;ii<=k;++ii)//第一个粒子处于第ii段
{
--jj;
while(p[i].f[ii].edt>p[j].f[jj].edt&&jj<k)//判断为了使有某时刻两段（ii，jj+1）同时存在，ii是否需要++
{
++jj;
double t1=p[i].f[ii-1].edt,t2=p[j].f[jj-1].edt;
double x1=p[i].f[ii].stx,x2=p[j].f[jj].stx,y1=p[i].f[ii].sty,y2=p[j].f[jj].sty;
double vx1=p[i].f[ii].vx,vx2=p[j].f[jj].vx,vy1=p[i].f[ii].vy,vy2=p[j].f[jj].vy;
if(vx1==vx2&&vy1==vy2)//如果速度，方向相等，进行特殊处理
{
double t=max(t1,t2),xx1=x1+vx1*(t-t1),xx2=x2+vx2*(t-t2),yy1=y1+vy1*(t-t1),yy2=y2+vy2*(t-t2);
ans=min(ans,dis(xx1,yy1,xx2,yy2));
}
else
{
double xx1=x1+vx1*(t1<t2?t2-t1:0),xx2=x2+vx2*(t1>t2?t1-t2:0),yy1=y1+vy1*(t1<t2?t2-t1:0),yy2=y2+vy2*(t1>t2?t1-t2:0);
double vx=vx1-vx2,vy=vy1-vy2;xx=xx1-xx2,yy=yy1-yy2;
double a=vx*vx+vy*vy,b=2*xx*vx+2*yy*vy,c=xx*xx+yy*yy;
double t=-(b/(2*a));
if(t<0)ans=min(ans,dis(xx1,yy1,xx2,yy2));//如果t<0，取最近的合法值(t=0)
else if(t+max(t1,t2)>=p[i].f[ii].edt+z||t+max(t1,t2)>=p[j].f[jj].edt+z)//如果该时刻超出任意一段范围，取最近的合法值
{
t=min(p[i].f[ii].edt,p[j].f[jj].edt),xx1=x1+vx1*(t-t1),xx2=x2+vx2*(t-t2),yy1=y1+vy1*(t-t1),yy2=y2+vy2*(t-t2);
ans=min(ans,dis(xx1,yy1,xx2,yy2));
}
else//否则正常计算
{
double d=sqrt((4*a*c-b*b)/(4*a));
ans=min(ans,d);
}
}//求距离
}
}
}
}
printf("%.3lf",ans);
return 0;
}