闲话 22.12.19

闲话

今天闲话水点就水点吧
也习惯了不是吗

今天在随机歌 看到情绪的一首歌《无法停止的白情》
曲调好熟悉 好听！作词作曲编曲：春卷饭
彳亍 已经在循环了

今日 SAM！

一阶微分方程

看到了 duls 的《A problem collection of ODE and differential technique》
深有感触 指一早读才看懂一阶微分方程怎么个解法

具体地，给定 \(f(t)\)，我们需要求出下式的解 \(x(t)\)：

\[\frac{\text{d}}{\text dt} x = f(x) \]

类似多项式牛顿迭代的做法，我们考虑倍增求解。令 \(x_n\) 为 \(x(t)\) 在模 \(t^n\) 意义下的值。

初值为 \(x_0 = 1\)，现在需要从 \(x_n\) 倍增出 \(x_{2n}\)。考虑将 \(f\) 在\(x_n\) 处展开，取前两项。（第三项及以后带系数 \((x_{2n} - x_n)^k\)，最低次项都达到了 \(t^{2n}\)）

\[\begin{aligned} \frac{\text{d}}{\text dt} x_{2n} & \equiv f(x_n) + f'(x_n)\left(x_{2n} - x_n\right) \pmod{t^{2n}} \\ & \equiv f(x_n) + f'(x_n)x_{2n} - f'(x_n)x_n \pmod{t^{2n}} \end{aligned}\]

考虑构造一个消元多项式

\[p(t) = \exp\left(\int - f'(x_n) \text dt \right) \]

\[p'(t) = - f'(x_n) p(t) \]

为什么这么叫呢？你看 \(\equiv\) 右边是不是有个 \(x_{2n}\) 项吗？我们两边乘 \(p\)，再加上 \(p'x_{2n}\) 凑一下试试？

\[\begin{aligned} p\times \frac{\text{d}}{\text dt} x_{2n} & \equiv (f(x_n) - f'(x_n)x_n)p + f'(x_n)x_{2n}p \pmod{t^{2n}} \\ \frac{\text{d}}{\text dt} (x_{2n}p) & \equiv (f(x_n) - f'(x_n)x_n)p + f'(x_n)x_{2n}p + x_{2n}p' \pmod{t^{2n}} \\ \frac{\text{d}}{\text dt} (x_{2n}p) & \equiv (f(x_n) - f'(x_n)x_n)p + f'(x_n)x_{2n}p - f'(x_n)x_{2n}p' \pmod{t^{2n}} \\ \frac{\text{d}}{\text dt} (x_{2n}p) & \equiv (f(x_n) - f'(x_n)x_n)p \pmod{t^{2n}} \\ x_{2n}p & \equiv \int \left((f(x_n) - f'(x_n)x_n)p\right) \text d t + 1 \pmod{t^{2n}} \\ x_{2n} & \equiv \frac{1} p \left(\int \left((f(x_n) - f'(x_n)x_n)p\right) \text d t + 1 \right) \pmod{t^{2n}} \\ \end{aligned}\]

这样就可以做迭代了。

模板题

这题的 \(f(x) = A\times e^{x - 1} + B\)。需要注意的是，\(f'(x) = A \times e^{x - 1}\)。

时间复杂度是 \(T(n) = T(\frac n2) + O(n\log n) = O(n\log n)\) 的。就是常数有点大。

code

int n;
poly A, B, F;
poly f(poly x_n) {
	return A * (x_n - 1).exp() + B;
} 

poly f2(poly x_n) {
	return A * (x_n - 1).exp();
}

poly r(poly x_n) {
	return (- f2(x_n)).intg().exp();
}

signed main() {
	iot;
	cin >> n; A.redegree(n); B.redegree(n);
	rep(i,0,n) cin >> A[i]; rep(i,0,n) cin >> B[i];
	F.resize(1); F[0] = 1;
	for (int i = 2; i <= (n << 1); i <<= 1) {
		F.resize(i);
		F = r(F).inv() * ((r(F) * (f(F) - f2(F) * F)).intg() + 1);
		F.resize(i);
	} F.redegree(n);
	cout << F << '\n';
}

可以常数变易法解出通解来，常数更小点。

code

int n;
poly A, B, A1, B1, expB1;

signed main() {
	iot;
	cin >> n;
	A.redegree(n); B.redegree(n);
	rep(i, 0, n) cin >> A[i]; rep(i, 0, n) cin >> B[i];
	A1 = A.intg().slice(n), B1 = B.intg().slice(n), expB1 = B1.exp();
	cout << (1 + B1 - (((A1 * B).slice(n - 1) * expB1).slice(n - 1).intg() - (A1 * expB1).slice(n) + 1).ln()) << '\n';
}

达成成就：\(14s \to 600ms\)。