专栏文章

题解：P12857 [NERC 2020 Online] Hit the Hay

llittlez_meow
P12857题解参与者 1已保存评论 0

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
0 条
当前快照
1 份
快照标识符
@mioy29yo
此快照首次捕获于
2025/12/03 03:01
3 个月前
此快照最后确认于
2025/12/03 03:01
3 个月前
查看原文
算上父母有五种状态：婴儿状态 00、父母入睡婴儿状态 1,21,2 、父母未入睡婴儿状态 1,21,2。我们称为状态 0,1,2,3,40,1,2,3,4。
对于时间离散的情况，我们可以直接设 fi(t)f_i(t) 表示当前处于状态 ii，距离夜晚结束还有 tt 时刻的最优策略下睡眠的期望时间。
但是非常可惜的是，这里时间是连续的。
不过我们还是可以考虑一个最小时间单位 ϵ→0\epsilon\to 0。
以 f1f_1 为例，写出转移：
f1(t+ϵ)=ϵ+p1ϵf1(t)+(1−p1ϵ)(q0f0(t)+(1−q0)f2(t))f_1(t+\epsilon)=\epsilon+p_1^{\epsilon}f_1(t)+(1-p_1^{\epsilon})(q_0f_0(t)+(1-q_0)f_2(t))
等式两边同时减掉 f1(t)f_1(t)：
f1(t+ϵ)−f1(t)=ϵ+(p1ϵ−1)f1(t)+(1−p1ϵ)(q0f0(t)+(1−q0)f2(t))f_1(t+\epsilon)-f_1(t)=\epsilon+(p_1^{\epsilon}-1)f_1(t)+(1-p_1^{\epsilon})(q_0f_0(t)+(1-q_0)f_2(t))
再同时除掉 ϵ\epsilon：
f1(t+ϵ)−f1(t)ϵ=1+p1ϵ−1ϵ−p1ϵ−1ϵ(q0f0(t)+(1−q0)f2(t))\dfrac{f_1(t+\epsilon)-f_1(t)}{\epsilon}=1+\dfrac{p_1^{\epsilon}-1}{\epsilon}-\dfrac{p_1^{\epsilon}-1}{\epsilon}(q_0f_0(t)+(1-q_0)f_2(t))
左边是导数定义，右边的 p1ϵ−1ϵ\dfrac{p_1^\epsilon-1}{\epsilon} 相当于是函数 g(x)=p1xg(x)=p_1^{x} 在 x=0x=0 处的导数 g′(0)=ln⁡p1g'(0)=\ln p_1。带入化简得到：
f1′(t)=1+f1(t)ln⁡p1−(q0f0(t)+(1−q0)f2(t))ln⁡p1f_1'(t)=1+f_1(t)\ln p_1-(q_0f_0(t)+(1-q_0)f_2(t))\ln p_1
同理我们还可以写出 f0,f2f_0,f_2 的转移：
f0′(t)=f0(t)ln⁡p0−f3(t)ln⁡p0f_0'(t)=f_0(t)\ln p_0-f_3(t)\ln p_0
f2′(t)=1+f2(t)ln⁡p2−f1(t)ln⁡p2f_2'(t)=1+f_2(t)\ln p_2-f_1(t)\ln p_2
然后是状态 44。现在已经是最佳状态，如果现在不决策入睡其他状态就更不会睡了，所以：
f4(t)=f2(t)f_4(t)=f_2(t)
最麻烦的是状态 33，因为可以选择是否入睡。具体地：
f3(t+ϵ)=max⁡{f1(t+ϵ),p1ϵf3(t)+(1−p1ϵ)(q0f3(t)+(1−q0)f4(t))}f_3(t+\epsilon)=\max\{f_1(t+\epsilon),p_1^\epsilon f_3(t)+(1-p_1^\epsilon)(q_0f_3(t)+(1-q_0)f_4(t))\}
先研究取到第一项的情况，即 f3(t)=f1(t)f_3(t)=f_1(t)。
那么上式即：
f3′(t)=max⁡{f1′(t),f1(t)ln⁡p1−(q0f1(t)+(1−q0)f2(t))ln⁡p1}f_3'(t)=\max\{f_1'(t),f_1(t)\ln p_1-(q_0f_1(t)+(1-q_0)f_2(t))\ln p_1\}
把第一项减去第二项得到 1+q0(f1(t)−f0(t))ln⁡p11+q_0(f_1(t)-f_0(t))\ln p_1。如果取到第二项，则该式 ≤0\le 0，即：
−q0(f1(t)−f0(t))ln⁡p1≥1-q_0(f_1(t)-f_0(t))\ln p_1\ge 1
感性上，f1(t)−f0(t)f_1(t)-f_0(t) 应该是单调递增的。所以在解出来 f3(t)=f1(t)f_3(t)=f_1(t) 的微分方程组之后，我们可以二分出满足上式最小的 tt，然后在 tt 大于这个值的时候解出 f3(t)f_3(t) 取第二项的微分方程组。
具体地，我们的算法流程应该为：
  1. 求解如下常系数一阶线性微分方程组：
{f0′(t)=f0(t)ln⁡p0−f3(t)ln⁡p0f1′(t)=1+f1(t)ln⁡p1−(q0f0(t)+(1−q0)f2(t))ln⁡p1f2′(t)=1+f2(t)ln⁡p2−f1(t)ln⁡p2f3(t)=f1(t)f4(t)=f2(t)f0(0)=f1(0)=f2(0)=f3(0)=f4(0)=0\left\{\begin{matrix}f_0'(t)=f_0(t)\ln p_0-f_3(t)\ln p_0 \\f_1'(t)=1+f_1(t)\ln p_1-(q_0f_0(t)+(1-q_0)f_2(t))\ln p_1 \\f_2'(t)=1+f_2(t)\ln p_2-f_1(t)\ln p_2 \\f_3(t)=f_1(t) \\f_4(t)=f_2(t) \\f_0(0)=f_1(0)=f_2(0)=f_3(0)=f_4(0)=0\end{matrix}\right.
  1. 二分找到满足 −q0(f1(t)−f0(t))ln⁡p1≥1-q_0(f_1(t)-f_0(t))\ln p_1\ge 1 的最小的 tt 为 vv，则上述解的定义域为 [0,l][0,l]。
  2. 求解如下常系数一阶线性微分方程组，其中 t≥lt\ge l，初值为 t=lt=l 时与第一步中的解数值相同：
{f0′(t)=f0(t)ln⁡p0−f3(t)ln⁡p0f1′(t)=1+f1(t)ln⁡p1−(q0f0(t)+(1−q0)f2(t))ln⁡p1f2′(t)=1+f2(t)ln⁡p2−f1(t)ln⁡p2f3′(t)=f3(t)ln⁡p1−(q0f3(t)+(1−q0)f4(t))ln⁡p1f4(t)=f2(t)\left\{\begin{matrix}f_0'(t)=f_0(t)\ln p_0-f_3(t)\ln p_0 \\f_1'(t)=1+f_1(t)\ln p_1-(q_0f_0(t)+(1-q_0)f_2(t))\ln p_1 \\f_2'(t)=1+f_2(t)\ln p_2-f_1(t)\ln p_2 \\f_3'(t)=f_3(t)\ln p_1-(q_0 f_3(t)+(1-q_0)f_4(t))\ln p_1 \\f_4(t)=f_2(t) \end{matrix}\right.
  1. 答案为 f0(k)f_0(k)。
最后的问题是，如何求解一阶线性微分方程组。
第一种方法是求精确解。我们有：
【常系数一阶线性微分方程组解的结构】 对于常系数一阶线性微分方程组 dX⃗dt=AX⃗+F\dfrac{\mathrm{d}\vec X}{\mathrm{d}t}=A\vec X+F，其中对向量求导定义为对向量每一维求导，其通解形式为 X⃗=eAtX⃗(0)+φ⃗\vec X=e^{At}\vec X(0)+\vec \varphi，其中向量函数定义为对每一维代入自变量，矩阵指数定义为泰勒级数，φ⃗\vec\varphi 为特解。
至于特解，官方题解说可以观察到 AA 每一行的和都是 00，所以 00 是 AA 的特征值，故特解可以被看成线性函数。
但是这种方法麻烦完了。既然这题概率不取模，我们完全可以用基于精度的方法。
考虑四阶 Runge–Kutta 法。设 y=y(x)y=y(x)，下面简述对于微分方程 y′=f(x,y)y'=f(x,y) 如何求，微分方程组改成向量就行。
我们取步长 hh，设 xi=ih+x0,yi=y(xi)x_i=ih+x_0,y_i=y(x_i)。
使用拉格朗日中值定理，存在 t∈[xi,xi+1]t\in[x_i,x_{i+1}]，使得 yi+1=yi+hy′(t)y_{i+1}=y_i+hy'(t)。
我们使用 [xi,xi+1][x_i,x_{i+1}] 区间上的平均斜率来近似 y′(t)y'(t)。
设 k1=f(xi,yi),k2=f(xi+h2,yi+h2k1),k3=f(xi+h2,yi+h2k2),k4=f(xi+h,yi+hk3)k_1=f(x_i,y_i),k_2=f(x_i+\frac h 2,y_i+\frac h 2k_1),k_3=f(x_i+\frac h 2,y_i+\frac h 2k_2),k_4=f(x_i+h,y_i+hk_3)。把它们加权平均得到 k∗=16(k1+2k2+2k3+k4)k^*=\frac 1 6(k_1+2k_2+2k_3+k_4)，直接估算 y′(t)=k∗y'(t)=k^*。
取步长 h=kCh=\dfrac k C，时间复杂度为 O(TC)O(TC)。我取的是 C=5×104C=5\times 10^4。

代码

CPP
#include<bits/stdc++.h>
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
using namespace std;
const int C=50000;
double k,p0,p1,p2,q0,h;
inline array<double,4> calc1(double x,const array<double,4>&y){
	array<double,4>res;
	res[0]=y[0]*p0-y[3]*p0;
	res[1]=1+y[1]*p1-(q0*y[0]+(1-q0)*y[2])*p1;
	res[2]=1+y[2]*p2-y[1]*p2;
	res[3]=res[1];
	return res;
}
inline array<double,4> calc2(double x,const array<double,4>&y){
	array<double,4>res;
	res[0]=y[0]*p0-y[3]*p0;
	res[1]=1+y[1]*p1-(q0*y[0]+(1-q0)*y[2])*p1;
	res[2]=1+y[2]*p2-y[1]*p2;
	res[3]=y[3]*p1-(q0*y[3]+(1-q0)*y[2])*p1;
	return res;
}
array<double,4> operator+(const array<double,4>&a,const double&b){return {a[0]+b,a[1]+b,a[2]+b,a[3]+b};}
array<double,4> operator+(const array<double,4>&a,const array<double,4>&b){return {a[0]+b[0],a[1]+b[1],a[2]+b[2],a[3]+b[3]};}
array<double,4> operator*(const double&b,const array<double,4>&a){return {a[0]*b,a[1]*b,a[2]*b,a[3]*b};}
double x[C*2];
array<double,4>y[C*2];
inline void solve1(){
	y[0]={0,0,0,0};
	array<double,4>k1,k2,k3,k4;
	F(i,1,C){
		k1=calc1(x[i-1],y[i-1]);
		k2=calc1(x[i-1]+h/2,y[i-1]+h/2*k1);
		k3=calc1(x[i-1]+h/2,y[i-1]+h/2*k2);
		k4=calc1(x[i-1]+h,y[i-1]+h*k3);
		y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);
	}
}
inline void solve2(int l){
	array<double,4>k1,k2,k3,k4;
	F(i,l,C){
		k1=calc2(x[i-1],y[i-1]);
		k2=calc2(x[i-1]+h/2,y[i-1]+h/2*k1);
		k3=calc2(x[i-1]+h/2,y[i-1]+h/2*k2);
		k4=calc2(x[i-1]+h,y[i-1]+h*k3);
		y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);
	}
}
signed main(){
	ios::sync_with_stdio(0);cin.tie(0),cout.tie(0);
	int T;
	cout<<fixed<<setprecision(12);
	for(cin>>T;T;--T){
		cin>>k>>p0>>p1>>p2>>q0;
		p0=log(p0),p1=log(p1),p2=log(p2);
		h=k/C;
		F(i,0,C) x[i]=i*h;
		solve1();
		int l=0,r=C+1;
		while(l<r){
			int mid=(l+r)>>1;
			if(-q0*(y[mid][1]-y[mid][0])*p1>=1) r=mid;
			else l=mid+1;
		}
		if(l!=C+1) solve2(l);
		cout<<y[C][0]<<"\n";
	}
	return 0;
}

评论

共 0 条评论，欢迎与作者交流。

正在加载评论...