注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

n+e

NewWeb:http://trinkle.is-programmer.com/

 
 
 

日志

 
 

[BZOJ1094][ZJOI2007]粒子运动  

2014-07-21 22:54:07|  分类: BZOJ |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

Description

阿Q博士正在观察一个圆形器皿中的粒子运动。不妨建立一个平面直角坐标系,圆形器皿的圆心坐标为(x0, y0),半径为R。器皿中有若干个粒子,假设第i个粒子在时刻0的位置为(xi, yi),速度为(vxi,vyi)(注:这是一个速度向量,若没有发生碰撞,t时刻的位置应该是(xi + t * vxi, yi + t * vyi) )。假设所有粒子的运动互不干扰;若某个粒子在某个时刻碰到了器皿壁,将发生完全弹性碰撞,即速度方向按照碰撞点的切线镜面反射,且速度大小不变(如图)。认为碰撞是瞬间完成的。
 [BZOJ1094][ZJOI2007]粒子运动 - Trinkle - n+e
尽管碰撞不会影响粒子的速率,但是粒子却会受到一定的伤害,所以若某一个粒子碰撞了k次器皿壁,那么在第k次碰撞时它便会消亡。出于研究的需要,阿Q博士希望知道从时刻0到所有粒子都消亡这段时间内,所有粒子之间的最近距离是什么。你能帮助他么?

Input

输入文件particle.in第一行包含三个实数,分别为x0, y0, R,即圆形器皿的圆心坐标及半径。第二行包含两个正整数N, k,分别表示粒子的总数与消亡碰撞次数。接下来N行每行四个实数,分别为xi, yi, vxi , vyi,保证(xi, yi)都在圆内且(vxi, vyi)非零。

Output

仅包含一个实数,即所有粒子的历史最近距离,精确到小数点后三位。

Sample Input

0 0 10
2 10
0 -5 0 1
5 0 1 0

Sample Output

7.071

HINT

对于所有的数据,2 ≤N ≤100。1≤k ≤100。
请注意实数精度问题。

Solution

好久以前写的题了,好像是上学期寒假。挖出来看了看,把之前的底层代码优化了一下,快了不少。

果然手写才是硬道理= =

题目大意:有若干个粒子在圆形器皿内运动,碰撞器皿k次后停止运动,求所有粒子间的历史最近距离。

算法分析:显然我们可以把求一堆粒子中的历史最近距离减少为求两个粒子的历史最近距离。让我们来推式子吧~

以下要做三件事情:

①已知当前点Px,y)和vp,q),求它会撞到圆上的点B

[BZOJ1094][ZJOI2007]粒子运动 - Trinkle - n+e
 


B=P+λv.(λ>0)B在圆上 ∴XB+YB=R^2 即 (x+λp)^2+(y+λq)^2=R^2

展开得λ^2*(p^2+q^2)+2λ(xp+yq)+x^2+y^2-R^2=0

是一个关于λ的一元二次方程,代入求根公式求正根即可,其中

a=dot(v,v),b=2dot(v,p),c=dot(P,P)-R^2


②已知Bv,求碰撞后的v

[BZOJ1094][ZJOI2007]粒子运动 - Trinkle - n+e

为了方便,我们在复平面上定义点和向量,这样就能够通过向量乘向量、向量除向量来旋转角度,则点B的切线向量为B*i。由复数的定义,我们可以得到v=(B*i)^2/v

注意,碰撞前后速度相同,所以length(v)=length(v)

 

③已知点P的速度为u,点Q的速度为v,它们保持这样运动的时间为time0,求距离的最小值。

[BZOJ1094][ZJOI2007]粒子运动 - Trinkle - n+e

设经过时间t后,|PtQt|=d.

Pt=(Px+tux,Py+tuy),Qt=(Qx+tvx,Qy+tvy)

