【SHOI2017】【期望dp】分手是祝愿

【描述】
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 取模之后的结果。

【输入】
第一行两个整数 n; k。

接下来一行 n 个整数，每个整数是 0 或者 1，其中第 i 个整数表示第 i 个灯的初始情况。

【输出】
输出一行，为操作次数的期望乘以 n 的阶乘对 100003 取模之后的结果。

• 对于 0% 的测试点，和样例一模一样；

• 对于另外 30% 的测试点， n ≤ 10；

• 对于另外 20% 的测试点， n ≤ 100；

• 对于另外 30% 的测试点， n ≤ 1000；

• 对于 100% 的测试点， 1 ≤ n ≤ 100000; 0 ≤ k ≤ n；

• 对于以上每部分测试点，均有一半的数据满足 k = n。

【思路】

除了预处理约数看起来无从下手。不过不管写不写代码，我们已经拿到了第一档部分分。
我们先来考虑考虑最优策略是怎么样的。简单分析得出下列性质：
1.按键是独特的，不能通过一些按键的组合达到和一个按键等价的效果。
2.对于灯的若干次操作等价于两种操作：操作一次，不操作。
这是两个很有趣的性质，这意味着：当每个灯都关上时，每个灯的等价意义下的操作是确定的。那么考虑贪心，从大到小依次考虑，如果灯亮则按灯，这一定是最优策略。我们这样就可以统计出我们需要对多少盏灯进行操作。由于第一条性质，我们知道在任何策略中，等价意义下这些灯是必须操作的且只能操作这些灯。我们考虑状态$f[i]$表示还剩i盏灯需要操作的状态下，转移到还剩i-1盏灯需要操作的期望操作次数。我们有$i/n$的概率进行正确的操作，此时操作次数为1，我们有$(n-i)/n$的概率错误操作，此时，根据性质1，我们必须在以后把这次错误操作修正回来，所以操作次数为f[i+1]+f[i]+1。移项化简得到f[i]的递推式即可进行递推。设need表示我们需要操作的灯的数量，答案显然为${\sum }_{i=k+1}^{need}f[i]+k$。
代码：

#include<bits/stdc++.h>
#define re register
using namespace std;
const int N=1e5+5;
const int mod=1e5+3;
inline int red(){
    int data=0;bool w=0; char ch=0;
    ch=getchar();
    while(ch!='-' && (ch<'0' || ch>'9')) ch=getchar();
    if(ch=='-') w=1,ch=getchar();
    while(ch>='0' && ch<='9') data=(data<<3)+(data<<1)+ch-'0',ch=getchar();
    return w?-data:data;
}
int n,m,b,c,num=0;
bool a[N];
vector<int>g[N];
double f[N];
int fac[N]={1,1};
int inv[N]={0,1};
int ans=0;
inline int mul(const int&a,const int &b){return 1ll*a*b%mod;}
void pre(){
	for(int re i=2;i<=n;i++){
		fac[i]=mul(fac[i-1],i);
		inv[i]=mul(inv[mod%i],mod-mod/i);
	}
}
int main()
{
	n=red();m=red();pre();
	for(int re i=1;i<=n;i++)a[i]=red();
	for(int re i=1;i<=n;i++)
		for(int re j=i;j<=n;j+=i)
			g[j].push_back(i);
	for(int re i=n;i;--i){
		if(a[i]){++num;
			for(int re j=g[i].size()-1;~j;--j)
				a[g[i][j]]^=1;
		}
	}
	if(num<m)return cout<<mul(num,fac[n])<<"\n",0;
	f[n]=1;
	for(int re i=n-1;i;--i)f[i]=(1+mul(mul(n-i,f[i+1]+1),inv[i]))%mod;
	for(int re i=num;i>m;--i)(ans+=f[i])%=mod;
	cout<<mul(ans+m,fac[n])<<"\n";
}