bzoj 4872: [Shoi2017]分手是祝愿

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

Input

第一行两个整数 n, k。

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

1 ≤ n ≤ 100000, 0 ≤ k ≤ n；

Output

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

预处理出每个开关应被操作奇数次还是偶数次，设当前有c（c>k）个开关的操作次数奇偶与目标不同，则有c/n概率转移到c-1，1-c/n概率转移至c+1

这相当于在一条链上移动，初始在c，求到达k的期望步数，走到k立即结束

设$p(a,b)$表示走了a步，走到b的概率

设$P(x)$,$E(x)$分别为x状态期望被到达的次数和每次到达的步数之和的期望值，即

$P(x)=sum_{i=0}^{+infty}p(i,x)$

$E(x)=sum_{i=0}^{+infty}i cdot p(i,x)$

则

$P(x)=[x>k+1] cdot frac{x-1}{n} cdot P(x-1)+[x<n] cdot frac{n-x-1}{n} cdot P(x+1)+[x==c]$
$E(x)=[x>k+1] cdot frac{x-1}{n} cdot (E(x-1)+P(x-1))+[x<n] cdot frac{n-x-1}{n} cdot (E(x+1)+P(x+1))$

得到$O(n)$条包含至多3个变量的方程，按正常顺序高斯消元恰好可以直接$O(n)$得到每个变量的值

#include<cstdio>
typedef long long i64;
const int P=100003;
int n,k,d[P],fac[P],fiv[P],iv[P];
int inv(int a){
    int v=1;
    for(int n=P-2;n;n>>=1,a=1ll*a*a%P)if(n&1)v=1ll*v*a%P;
    return v;
}
int p[P];
long long xs[P][4];
int main(){
    scanf("%d%d",&n,&k);
    int c=0;
    for(int i=1;i<=n;++i)scanf("%d",d+i);
    for(int i=fac[0]=1;i<=n;++i)fac[i]=i64(i)*fac[i-1]%P;
    fiv[n]=inv(fac[n]);
    for(int i=n;i;--i)fiv[i-1]=1ll*i*fiv[i]%P;
    for(int i=1;i<=n;++i)iv[i]=1ll*fiv[i]*fac[i-1]%P;
    for(int i=n;i;--i)if(d[i]){
        for(int j=1;j*j<=i;++j)if(i%j==0){
            d[j]^=1;
            if(j*j!=i)d[i/j]^=1;
        }
        ++c;
    }
    if(c<=k){
        printf("%lld",1ll*c*fac[n]%P);
        return 0;
    }
    for(int i=k;i<=n;++i){
        if(i>k+1)xs[i][0]=-(1-1ll*(i-1)*iv[n])%P;
        xs[i][1]=1;
        if(i<n)xs[i][2]=-1ll*(i+1)*iv[n]%P;
        xs[i][3]=0;
    }
    xs[c][3]=1;
    
    for(int i=n;i>k;--i){
        int c=xs[i-1][2]*inv(xs[i][1])%P;
        xs[i-1][3]=(xs[i-1][3]-xs[i][3]*c)%P;
        xs[i-1][2]=0;
        xs[i-1][1]=(xs[i-1][1]-xs[i][0]*c)%P;
    }
    p[k]=1;
    for(int i=k+1;i<=n;++i){
        p[i]=(xs[i][3]-xs[i][0]*p[i-1])%P*inv(xs[i][1])%P;
    }
    for(int i=k;i<=n;++i){
        xs[i][0]=0;
        xs[i][1]=1;
        xs[i][2]=0;
        xs[i][3]=0;
        if(i>k+1){
            xs[i][0]=-(1-(i-1)*1ll*iv[n])%P;
            xs[i][3]=(xs[i][3]-xs[i][0]*p[i-1])%P;
        }
        if(i<n){
            xs[i][2]=-(i+1)*1ll*iv[n]%P;
            xs[i][3]=(xs[i][3]-xs[i][2]*p[i+1])%P;
        }
    }
    for(int i=n;i>k;--i){
        int c=xs[i-1][2]*inv(xs[i][1])%P;
        xs[i-1][3]=(xs[i-1][3]-xs[i][3]*c)%P;
        xs[i-1][2]=0;
        xs[i-1][1]=(xs[i-1][1]-xs[i][0]*c)%P;
    }
    int res=xs[k][3]*inv(xs[k][1])%P;
    int ans=1ll*(res+k)%P*fac[n]%P;
    printf("%d
",(ans+P)%P);
    return 0;
}