d=sqrt(t^2[(ux-vx)^2+(uy-vy)^2] + 2t[(ux-vx)(Px-Qx)+(uy-vy)(Py-Qy)] +(Px-Qx)^2+(Py-Qy)^2 )

除去sqrt不看,这是一个关于t的二次函数,其中

a=dot(u-v,u-v),b=2*dot(u-v,p-q),c=dot(p-q,p-q)

0-b/2atime0时,能取到最小值c-b^2/4a

 

有了这三个函数,我们就可以很方便地跟踪两个点在一段连续时刻的运动轨迹,不断更新答案就是最终答案。

复杂度:O(n^2*k)

注意输入数据除n,k外均是实数(就因为这我在横八上wa了五六次,还以为是精度问题......


/**************************************************************
    Problem: 1094
    User: wjy1998
    Language: C++
    Result: Accepted
    Time:508 ms
    Memory:820 kb
****************************************************************/
 
#include<cstdio>
#include<cmath>
typedef double ld;
int n,k;ld ans=1e10,X0,Y0,R;
inline ld min(const ld&x,const ld&y){return x<y?x:y;}
struct P{ld x,y;}I=(P){0,1};
P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
P operator*(const P&a,const ld&p){return (P){a.x*p,a.y*p};}
P operator*(const P&a,const P&b){double d=a.x*b.x,e=a.y*b.y,f=(a.x+a.y)*(b.x+b.y);return (P){d-e,f-e-d};}
P operator/(const P&a,const ld&p){return (P){a.x/p,a.y/p};}
P operator/(const P&a,const P&b){double x=b.x*b.x+b.y*b.y;return (P){(a.x*b.x+a.y*b.y)/x,(a.y*b.x-a.x*b.y)/x};}
ld dot(const P&a,const P&b){return a.x*b.x+a.y*b.y;}
ld cross(const P&a,const P&b){return a.x*b.y-a.y*b.x;}
ld length(const P&a){return sqrt(dot(a,a));}
ld solve_x2(ld a,ld b,ld c){return(-b+sqrt(b*b-4*a*c))/(a+a);}
ld lim_x2(ld a,ld b,ld c,ld t){return min(min(a*t*t+b*t+c,c),(2*a*t+b>=0&&a*b<0?(c-b*b/(4*a)):1e10));}
ld get_point(const P&p,const P&v){return solve_x2(dot(v,v),2*dot(v,p),dot(p,p)-R*R);}
P set(P a,ld len){ld length=len/sqrt(a.x*a.x+a.y*a.y);a.x*=length;a.y*=length;return a;}
P make(P B,const P&v){B=B*I;return set(B*B/v,length(v));}
struct C{P p,v;}c[110];
double solve(int c1,int c2){
    P p=c[c1].p,q=c[c2].p,u=c[c1].v,v=c[c2].v;
    double time,ans=1e10,lmd1,lmd2;
    for(int t1=1,t2=1;t1<=k&&t2<=k;)
    {
        lmd1=get_point(p,u),lmd2=get_point(q,v),time=min(lmd1,lmd2);
        ans=min(ans,lim_x2(dot(u-v,u-v),2*dot(u-v,p-q),dot(p-q,p-q),time));
        if(time==lmd1)t1++,p=p+u*lmd1,u=make(p,u),q=q+v*lmd1;
        else t2++,p=p+u*lmd2,q=q+v*lmd2,v=make(q,v);
    }
    return ans;
}
int main(){
    int i,j;double x,y,vx,vy;
    scanf("%lf%lf%lf%d%d",&X0,&Y0,&R,&n,&k);
    for(i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf%lf",&x,&y,&vx,&vy),x-=X0,y-=Y0;
        c[i].p=(P){x,y},c[i].v=(P){vx,vy};
    }
    for(i=1;i<n;i++)
    for(j=i+1;j<=n;j++)ans=min(ans,solve(i,j));
    if(ans<0)ans=0;printf("%.3lf\n",sqrt(ans));
}
  评论这张
 
阅读(452)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018