参考博客:https://www.90yang.com/2019ccpc-xiangtan-f-neko-and-function/
思维比较巧妙,首先从贡献的角度计算 ,枚举 n 的每一个质因子,设其幂次为 ,问题等价于将这 个 放到 个盒子里。
由于盒子不允许为空,这个贡献很难计算,考虑盒子可以为空,用容斥得到最后的答案:令 为 去掉 的限制。
没有这个限制,每一个素数的贡献独立,即 是一个积性函数
考虑容斥得到
:枚举有几个盒子为空,容斥一下即可:
由于
是积性函数,考虑用min_25筛得到
由于这题卡常,我最终也没有卡过去,待update
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e6 + 10;
typedef long long ll;
const int mod = 1e9 + 7;
int num,tot,k;
bool ispri[maxn];
int sum[maxn],pri[maxn],w[maxn],g[maxn],id1[maxn],id2[maxn];
int n,sqr,fact[maxn],ifact[maxn],ans[maxn];
inline void sieve(int x) {
num = 0; ispri[0] = ispri[1] = true;
for (int i = 2; i <= x; i++) {
if (!ispri[i]) {
pri[++num] = i;
sum[num] = (sum[num - 1] + 1) % mod;
}
for (int j = 1; j <= num && i * pri[j] <= x; j++) {
ispri[i * pri[j]] = true;
if (i % pri[j] == 0) break;
}
}
}
int fpow(int a,int b) {
int r = 1;
while (b) {
if (b & 1) r = 1ll * r * a % mod;
b >>= 1;
a = 1ll * a * a % mod;
}
return r;
}
inline int C(int x,int y) {
if (y > x || y < 0) return 0;
return 1ll * fact[x] * ifact[y] % mod * ifact[x - y] % mod;
}
inline int S(ll x,int y,int z) {
if (pri[y] >= x) return 0;
int k = x <= sqr ? id1[x] : id2[n / x];
int ans = 1ll * z * (g[k] - sum[y] + mod) % mod;
for (int i = y + 1; i <= num && pri[i] <= x / pri[i]; i++)
for (ll e = 1, pe = pri[i]; pe <= x; e++, pe = pe * pri[i])
ans = (ans + 1ll * (S(x / pe,i,z) + (e != 1)) % mod * C(e + z - 1,z - 1) % mod) % mod;
return ans;
}
int main() {
int N = maxn - 10;
sieve(N);
fact[0] = 1;
for (int i = 1; i <= N; i++)
fact[i] = 1ll * fact[i - 1] * i % mod;
ifact[N] = fpow(fact[N],mod - 2);
for (int i = N - 1; i >= 0; i--)
ifact[i] = 1ll * ifact[i + 1] * (i + 1) % mod;
while(scanf("%d%d",&n,&k) != EOF) {
if (k == 1) {
printf("%d\n",(n - 1) % mod);
continue;
}
else if (n < (1<<k)) {
puts("0");
continue;
}
else if (n < (1 << k) + (1 << (k - 1))) {
puts("1");
continue;
}
tot = 0;
sqr = sqrt(n);
for (int i = 1,j; i <= n; i = j + 1) {
j = n / (n / i);
w[++tot] = n / i;
int p = n / i % mod;
g[tot] = (p - 1) % mod;
if (n / i <= sqr) id1[n / i] = tot;
else id2[j] = tot;
}
for (int i = 1; i <= num; i++) {
for (int j = 1; j <= tot && pri[i] <= w[j] / pri[i]; j++) {
int x = w[j] / pri[i];
x = x <= sqr ? id1[x] : id2[n / x];
g[j] -= (g[x] - sum[i - 1] + mod) % mod;
if (g[j] < 0) g[j] += mod;
}
}
int res = 0, t = 1;
for (int i = k; i >= 0; i--) {
if (k == 0) res = (res + C(k,i) * t % mod) % mod;
else if (k == 1) res = (res + 1ll * n * C(k,i) % mod * t % mod) % mod;
else res = (res + 1ll * S(n,0,i) * C(k,i) % mod * t % mod) % mod;
if (res < 0) res += mod;
t *= -1;
}
printf("%d\n",res);
}
return 0;
}