ZOJ Problem Set – 3409
KKV
Time Limit: 10 Seconds Memory Limit: 65536 KB Special Judge
原题链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3409
算法：模拟题、高中动量守恒公式套上三维模板、计算机模拟手解一元二次方程
题目大意：你从(0,0,0)出发，途中经过N-1个点，最后到达第N个点。(N<=40)
这些点均在三维空间中，这N个 点的三维坐标均给出。题目背景是太空中的宇宙飞行，故无重力影响，每个点只靠火箭助推器的喷射来改变速度大小和方向。（而且每次动量守恒求解过程中，保证 助推器的速度大小均达到最大。速度根据情况而定）故每个拐点的运动改变方式完全符合动量守恒定律（只不过是三维空间的）。
现给你M，总质量，m1，不含所有燃料的火箭剩余质量（1<=M1<M-n*m0）,m0，单块燃料（火箭助推器）的质量，v0，每次喷射时火箭助推器的速度。N个点的坐标。
现求解，从第0个点到第N-1个点都有火箭助推器参与的速度改变过程，故问到第N点时所花费的总时间为多少。若无解，则输出无解时对应的报错语句。
现给出代码及注释：
//解法在注释中给出，分别求解一个一元一次方程和一元二次方程即可
/*
还是在之前模拟一下方程吧。防止看不懂
动量守恒定律基本方程。
m1v1+m2v2=m1v1’+m2v2’
此处分别有两次应用（实际上可以规约为一种应用，但是怕挫，就分成两种了）
第一种是初始点0=(m-m0) * v(x,y,z)+m0 * v0(x,y,z)//v和v0分别为三维速度矢量，v为第一次喷射后剩余火箭箭体的速度矢量，v0为第一次喷射燃料时燃料的速度矢量（所以这其实是三个方程联立，解法为同时两边平方再相加）
式中各参量分别满足v0(x)^2+v0(y)^2+v0(z)^2=v0^2//v0已知
v(x,y,z) = k * vec(x,y,z)//vec(x,y,z)为第0点到第一个点的位移矢量
解得K^2= M0^2*V0^2/(Vec^2*(m-m0)^2), k= sqrt(K^2)
所有值均已知，代入即可解得//vec为vec矢量的模长，v0为标量
第二种为一般情况，针对第i个点（1<=I<=n-1）
M’v’(x,y,z)=(m’-m0)v(x,y,z)+m0v0(x,y,z)//实际上三方程联立
v’即为第i-1个点的速度矢量，v为第I个点的速度矢量，m’为第i-1个点剩余的火箭质量
v(x,y,z)=k * vec(x,y,z)//vec为第i个点到第I+1个点的位移矢量
v0x^2+v0Y^2+v0z^2=|v0|^2
将联立方程组左右平方并加和将常数项移到一侧，从而转化为求解k为主元的一元二次方程问题，化归式如下：
(m’- m0) ^ 2 * |vec| ^ 2 * k ^ 2 – 2 * M’ * (m’- m0) * scalar(v’,vec) * k + m’ ^ 2 * |v’| ^ 2 - m0 ^ 2 * |V0| ^ 2 = 0
一元二次方参数：
A=(M’-m0)^2 * |vec| ^ 2
B=-2*M’*(m’-m0)*scalar(v’,vec)
C=M’^2*|v’|^2-M0^2*|v0|^2
det = b^2 – 4 * a * c
判断det有解与否即可
有解则可由v(x,y,z)=k*vec(x,y,z)知道v的矢量及大小
从而算出ti = |vec| / |v|
将第一段与之后的n-1段的时间即可
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double eps = 1e-10;
inline double sqr(double x){return x * x;}
struct point
{
double x,y,z;
void input()
{
scanf("%lf%lf%lf",&x,&y,&z);
}
}p[50];
struct vect//矢量定义
{
double x,y,z,mod;
void init(point p,point q) //矢量初始化并计算模长
{
x = q.x - p.x;
y = q.y - p.y;
z = q.z - p.z;
mod = sqrt(sqr(x) + sqr(y) + sqr(z));
}
double cal()//对于速度矢量进行计算速率
{
return mod = sqrt(sqr(x) + sqr(y) + sqr(z));
}
}vec[50],v[50]; //位移矢量和速度矢量
double scalar(vect a,vect b) //点积
{
return a.x * b.x + a.y * b.y + a.z * b.z;
}
void function(double a,double b,double c,bool &mark,double &k1,double &k2)
{ //求解二次方程专用函数
double det = sqr(b) - 4.0 * a * c; //delta
if(det < eps && fabs(det) > eps) //delta<0则无解，无解即对应题目中无解情况
{
mark = false;
return ;
}
else if(fabs(det) < eps) //唯一解
{
mark = true;
k1 = k2 = -b / (2.0 * a);
return ;
}
else //二解的话一正解一负解舍弃负解即可（貌似最后算出来delta都等于0？我没测试，大概吧。。。我甚至没去手工模拟一下负解所对应的实际情况是怎样的，就直接否决了，因为负解一定不会符合题意，根据我方程的参数即可知道）
{
mark = true;
k1 = (-b - sqrt(det)) / (2.0 * a);
k2 = (-b + sqrt(det)) / (2.0 * a);
return ;
}
}
void solve(int n)
{
double m,m1,mt,m0,v0,k,t = 0.0,temp;
int i;
scanf("%lf%lf%lf%lf",&m,&m1,&m0,&v0);
p[0].x = p[0].y = p[0].z = 0.0;
for(i = 1;i <= n;i++)
{
p[i].input();
vec[i - 1].init(p[i - 1],p[i]); //位移矢量初始化
}
//p[0] //对于原点处速度矢量和第一段距离的时间的求解
/*
0 = (M-m0)*Vx - m0*V0x //其实可以列一个三元一次方程组用行列式解，我嫌麻烦
0 = (M-m0)*Vy - m0*V0y //也可以用纯几何的点积叉积做，我不会~算了~
0 = (M-m0)*Vz - m0*V0z
(Vx,Vy,Vz) = k * (Vecx,Vecy,Vecz) //求解对象即为k，因为速度矢量方向确定，但是大小不确定，解出大小即可
V0x^2 + V0y^2 + V0z^2 = V0^2 //V0作为已知量，要做好未知量的化归转化，才方便求出可行解
t0 = |Vec0| / |V[0]| //标量除法
t += t0
*/
mt = m - m0; //第一次后剩余质量
k = sqrt(sqr(m0) * sqr(v0) / sqr(mt) / sqr(vec[0].mod));
v[0].x = k * vec[0].x;
v[0].y = k * vec[0].y;
v[0].z = k * vec[0].z;
temp = vec[0].mod / v[0].cal();
t += temp;
//对于一般点处，求解一般方程
/*
ti = |Veci| / |V[i]|
Mt * Vx[i - 1] = (Mt - m0) * Vx[i] + m0 * V0x
Mt * Vy[i - 1] = (Mt - m0) * Vy[i] + m0 * V0y
Mt * Vz[i - 1] = (Mt - m0) * Vz[i] + m0 * V0z
V[i](x,y,z) = k * Vec[i](x,y,z)
V0x^2 + V0y^2 + V0z^2 = V0^2
k^2*|Vec[i]|^2*(Mt-m0)^2 - 2*Mt(Mt-m0)*k*(scalar(Vec[i],V[i-1])) + Mt^2*|V[i-1]|^2 - m0^2*|V0|^2 = 0
*/
for(i = 1;i < n;i++)
{
double k1,k2;
double a,b,c;
bool mark = true;
a = sqr(vec[i].mod) * sqr(mt - m0);
b = -2.0 * mt * (mt - m0) * scalar(vec[i],v[i - 1]);
c = sqr(mt) * sqr(v[i - 1].cal()) - sqr(m0) * sqr(v0);
function(a,b,c,mark,k1,k2);
if(mark)
{
if(k1 > 0)
k = k1;
if(k2 > 0)
k = k2;
v[i].x = k * vec[i].x;
v[i].y = k * vec[i].y;
v[i].z = k * vec[i].z;
temp = vec[i].mod / v[i].cal();
t += temp;
}
else
{
printf("Another kind of KKV.\n");
return ;
}
mt = mt - m0;
}
printf("%.2lf\n",t);
}
int main()
{
//freopen("3409.in","r",stdin);
int n;
while(scanf("%d",&n) != EOF)
{
solve(n);
}
}