【SHOI&SXOI2017】bzoj4872 分手是祝愿

不难从大到小贪心求出对于当前状态，最优解需要几步。
记$f_i$为如果最优解需要$i$步，期望需要的步数。可以发现

$$f_i= \begin{cases} i & i\le k \\ \frac in f_{i-1}+\frac{n-i}nf_{i+1}+1 & i\gt k \\ \end{cases}$$
这样是不能直接求的，不妨记$g_i=f_i-f_{i-1}$，那么

$$g_i= \begin{cases} 1 & i\le k \\ \frac{n-i}i g_{i+1} +\frac ni & i\gt k\\ \end{cases}$$
又因为 $g_n=1$，就可以递推求出所有 $g$，进而求出$f$。

#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
#define LL long long
const int maxn=100010,p=100003;
vector<int> v[maxn];
vector<int>::iterator it;
int a[maxn],f[maxn],g[maxn],inv[maxn],prm[maxn],
n,m,k,tot;
int pow(int base,int k)
{
    int ret=1;
    for (;k;k>>=1,base=(LL)base*base%p)
        if (k&1) ret=(LL)ret*base%p;
    return ret;
}
int main()
{
    int ans;
    scanf("%d%d",&n,&k);
    for (int i=1;i<=n;i++) scanf("%d",&a[i]);
    for (int i=1;i<=n;i++)
    {
        if (!inv[i])
        {
            inv[i]=pow(i,p-2);
            prm[++tot]=i;
        }
        for (int j=1;j<=tot&&(LL)i*prm[j]<=n;j++)
        {
            inv[i*prm[j]]=inv[i]*inv[prm[j]];
            if (i%prm[j]==0) break;
        }
    }
    for (int i=1;i<=n;i++)
        for (int j=i;j<=n;j+=i)
            v[j].push_back(i);
    for (int i=n;i;i--)
        if (a[i])
        {
            m++;
            for (it=v[i].begin();it!=v[i].end();++it) a[*it]^=1;
        }
    g[n]=1;
    for (int i=1;i<=k;i++) g[i]=1;
    for (int i=n-1;i>k;i--) g[i]=((LL)g[i+1]*(n-i)%p*inv[i]+(LL)n*inv[i])%p;
    for (int i=1;i<=k;i++) f[i]=i;
    for (int i=k+1;i<=m;i++) f[i]=f[i-1]+g[i];
    ans=f[m];
    for (int i=1;i<=n;i++) ans=(LL)ans*i%p;
    printf("%d\n",ans);
}