【每日一题】Steps to One (容斥+错位相减)cf1139D

题目链接
前置推论:
1.长度为1的序列只可能是一个1。
2.假设当前的序列长度为 k + 1 k+1 k+1,且前 k k k个数的 g c d = a gcd=a gcd=a不为1的话,那么第 k + 1 k+1 k+1位的数必然与a互质,所以第 k + 1 k+1 k+1个数选择的方案为[1,m]中与a互质的数的个数,记 g e t ( a ) get(a) get(a) [ 1 , m ] [1,m] [1,m]中与a互质的数的个数, g e t ( x ) get(x) get(x)可以通过容斥计算。
那么我们需要统计的就是所有长度>=2的期望。
我们从大到小遍历 [ 1 , m ] [1,m] [1,m],设当前枚举到的数为x,我们计算所有长度为 i + 1 ≥ 2 i+1\geq2 i+12且序列的前i项的 g c d gcd gcd为x的期望,那么可以前i项必然都是x的倍数,显然通过求和我们可以得到一个较大的答案

∑ i = 1 ∞ ⌊ m x ⌋ i g e t ( x ) m ( i + 1 ) \sum_{i=1}^{\infty}\lfloor\frac{m}{x}\rfloor^i\frac{get(x)}{m}(i+1) i=1xmimget(x)(i+1)
其中 g e t ( x ) m \frac{get(x)}{m} mget(x)为常数项可以拎出来单独计算,那么后面的 ∑ i = 1 ∞ ⌊ m x ⌋ i ( i + 1 ) \sum_{i=1}^{\infty}\lfloor\frac{m}{x}\rfloor^i(i+1) i=1xmi(i+1)就是我们要求的,这样求和的东西就是高中常说的错位相减法搞一下就能算出来了。
但是为什么说这个答案较大呢?仔细想一下,这样生成的序列的前i项的gcd必然包含所有x的倍数。
P ( d x , d y ) P(dx,dy) P(dx,dy)为序列长度为dx且gcd为dy的概率
那么后半部分真正的答案就是

∑ i = 1 ∞ ( ⌊ m x ⌋ i − ∑ j = i + 1 m P ( i , j ) [ j % i = 0 ] ) ( i + 1 ) = ∑ i = 1 ∞ ⌊ m x ⌋ i ( i + 1 ) − ∑ j = x + 1 m ∑ i = 1 ∞ P ( i , j ) ( i + 1 ) [ j m o d    x = 0 ] \sum_{i=1}^{\infty}(\lfloor\frac{m}{x}\rfloor^i-\sum_{j=i+1}^mP(i,j)[j\%i=0])(i+1)= \sum_{i=1}^\infty\lfloor\frac{m}{x}\rfloor^i(i+1)-\sum_{j=x+1}^m\sum_{i=1}^\infty P(i,j)(i+1)[j\mod x=0] i=1(xmij=i+1mP(i,j)[j%i=0])(i+1)=i=1xmi(i+1)j=x+1mi=1P(i,j)(i+1)[jmodx=0]
我们设 a n [ c ] = ∑ i = 1 ∞ P ( i , c ) ( i + 1 ) an[c]=\sum_{i=1}^{\infty}P(i,c)(i+1) an[c]=i=1P(i,c)(i+1),那么 a n [ x ] = ∑ i = 1 ∞ ⌊ m x ⌋ i ( i + 1 ) − ∑ j = i + 1 m a n [ j ] [ j % i = 0 ] an[x]=\sum_{i=1}^{\infty}\lfloor\frac{m}{x}\rfloor^i(i+1)-\sum_{j=i+1}^man[j][j\% i=0] an[x]=i=1xmi(i+1)j=i+1man[j][j%i=0],那么长度为i+1项且前i项的gcd为x的真正期望就是 a n [ x ] g e t ( x ) m an[x]\frac{get(x)}{m} an[x]mget(x)
综上,对每个数进行讨论并计算答案即可。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2e5 + 10;
#define fi first
#define se second
#define pb push_back
LL m;
const LL mod=1e9+7;
LL pm(LL x,LL y){
  LL z=1;
  while(y){
    if(y&1)z=z*x%mod;
    x=x*x%mod;
    y>>=1;
  }
  return z;
}
int a[N],p[N],cnt;
void P(){
  for(int i=2;i<N;i++){
    if(!p[i])a[++cnt]=i;
    for(int j=1;j<=cnt&&1ll*a[j]*i<N;j++){
      p[a[j]*i]=1;
      if(i%a[j]==0)break;
    }
  }
}
LL get(int y){
  vector<int>x;
  for(int i=1;i<=cnt&&1ll*a[i]*a[i]<=y;i++){
    if(y%a[i]==0){
      x.pb(a[i]);
      while(y%a[i]==0)y/=a[i];
    }
  }
  if(y>1)x.pb(y);
  int sz=x.size();
  int ans=m;
  for(int i=1;i<1<<sz;i++){
    int l=0,r=1;
    for(int j=0;j<sz;j++){
      if(1<<j&i){
        l++;
        r*=x[j];
      }
    }
    if(l&1)ans-=m/r;
    else ans+=m/r;
  }
  return ans;
}
LL _get(LL x){
  return (x*(m-x)+x*m)%mod*pm(1ll*(m-x)*(m-x)%mod,mod-2)%mod;
}
LL an[N];
int main() {
  ios::sync_with_stdio(false);
  cin>>m;P();
  LL re_m=pm(m,mod-2);
  LL ans=re_m;
  for(int i=m;i>=2;i--){
    LL x=get(i);
    LL res=_get(m/i)%mod;
    for(int j=i+i;j<=m;j+=i)res-=an[j],res=(res+mod)%mod;
    an[i]=res;
    ans+=res*x%mod*re_m%mod;
    ans%=mod;
  }
  cout<<ans<<'\n';
  return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值