md我连exgcd都不会了我好菜啊
考虑首先 n n n为质数时
考虑如果前
k
−
1
k-1
k−1个都随便取
那最后一个就固定了选哪个
但是这样有可能和前面重复
考虑容斥
a
n
s
=
1
k
∑
i
=
1
k
(
−
1
)
i
−
1
C
n
k
−
i
ans=\frac 1 k\sum_{i=1}^{k} (-1)^{i-1}C_n^{k-i}
ans=k1∑i=1k(−1)i−1Cnk−i
由于
n
n
n为质数,所以我们枚举相同的最后选了
i
i
i次,前面
k
−
i
k-i
k−i随便选
相当于解同余方程,有唯一解
考虑当
n
n
n不为质数时
相当于解方程
x
i
+
s
=
p
y
xi+s=py
xi+s=py
要求
s
∣
g
c
d
(
i
,
p
)
s|gcd(i,p)
s∣gcd(i,p)
现在要求的是有多少个
x
x
x满足这个方程
相当于求
x
x
x的解的个数
考虑
x
x
x变
p
g
c
d
(
i
,
p
)
\frac p {gcd(i,p)}
gcd(i,p)p,
y
y
y变
i
g
c
d
(
i
,
p
)
\frac i {gcd(i,p)}
gcd(i,p)i
恰好有一组正整数解
x
∈
[
0
,
n
−
1
]
x\in[0,n-1]
x∈[0,n−1]
有
n
n
n个取值
所以方案为
n
p
g
c
d
(
i
,
p
)
=
n
/
p
∗
g
c
d
(
i
,
p
)
\frac {n}{\frac{p}{gcd(i,p)}}=n/p*gcd(i,p)
gcd(i,p)pn=n/p∗gcd(i,p)
设 f [ p ] [ k ] f[p][k] f[p][k]表示选 k k k个为 p p p的倍数的方案数
则 f [ p ] [ k ] = ∑ i = 1 k ( − 1 ) k − i n / p ∗ g c d ( i , p ) f [ g c d ( i , p ) ] [ k − i ] f[p][k]=\sum_{i=1}^{k}(-1)^{k-i}n/p*gcd(i,p)f[gcd(i,p)][k-i] f[p][k]=∑i=1k(−1)k−in/p∗gcd(i,p)f[gcd(i,p)][k−i]
递归求解即可
#include<bits/stdc++.h>
using namespace std;
#define gc getchar
inline int read(){
char ch=gc();
int res=0,f=1;
while(!isdigit(ch))f^=ch=='-',ch=gc();
while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
return f?res:-res;
}
#define re register
#define pb push_back
#define cs const
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
const int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Dec(int &a,int b){(a-=b)<0?(a+=mod):0;}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){
for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;
}
cs int N=5005;
int inv[N],fac[N],ifac[N],c[N],n,k;
inline void init(){
inv[1]=1;
for(int i=2;i<N;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
fac[0]=ifac[0]=1;
for(int i=1;i<N;i++)fac[i]=mul(fac[i-1],i);
ifac[N-1]=ksm(fac[N-1],mod-2);
for(int i=N-2;i;i--)ifac[i]=mul(ifac[i+1],i+1);
c[0]=1,c[1]=n;
for(int i=2;i<=k;i++)c[i]=mul(c[i-1],mul(n-i+1,inv[i]));
}
int f[N][N];
inline int gcd(int x,int y){
return y?gcd(y,x%y):x;
}
int dfs(int p,int k){
if(k==0)return 1;
if(p==1)return c[k];
if(p<N&&f[p][k]!=-1)return f[p][k];
int res=0;
for(int i=1;i<=k;i++){
int g=gcd(i,p);
if(i&1)Add(res,mul(n/p,mul(g,dfs(g,k-i))));
else Dec(res,mul(n/p,mul(g,dfs(g,k-i))));
}
Mul(res,inv[k]);
if(p>=N)return res;
return f[p][k]=res;
}
int main(){
n=read(),k=read();
init(),memset(f,-1,sizeof(f));
cout<<dfs(n,k);
}