BZOJ4872 [SHOI2017]分手是祝愿

题意

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

题解

先考虑对于某一个状态,如何求出最小的操作次数使得灯全部熄灭。我们可以想到这样一个贪心策略:从编号大的往编号小的走,如果没有灭灯就把它灭了。容易发现这样一定是最小步数。然后来考虑怎么DP。设$dp[i]$表示最少从仍需要$i$步才能将灯全部熄灭到还有$i-1$步的步数期望。然后就有这么个转移方程:
$$dp[i]=\frac{i}{n}+(1-\frac{i}{n}) \times (dp[i+1]+dp[i]+1)$$
式子不难理解,前面就是这一步改对了(即没有让期望步数增加),后面就是改错的期望,就需要再花$dp[i+1]+dp[i]$才能改回来。然后就可以DP了。
时间复杂度的瓶颈在于第一步求步数,复杂度为$O(n\ln n)$。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int INF=1e9;
const long double eps=1e-9;
const int maxn=1e5+10;
const int mod=1e5+3;
vector <int> G[maxn];
int a[maxn],dp[maxn],inv[maxn];
inline int read(){
int x=0,flag=1;
char ch=getchar();
while(!isdigit(ch) && ch!='-')ch=getchar();
if(ch=='-')flag=-1,ch=getchar();
while(isdigit(ch))x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return x*flag;
}
inline int Mod(int val){
return val>=mod ? val-mod : val;
}
int main(){
int i,j,k,m,n;
#ifndef ONLINE_JUDGE
freopen("BZOJ4872.in","r",stdin);
freopen("BZOJ4872.out","w",stdout);
#endif
n=read();k=read();
inv[1]=1;
for(i=2;i<=n;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
for(i=1;i<=n;i++)a[i]=read();
for(i=1;i<=n;i++)
for(j=i;j<=n;j+=i)
G[j].push_back(i);
int cnt=0;
for(i=n;i>=1;i--){
if(!a[i])continue;
cnt++;
for(j=0;j<(int)G[i].size();j++)
a[G[i][j]]^=1;
}
int ans=0;
if(cnt<=k)ans=cnt;
else{
dp[n]=1;
for(i=n-1;i>k;i--)dp[i]=1ll*(n-i)*(dp[i+1]+1)%mod*inv[i]%mod+1;
for(i=cnt;i>k;i--)ans=Mod(ans+dp[i]);
ans=Mod(ans+k);
}
for(i=1;i<=n;i++)ans=1ll*ans*i%mod;
printf("%d\n",ans);
return 0;
}