[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]=iN×f[i−1]+N−1N×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+−(N−i)∗R−N+NQi，分母就是乘上逆元。最后把期望乘上阶乘即为答案

至于为什么得到的每个f[i]都可以表示成 x+常数 的形式 ，me写了一页草稿纸= =，有兴趣的可以自己推导。像这种一项与相邻两项有关系的情况，一般方法是两项差分，也就是设b[i]=f[i]−f[i−1]，有 b[i]=in+(1−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 ; 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
%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
;
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() ;
}