[六省联考2017]分手是祝愿

题面在这里

题意

有n盏灯,当前状态为亮或者不亮,当改变第x盏灯的开关状态时(由亮变暗，由暗变亮)，为x约数编号的灯也会改变开关状态
B先生先随机操作，当存在一种能使用小于k的操作次数让灯全部熄灭的方案的时候，B先生会使用操作次数最小的操作方法
询问期望步数。n<=100000

sol

直接看这道题显然非常头疼

先考虑如何找到熄灭所有灯的最小步数
考虑单次操作:对于当前操作的灯,编号比它大的灯显然不会改变任何状态
并且在最优方案中，对同一个灯进行两次操作也显然是不合理的
那么我们就有了这样一种办法:从n到1考察所有灯，当这盏灯的当前状态为亮时对其进行一次操作,如果不亮则不操作;
求出x的约数序列，只需要在埃氏筛法预处理的基础上加n个vector就可以了(具体参见代码)
于是我们得到了最小操作次数\(tot\)

继续思考,我们可以了解到对于每一盏灯的操作都是独立的，因为任何一盏灯的控制范围都不能被任意两盏灯的控制范围所代替
于是B先生要想熄灭所有的灯,必须在最优方案包含的灯上操作奇数次
那么可以转换模型:
现在我们有\(tot\)个1,\((n-tot)\)个0，B先生相当于有\(\frac{tot}{n}\)的概率把一个1变成0，
有\(\frac{n-tot}{n}\)的概率把一个0变成1，B先生的目的是要把所有1变成0。

于是答案就只和\(tot\)，\(n\)和\(k\)有关，当前所有灯状态可以囊括为\(f[x]\)，即最少需要x步将所有灯熄灭(有x个1)，那么对于\(x<=k\)显然有\(f[x]=x\)
对于\(k<x<n\)，我们有这样一个递推式:$$f[x]=f[x-1]\times\frac{x}{n}+f[x+1]\times\frac{n-x}{n}+1$$
最后\(f[n]=f[n-1]+1\)
把\(f[n]=f[n-1]+1\)代入上面的f[x]可以得到下面这个奇怪的式子:

\[f[x]=f[x-1]+\frac{n}{x}+\frac{n-x}{x}\times(\frac{n}{x+1}+\frac{n-(x+1)}{x+1}\times(...\times(\frac{n}{n-2}+\frac{2}{n-2}\times(\frac{n}{n-1}+\frac{1}{n-1}))...)) \]

可以看到这个式子是递归的，所以本蒟蒻使用了一个递归函数求解。。。
注意逆元和取模

代码

//update 10.8:感谢@ve_2021的提醒，稍微简化了一下代码
#include<bits/stdc++.h>
#define pb push_back
#define RG register
#define il inline
using namespace std;
const int mod=1e5+3;
const int N=100010;
typedef unsigned long long ull;
typedef vector<int>VI;
typedef long long ll;
typedef double dd;
il ll read(){
  RG ll data=0,w=1;RG char ch=getchar();
  while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
  if(ch=='-')w=-1,ch=getchar();
  while(ch<='9'&&ch>='0')data=data*10+ch-48,ch=getchar();
  return data*w;
}

il ll fac(ll x){if(!x)return 1;return 1ll*fac(x-1)*x%mod;}
il ll poww(ll a,ll b){
  RG ll ret=1;for(;b;b>>=1,a=a*a%mod)if(b&1)ret=ret*a%mod;return ret;}
int n,k,tot,rev[N],f[N];
bool light[N];

VI cz[N];
il void sieve(){
  for(RG int i=1;i<=n;i++)cz[i].pb(1);
  for(RG int i=2;i<=n;i++){
    for(RG int s=1;1ll*s*i<=n;s++)
      cz[s*i].push_back(i);
  }
  for(RG int i=1;i<=n;i++)rev[i]=poww(i,mod-2);
}

il ll search(int x,ll sum){
  //核心部分,递归计算f[x]
  if(x<=k)return f[x]=x;
  sum=(1ll*n*rev[x]%mod+sum*(n-x)%mod*rev[x]%mod)%mod;
  return f[x]=(search(x-1,sum)+sum)%mod;
}

int main()
{
  n=read();k=read();sieve();
  for(RG int i=1;i<=n;i++)light[i]=read();
  for(RG int i=n,sz;i;i--)
    if(light[i]){
      //有亮着的灯就关
      tot++;sz=cz[i].size();
      for(RG int j=0;j<sz;j++)
	light[cz[i][j]]^=1;
    }
    
  search(n,1);
  printf("%lld\n",tot>k?f[tot]*fac(n)%mod:tot*fac(n)%mod);
  //如果总步数<=k则直接一步到位
  return 0;
}