算法原理

本文参考了 zzq's blog 。

\(\text{powerful number}\) 的定义是每个质因子次数都 \(\ge 2\) 的数，有个结论是 \(\ge n\) 的 \(\text{powerful number}\) 只有 \(\mathcal O(\sqrt n)\) 个，如何找这些数呢？用暴力 \(\text{dfs}\) 从小到达枚举质因子及其幂次即可（类似于 \(\text{min_25}\) 第二部分）。

比如对于函数 \(F(p^q) = p^k\) 其中 \(p\) 为素数且 \(k\) 为定值，且 \(f(x)\) 是积性函数。我们需要求 \(\sum_{i = 1}^{n} F(i)\)。

考虑令 \(G(x) = x^k\) ，令 \(\displaystyle H = \frac F G\) （其中除法为狄利克雷卷积的逆运算），由于 \(F, G\) 都为积性函数，所以 \(H\) 也为积性函数。

我们考虑求出 \(H\) ，有 \(F(p) = G(p)H(1) + H(p)G(1)\) 由于 \(F(p) = G(p)\) 且 \(H(1) = 1, G(1) = 1\) 易求（我们通常令 \(F(1) = 1\) ），那么有 \(H(p) = 0\) ，又由于 \(H\) 为积性函数，所以 \(H(x)\) 只有当 \(x\) 为 \(\text{powerful number}\) 时有值。

有什么用呢？我们考虑原来的式子 \(\displaystyle\sum_{i = 1}^{n} F(i) = \sum_{ij \le n} H(i) G(j) = \sum_{i = 1}^{n} H(i) \sum_{j = 1}^{\lfloor \frac n i\rfloor} G(j)\) 。

现在只剩下一个问题，如何求 \(H(x)\) ，由于是积性函数，我们只需要求出 \(H(p^q)\) ，可以归纳出 \(H(p^q) = p^{k} - p^{2k}\) （读者自证不难）。

但是对于通用的 \(H(x)\) 如何求呢？我们考虑对于 \(p\) 的指数 \(q\) 来说，等价于多项式求逆，可以 \(\mathcal O(q)\) 递推一项。

然后我们可以利用插值等求自然数幂和的方式在 \(\mathcal O(k\sqrt n)\) 的时间求出对应的前缀和，比 \(\text{min_25}\) 优秀许多。

算法特点

1. 复杂度是 \(\mathcal O(\sqrt n) \times \mathcal O(\text{calcG})\) ，所以大部分时候有显著的时间优势，而且十分简短好记。
2. 但是 \(G(x)\) 通常十分难以构造，注意是我们要令 \(G(p) = F(p)\) 且 \(G(x)\) 为积性函数，有十分大的局限性。

实现

#include <bits/stdc++.h>#define For(i, l, r) for(register int i = int(l), i##end = int(r); i <= i##end; ++ i)
#define Fordown(i, r, l) for(register int i = int(r), i##end = int(l); i >= i##end; -- i)
#define Rep(i, r) for(register int i = int(0), i##end = int(r); i < i##end; ++ i)
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define Set(a, b) memset(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << x << endlusing namespace std;typedef long long ll;template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }template<typename T = int>
inline T read() {T x(0), sgn(1); char ch(getchar());for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);return x * sgn;
}void File() {freopen ("function.in", "r", stdin);freopen ("function.out", "w", stdout);
}const int N = 1e7 + 1e3, Mod = 1e9 + 7, K = 25;ll n; int k;inline int fpm(int x, int power) {int res(1);for (; power; power >>= 1, x = 1ll * x * x % Mod)if (power & 1) res = 1ll * res * x % Mod;return res;
}bool is_prime[N]; int prime[N], prep[N], powp[N], pcnt;void Linear_Sieve(int maxn) {Set(is_prime, true); is_prime[0] = is_prime[1] = false;For (i, 2, maxn) {if (is_prime[i]) {prime[++ pcnt] = i;prep[pcnt] = (prep[pcnt - 1] + (powp[pcnt] = fpm(i, k))) % Mod;}for (int j = 1, res; j <= pcnt && (res = prime[j] * i) <= maxn; ++ j) {is_prime[res] = false; if (!(i % prime[j])) break;}}
}int pre[K], suf[K], fac[K], ifac[K], y[K];void Fac_Init(int maxn) {fac[0] = ifac[0] = 1;For (i, 1, maxn) y[i] = (y[i - 1] + fpm(i, k)) % Mod;For (i, 1, maxn) fac[i] = 1ll * fac[i - 1] * i % Mod;ifac[maxn] = fpm(fac[maxn], Mod - 2);Fordown (i, maxn - 1, 1) ifac[i] = ifac[i + 1] * (i + 1ll) % Mod;
}inline int Sumk(int m) {int maxn = k + 2, ans = 0;pre[0] = suf[maxn + 1] = 1;For (i, 1, maxn) pre[i] = 1ll * pre[i - 1] * (m - i) % Mod;Fordown (i, maxn, 1) suf[i] = 1ll * suf[i + 1] * (m - i) % Mod;For (i, 1, maxn) {int coef = 1ll * pre[i - 1] * suf[i + 1] % Mod, inv = ((maxn - i) & 1 ? -1ll : 1ll) * ifac[i - 1] * ifac[maxn - i] % Mod;ans = (ans + 1ll * coef * y[i] % Mod * inv) % Mod;}return ans;
}ll val[N]; int sum[N];int cnt, id1[N], id2[N], d, res[N];inline int &id(ll x) {return x <= d ? id1[x] : id2[n / x];
}int dcnt = 0;void Dfs(ll x, int cur, int coef) {(sum[id(n / x)] += coef) %= Mod;for (int i = cur + 1; i <= pcnt && x <= n / prime[i] / prime[i]; ++ i) {ll y = prime[i];do {y *= prime[i];Dfs(x * y, i, (powp[i] - 1ll * powp[i] * powp[i]) % Mod * coef % Mod);} while (y <= n / x / prime[i]);}
}int main() {File();n = read<ll>(); k = read();Fac_Init(k + 2); Linear_Sieve(d = sqrt(n + .5));for (ll i = 1; i <= n; i = n / (n / i) + 1) val[id(n / i) = ++ cnt] = n / i;Dfs(1, 0, 1);int ans = 0;For (i, 1, cnt) if (sum[i])ans = (ans + 1ll * sum[i] * Sumk(val[i] % Mod)) % Mod;printf ("%d\n", (ans + Mod) % Mod);return 0;}