题目链接
前置推论:
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+1≥2且序列的前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=1∑∞⌊xm⌋imget(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=1∞⌊xm⌋i(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∑∞(⌊xm⌋i−j=i+1∑mP(i,j)[j%i=0])(i+1)=i=1∑∞⌊xm⌋i(i+1)−j=x+1∑mi=1∑∞P(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=1∞P(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=1∞⌊xm⌋i(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;
}