算法原理
本文参考了 zzq's blog 。
( ext{powerful number}) 的定义是每个质因子次数都 (ge 2) 的数,有个结论是 (le n) 的 ( ext{powerful number}) 只有 (mathcal O(sqrt n)) 个,如何找这些数呢?用暴力 ( ext{dfs}) 从小到达枚举质因子及其幂次即可(类似于 ( ext{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) 为 ( ext{powerful number}) 时有值。
有什么用呢?我们考虑原来的式子 (displaystylesum_{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 floor} G(j)) 。
现在只剩下一个问题,如何求 (H(x)) ,由于是积性函数,我们只需要求出 (H(p^q)) ,可以归纳出 (H(p^q) = p^{k} - p^{2k}) (读者自证不难)。
但是对于通用的 (H(x)) 如何求呢?我们考虑对于 (p) 的指数 (q) 来说,等价于多项式求逆,可以 (mathcal O(q)) 递推一项。
然后我们可以利用插值等求自然数幂和的方式在 (mathcal O(ksqrt n)) 的时间求出对应的前缀和,比 ( ext{min_25}) 优秀许多。
算法特点
- 复杂度是 (mathcal O(sqrt n) imes mathcal O( ext{calcG})) ,所以大部分时候有显著的时间优势,而且十分简短好记。
- 但是 (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 << endl
using 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
", (ans + Mod) % Mod);
return 0;
}