[2020HDU多校第一场][HDU 6751][A. Avian Darts]

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=6751

题目大意：有一初始在原点的动点(P)，初始速度方向为((1,0,0))，该动点在运动时会同时携带它的标架进行移动（看题目链接中的样例解释应该能够理解），求在(n)次运动后点(P)所在的位置以及标架的方向。每次运动，数据会给出在这段运动中，点(P)在各个方向的角加速度、移动速度以及运动时长。

题解：前段时间刚学的四元数，没想到补这题就用上了_(:з」∠)_（虽然不用四元数也能做）

　　　首先通过分析（看样例解释）可以得出，(wx)代表着绕当前标架的(X)轴进行逆时针旋转，(wy)就代表绕(Y)轴旋转，(wz)同理。设点(Q(wx,wy,wz))，那么可以发现点(P)是在绕着(vec{OQ})进行着旋转，角速度为(|vec{OQ}|)，这个结论可以通过微分几何相关知识得出（可惜本弱微分几何差点挂科，学过的也基本忘了，证明完全不会，就当结论记吧）。由于每段运动中还有个沿(x)轴的初速度，于是在每段运动中，可以看成是点(P)在绕着这个轴进行螺旋式前进。

　　　接着我们可以利用学过的高中物理知识进行求解。在每段运动中，点(P)的初始运动方向为(vec{X})，即当前标架(X)轴方向，速度为(v)，则在不旋转的情况下，(P)的速度为(vec{V}=vcdot vec{X})。由于(P)在前进的同时需要绕(vec{OQ})旋转，因此我们可以设(vec{V})在(vec{OQ})上的投影为(vec{V1})，并令(vec{V2}=vec{V}-vec{V1})，这样点(P)的初速度就被分解成了两个方向，(vec{V1})是与(vec{OQ})平行的向量，而(vec{V2})与旋转轴垂直。因此在(vec{V1})方向上，点(P)是在做着匀速直线运动，而在另一方向，则是在绕着(vec{OQ})做匀速圆周运动。

　　　接下来我们需要求出这个匀速圆周运动的圆心与半径。很明显，这个圆周运动是在垂直于(vec{OQ})的平面上进行的，而高中学过的知识告诉我们，圆周运动时点到圆心连线与运动方向垂直。于是我们可以设这个圆心为(O')，并得出(vec{PO'})与(vec{V2})、(vec{V1})均垂直，不妨设(vec{V3}=vec{PO'})，可以发现(vec{V1})、(vec{V2})、(vec{V3})组成了一个标架，于是我们可以通过求(vec{V1})与(vec{V2})的叉积计算出(vec{V3})。计算出(vec{V3})后，我们就能得到(O')的坐标，并由此计算出这次圆周运动对答案的贡献。这样将两个运动叠加到一起，我们可以得到点(P)在这次运动结束后的坐标。同时，由于标架的运动不需要考虑动点前进的速度，只需要考虑角速度，所以只要将标架的每个方向绕(vec{OQ})进行相应旋转就好了。

　　　而关于绕轴旋转的操作如何实现，如果不会可以借助四元数这一工具，详情可以点击https://codeforces.com/blog/entry/46744进行学习。而每次旋转的角度是多少以及匀速圆周运动的半径具体怎么算本弱也不太清楚，因为本人是对着样例调系数算出来的（草），不过相信聪明的你一定能算出来√

#include<bits/stdc++.h>
using namespace std;
#define LL long long 
#define ld double
const ld eps=1e-9;
const ld pi=acos(-1.0);
struct Point
{
    ld x,y,z;
    //void read(){scanf("%Lf%Lf%Lf",&x,&y,&z);}
    Point operator +(const Point &t)const{return {x+t.x,y+t.y,z+t.z};}
    Point operator -(const Point &t)const{return {x-t.x,y-t.y,z-t.z};}
    Point operator *(const Point &t)const{return {y*t.z-z*t.y,z*t.x-x*t.z,x*t.y-y*t.x};}
    ld operator %(const Point &t)const{return x*t.x+y*t.y+z*t.z;}
    Point operator *(const ld &t)const{return {t*x,t*y,t*z};}
    Point operator /(const ld &t)const{return {x/t,y/t,z/t};}
    ld norm(){return x*x+y*y+z*z;}
    ld len(){return sqrt(norm());}
    //void print(){printf("%.6Lf %.6Lf %.6Lf
",x,y,z);}
    //void print(){cout<<x<<" "<<y<<" "<<z<<endl;}
    void print(){printf("%.9f %.9f %.9f
",x,y,z);}
    void printZ(){printf("%lld %lld %lld
",(LL)round(x),(LL)round(y),(LL)round(z));}
}P,X,Y,Z,dx,dy,dz,V,rev;
struct Q
{
    ld s;Point v;
    void set(const Point &t){s=0,v=t;}
    Q operator *(const ld &t)const{return {s*t,v*t};}
    Q operator /(const ld &t)const{return {s/t,v/t};}
    Q operator +(const Q &t)const{return {s+t.s,v+t.v};}
    Q operator -(const Q &t)const{return {s-t.s,v-t.v};}
    Q operator *(const Q &t)const{return {s*t.s-v%t.v,t.v*s+v*t.s+v*t.v};}
    Q conj(){return {s,v*(-1)};}
    ld norm(){return s*s+v.norm();}
    ld len(){return sqrt(norm());}
    Q inv(){return conj()/norm();}
    Q rotate(Point i,ld phi)
      {
      Q g={cos(phi*0.5),i*sin(phi*0.5)};
      return g*(*this)*g.inv();
      }
    Q rotate(Q i,ld phi)
      {
      Q g={cos(phi*0.5),i.v*sin(phi*0.5)};
      return g*(*this)*g.inv();
      }
    //void print(){printf("%.6Lf ",s);v.print();}
};
Point rotate(Point A,Point B,ld phi)
{
    if(B.len()<eps)return A;
    Q x;
    x.set(A);
    return x.rotate(B/B.len(),phi).v;
}
int n;
ld wx,wy,wz,v,t;
void init()
{
    P={0,0,0};
    X={1,0,0};
    Y={0,1,0};
    Z={0,0,1};
    for(int i=1;i<=n;i++)
      {
      scanf("%lf%lf%lf%lf%lf",&wx,&wy,&wz,&v,&t);
      dx=X*(wx*pi/180.0),dy=Y*(wy*pi/180.0),dz=Z*(wz*pi/180.0);
      V=X*v;rev=dx+dy+dz;
      if(rev.norm()<eps)P=P+V*t;
      else
        {
        X=rotate(X,rev,rev.len()*t);
        Y=rotate(Y,rev,rev.len()*t);
        Z=rotate(Z,rev,rev.len()*t);
        Point V1=rev*((V%rev)/rev.norm());
        Point V2=V-V1;
        Point V3=V2*rev/rev.norm();
        P=(P-V3)+rotate(V3,rev,rev.len()*t)+V1*t;
        }
      }
    P.print();
    X.print();
    Y.print();
    Z.print();
    printf("
");
}
int main()
{
    while(~scanf("%d",&n))init();
}