[BZOJ]4872: [Shoi2017]分手是祝愿 期望DP

Description

Zeit und Raum trennen dich und mich.
时空将你我分开。B 君在玩一个游戏，这个游戏由 n 个灯和 n 个开关组成，给定这 n 个灯的初始状态，下标为从 1 到 n 的正整数。每个灯有两个状态亮和灭，我们用 1 来表示这个灯是亮的，用 0 表示这个灯是灭的，游戏的目标是使所有灯都灭掉。但是当操作第 i 个开关时，所有编号为 i 的约数（包括 1 和 i）的灯的状态都会被
改变，即从亮变成灭，或者是从灭变成亮。B 君发现这个游戏很难，于是想到了这样的一个策略，每次等概率随机操作一个开关，直到所有灯都灭掉。这个策略需要的操作次数很多， B 君想到这样的一个优化。如果当前局面，可以通过操作小于等于 k 个开关使所有灯都灭掉，那么他将不再随机，直接选择操作次数最小的操作方法（这个策略显然小于等于 k 步）操作这些开关。B 君想知道按照这个策略（也就是先随机操作，最后小于等于 k 步，使用操作次数最小的操作方法）的操作次数的期望。这个期望可能很大，但是 B 君发现这个期望乘以 n 的阶乘一定是整数，所以他只需要知道这个整数对 100003 取模之后的结果。

Solution

最优策略显然是从大到小看每个灯，需要操作就操作。当$n=k$时就直接用这个来做，期望得分$50$，在LOJ上得到了$80$的好成绩。
然后可以发现，每次操作本质上是增加一次操作次数或减少一次操作次数，问题可以转化为，有$n+1$个格子，编号为$0$~$n$，设一开始的状态需要的最优步数为$cnt$，也就是站在编号为$cnt$的格子，在编号为$x$的格子有$\frac{x}{n}$的概率到$x-1$，有$\frac{n-x}{n}$的概率到$x+1$，问到达$k$的期望步数。
状态之间互相影响，显然可以高斯消元，但复杂度不是$O(n^3)$，因为在$n$下一步一定到$n-1$，所以可以直接递推。但据大神blog，某些系数可能在模意义下为$0$，这个做法只能获得$95$分。（但是ccz大爷好像AC了）
考虑设计另外一种状态，使得状态间不会互相影响。$f_i$表示$i$走到$i-1$的期望步数，每次要么走到$i-1$，要么回到$i+1$，所以有转移$f_i=\frac{i}{n}+\frac{n-i}{n}(f_i+f_{i-1}+1)$，直接转移就行了，最后答案是一段$f$的和。

Code

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define pa pair<int,int>
const int Maxn=100010;
const int inf=2147483647;
const int mod=100003;
int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return x*f;
}
int n,k,a[Maxn],inv[Maxn],f[Maxn];
vector<int>h[Maxn];
int main()
{
	n=read(),k=read();
	inv[0]=inv[1]=1;
	for(int i=2;i<=n;i++)inv[i]=(LL)(mod-mod/i)*inv[mod%i]%mod;
	for(int i=1;i<=n;i++)
	for(int j=i;j<=n;j+=i)
	h[j].push_back(i);
	for(int i=1;i<=n;i++)a[i]=read();
	int ans=0;
	for(int i=n;i;i--)
	if(a[i])
	{
		ans++;
		for(int j=0;j<h[i].size();j++)a[h[i][j]]^=1;
	}
	if(ans<=k)
	{
		ans%=mod;
		for(int i=2;i<=n;i++)ans=(LL)ans*i%mod;
		return printf("%d",ans),0;
	}
	int Ans=k;f[n]=1;
	for(int i=n-1;i>k;i--)
	f[i]=(LL)n*inv[i]%mod*(1+(LL)(n-i)*inv[n]%mod*f[i+1]%mod)%mod;
	for(int i=ans;i>k;i--)Ans=(Ans+f[i])%mod;
	for(int i=2;i<=n;i++)Ans=(LL)Ans*i%mod;
	printf("%d",Ans);
}