[SHOI2017]分手是祝愿

题目：
https://ac.nowcoder.com/acm/problem/20437

有$n$盏灯，$1$表示开的，$0$表示关的。每次操作随机选一个开关$i(1\le i\le n)$，对于所有的$j(j|i)$，第$j$盏灯的状态反转。如果当前在最优的情况下只需要$x(x\le k)$个开关就能使所有的灯关掉，则按照最优的情况进行操作。求关掉所有灯的期望次数。

思路：

• 首先考虑如何对于一个状态求最小步数。
  由于编号小的开关不会影响到后面的灯，所以应该按编号从大到小决策。对于当前枚举到的灯$i$，如果它是开着的，那么必须使用开关$i$。所以可以很快求出对于初始状态的步数。

1. 最优情况每个开关最多操作一次，所以最优情况最多为$n$。
2. 最优操作的顺序没有关系。
3. 最优操作唯一(不考虑顺序)

• $dp[i]$表示在最优情况下还需要$i$次操作才能把所有灯关掉的状态到完成态的期望操作次数。(因为关掉一盏灯会影响还需要的最优步数)那么接下来有$\frac{i}{n}$的概率到达最优情况下还要$i-1$次操作的状态，$\frac{n-i}{n}$的概率到达$i+1$的状态。(如果使用的就是最优策略的$i$个开关中的一个，那么最优次数减$1$，否则加$1$)
  $$dp[i]=\{\begin{array}{l}\frac{i}{n}\cdot dp[i-1]+\frac{n-i}{n}\cdot dp[i+1]+1(k<i\le n)①\\ i(i\le k)②\end{array}$$
  设初始最优要$num$次，则$dp[0]=0$，$dp[num]$为答案。

1. 解方程组，时间复杂度太高
2. 看成数列，构造
   $$\begin{array}{rl}dp[i]& =\frac{i}{n}\cdot dp[i-1]+\frac{n-i}{n}\cdot dp[i+1]+1\\ \frac{i}{n}\cdot (dp[i]-dp[i-1])& =\frac{n-i}{n}\cdot (dp[i+1]-dp[i])+1令dp[i]-dp[i-1]=f[i]\\ f[i]& =\frac{n-i}{i}\cdot f[i+1]+\frac{n}{i}f[n]=1k<i\le n\end{array}$$
   从后往前递推$f[i]$。如果$num>k,dp[num]=dp[k]+{\sum }_{i=k+1}^{num}f[i]$

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int N=100005,mo=100003;

typedef long long LL;

int n,Fac[N],C[N],m,k,bz[N],Inv[N],ans;

char c;

int read()
{
    int x=0,sig=1;
    for (c=getchar();c<'0' || c>'9';c=getchar()) if (c=='-') sig=-1;
    for (;c>='0' && c<='9';c=getchar()) x=x*10+c-48;
    return x*sig;
}

int main()
{
    n=read(); k=read();
    for (int i=1;i<=n;i++) bz[i]=read();
    for (int i=n;i;i--) if (bz[i])
    {
        m++;
        for (int j=1;j*j<=i;j++) if (i%j==0)
        {
            bz[j]^=1;
            if (i/j>j) bz[i/j]^=1;
        }
    }
    Fac[0]=Fac[1]=Inv[1]=1;
    for (int i=2;i<=n;i++)
    {
        Fac[i]=(LL)Fac[i-1]*i%mo;
        Inv[i]=(LL)Inv[mo%i]*(mo-mo/i)%mo;
    }
    if (k>=m)
    {
        printf("%d\n",(LL)Fac[n]*m%mo);
        return 0;
    }
    C[0]=1;
    for (int i=0;i<n;i++) C[i+1]=(LL)C[i]*(n-i)%mo*Inv[i+1]%mo;
    for (int i=1;i<=n;i++) C[i]=(C[i]+C[i-1])%mo;
    for (int i=k+1;i<=m;i++) ans=(ans+(LL)Fac[n-i]*Fac[i-1]%mo*(C[n]-C[i-1]))%mo;
    if (ans<0) ans+=mo;
    ans=(LL)ans*n%mo;
    ans=(ans+(LL)Fac[n]*k)%mo;
    printf("%d\n",ans);
    return 0;
}