[BZOJ4872]-[Shoi2017]分手是祝愿-期望DP+画柿子

说在前面

啊啊啊啊啊啊= =
模拟考试考了这套题，考场上推导了一个多小时终于搞出来了，然而没注意到模数是质数，以为求不出逆元…n=k的情况me又没有乘上阶乘，最后只拿到了5分
心中有100…0000万句mmp= =

---

题目

BZOJ4872传送门

题面

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$ 个灯的初始情况。
$1 ≤ N ≤ 100000\ , 0 ≤ k ≤ N$

---

解法

考虑在某个局面，要把所有灯关掉的最少步数是什么

如果先关编号较小的灯，对于编号比他大的灯是没有影响的，然而在这样的策略下，关掉编号较大的灯有可能又把编号较小的灯打开了，从而浪费很多操作步数。
所以正确的策略是应该先关编号较大的灯，这样的话首先每个灯最多只会被关一次（显然）。其次就是，假设比当前灯编号大的灯都被关掉了，如果现在灯是开的，那么无论如何这个开关都必须要按下1次（奇数次），才能把这个灯关掉（因为关掉它的倍数的灯的时候总会影响到它）
同时可以发现，那些需要被摁下的开关都得去摁下1次（奇数次），所以与摁下的顺序无关（可以用异或关系理解）

具体的做法

因此可以先从大到小用 根号枚举约数 的方式，计算出初始局面关掉所有灯的最小步数，如果这个数小于等于$K$就直接输出了（这样的数据有$80$分，简直一点区分度都没有= =）
现在设出$f[i]$ 表示在最优策略下还需要摁下$i$个开关才能关掉所有灯，明显有$f[i] = \frac{i}{N}\times f[i-1]+\frac{N-1}{N}\times f[i+1]$，并且可以轻易的发现边界情况$f[K]=K$以及$f[N] = f[N-1]+1$。

根据最后一个条件带回到f[N-1]的式子里，发现 f[N-2]也可以用f[N-1]表示出来。一直倒推回去最后就可以求出 f[K]关于f[N-1]的表达式，然后同样可以得到 f[当前最优步数]关于f[N-1]的表达式，然后就可以求出答案了。
令f[N-1]为x，则有f[i]=x+Q，f[i+1]=x+R，带入f[i]的式子，那么有$f[i-1]=x+\frac{-(N-i)*R-N+NQ}{i}$，分母就是乘上逆元。最后把期望乘上阶乘即为答案

至于为什么得到的每个f[i]都可以表示成 x+常数 的形式 ，me写了一页草稿纸= =，有兴趣的可以自己推导。像这种一项与相邻两项有关系的情况，一般方法是两项差分，也就是设$b[i]=f[i]-f[i-1]$，有 $b[i]=\frac in+(1-\frac in)·(1+b[i+1]+b[i])$，边界$b[i]=1$，倒推即可，答案就是b[K+1]到b[当前最优步数]之和

---

下面是自带大常数的代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

int N , K , a[100005] , cnt ;
int fx[100005] , fac[100005] , inv[100005] , mmod = 100003 ;

void preWork(){
    for( int i = N , j ; i >= 1 ; i -- )
        if( a[i] ){
            for( j = 1 , cnt ++ ; j * j < i ; j ++ )
                if( i % j == 0 ) a[j] ^= 1 , a[i/j] ^= 1 ;
            if( j * j == i ) a[j] ^= 1 ;
        }
    fx[N] = N ; fac[0] = fac[1] = 1 ;
    for( int i = N - 1 ; i >= 1 ; i -- )
        fx[i] = 1LL * fx[i+1] * i % mmod ; 
    for( int i = 2 ; i <= N ; i ++ )
        fac[i] = 1LL * fac[i-1] * i %mmod ;

    inv[1] = 1 ;
    for( int i = 2 ; i <= N ; i ++ ){
        inv[i] = -1LL * ( mmod / i ) * inv[ mmod%i ] %mmod ;
        inv[i] = ( inv[i] %mmod + mmod )%mmod ;
    }
}

void solve(){
    if( cnt <= K ){
        printf( "%d" , 1LL * cnt * fac[N] %mmod ) ;
        return ;
    }
    long long Rk = 0 , Rcnt = 0 , ans , Q = 0 , R = 1 ;
    for( int i = N - 1 ; i > K ; i -- ){
        Rk = ( -1LL * N * R %mmod + i * R %mmod - N ) * inv[i] ;
        Rk = ( Rk %mmod + N * Q %mmod * inv[i] %mmod ) ;
        Rk = ( Rk %mmod + mmod )%mmod ;
        R = Q ; Q = Rk ;
    }

    if( cnt == N ) Rcnt = 1 ;
    Q = 0 , R = 1 ;
    for( int i = N - 1 ; i > cnt ; i -- ){
        Rcnt = ( -1LL * N * R %mmod + i * R %mmod - N ) * inv[i] ;
        Rcnt = ( Rcnt %mmod + N * Q %mmod * inv[i] %mmod ) ;
        Rcnt = ( Rcnt %mmod + mmod )%mmod ;
        R = Q ; Q = Rcnt ;
    }
    ans = ( K - Rk + Rcnt ) * fac[N] ;
    printf( "%lld" , ( ans %mmod + mmod )%mmod ) ;
}

int main(){
//  freopen( "trennen.in" , "r" , stdin ) ;
//  freopen( "trennen.out", "w" , stdout) ;
    scanf( "%d%d" , &N , &K ) ;
    for( int i = 1 ; i <= N ; i ++ )
        scanf( "%d" , &a[i] ) ;
    preWork() ;
    solve() ;
}