首页 > 技术文章 > U179915 关于分级火箭的一点理想化的计算

Square-Circle 2021-10-09 21:46 原文

题目地址

本题是一道疯狂推式子的玄学复杂度sb题。

解题思路

1.数学部分

​ 首先假定已经将火箭分成了 \(n+1\) 级,记使用了 \(n\) 个分级器。记各级的开始时间点为:

\[0=t_0<t_1<t_2<\cdots<t_n<t_{n+1}=T \]

其中 \(t_0\) 总的开始,\(t_{n+1}\) 为总的结束。称时间段 \([t_k,t_{k+1}]\) 为第 \(k\) 阶段, \(0\leqslant k \leqslant n\)

​ 为了方便,记第 \(k\) 阶段除燃料之外的质量总共为:

\[M_{\times}(k)=M_0+\frac{n+1-k}{n+1}M_1+(n-k)M_2 \]

设火箭该阶段燃料质量关于时间的函数为 \(m_k(t)\) ,根据“火箭瞬时的燃料消耗速度与火箭的瞬时总质量成正比”,以及每一阶段结束时剩余的燃料质量,列出方程:

\[\left\{ \begin{align*} & \frac{\text{d}m_k}{\text{d}t}=-\frac{1}{\eta}(m_k+M_{\times}(k)) \\ & m_k(t_{k+1})=\frac{n-k}{n+1}M \end{align*} \right. \]

解得:

\[m_k(t)=\left(M_{\times}(k)+\frac{n-k}{n+1}M\right){e}^{\frac{t_{k+1}-t}{\eta}}-M_{\times}(k) \]

​ 考虑第 \(k\) 阶段开头,有:

\[m_k(t_k)=\frac{n+1-k}{n+1}M \]

带入化简可得:

\[\frac{t_{k+1}-t_k}{\eta} =\ln{\left(M_{\times}(k)+\frac{n+1-k}{n+1}M\right)}-\ln{\left(M_{\times}(k)+\frac{n-k}{n+1}M\right)} \]

累加,得到:

\[\frac{T}{\eta}=\frac{t_{n+1}-t_0}{\eta}=\sum_{k=0}^{n}{\frac{t_{k+1}-t_k}{\eta}} =\sum_{k=0}^{n}{\left(\ln{\left(M_{\times}(k)+\frac{n+1-k}{n+1}M\right)}-\ln{\left(M_{\times}(k)+\frac{n-k}{n+1}M\right)}\right)} \]

2.代码部分

​ 根据以上推导,我们列出了关于未知数 \(M\) 的方程:

\[\sum_{k=0}^{n}{\left(\ln{\left(M_{\times}(k)+\frac{n+1-k}{n+1}M\right)}-\ln{\left(M_{\times}(k)+\frac{n-k}{n+1}M\right)}\right)}=\frac{T}{\eta} \]

由于题目仅要求整数级别的复杂度,我们可以通过二分(通过“简单“的求导就能证明其单调性)近似求解 \(M\) 的值,每次通过循环计算出上式中等号左侧的值,若不小于右侧则合法。其复杂度为 \(n\log M\)

​ 然后我们需要对 \(n\) 进行遍历。若采取直接遍历的方法,总复杂度为 \(n^2\log M\) ,会TLE。注意到 \(M\) 关于 \(n\) 的函数先减后增,于是可以采用爬山算法,随机地进行求解。

代码\(_{_{\text{(数据就是它造的)}}}\)

#include<iostream>
#include<cmath>
using namespace std;
typedef unsigned long long ull;
int M0, M1, M2, T, eta;		//五个参数
const int steps = 5, step[steps] = { 10000,1000,100,10,1 }, maxn = 100000;
		//步数大小和最大层数,这会影响算法正确性,请“适当”地自行选择
constexpr double Mx(int n, int k) {
	return (double)(n + 1 - k) * M1 / (n + 1) + M0 + (double)(n - k) * M2;
}	//第k阶段除燃料之外的质量
bool check(int n, ull m) {
	double sum = 0;
	for (int k = 0; k <= n; k++)
		sum += log(Mx(n, k) + (double)(n + 1 - k) / (n + 1) * m)
		- log(Mx(n, k) + (double)(n - k) / (n + 1) * m);
	if (sum < (double)T / eta) return false;
	return true;
}	//检测合法性
int main() {
	cin >> M1 >> M2 >> eta >> M0 >> T;
	int ansn = 0;
	ull ansm = (ull)ceil((exp((double)T / eta) - 1) * ((double)M0 + M1));
	for (int i = 0; i < steps; i++) {		//枚举步数大小
		for (int n = ansn + step[i]; n < maxn; n += step[i]) {
			ull lm = (ull)floor(Mx(n, 0) * eta), rm = ansm;
			if (lm > rm || !check(n, rm)) break;
			while (lm < rm) {
				ull mid = (lm + rm) / 2;
				if (check(n, mid)) rm = mid;
				else lm = mid + 1;
			}		//易懂的二分答案
			if (rm < ansm) {		//为使n尽量小,正向时是小于
				ansm = rm; ansn = n;
			}
		}		//正向爬山
		for (int n = ansn - step[i]; n < maxn; n -= step[i]) {
			ull lm = (ull)floor(Mx(n, 0) * eta), rm = ansm;
			if (lm > rm || !check(n, rm)) break;
			while (lm < rm) {
				ull mid = (lm + rm) / 2;
				if (check(n, mid)) rm = mid;
				else lm = mid + 1;
			}
			if (rm <= ansm) {		//为使n尽量小,反向时是小于等于
				ansm = rm; ansn = n;
			}
		}		//反向爬山
	}
	cout << ansn << endl;
	return 0;
}

推荐阅读