学习笔记：powerful number求积性函数前缀和

算法原理

本文参考了 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}) 优秀许多。

算法特点

1. 复杂度是 (mathcal O(sqrt n) imes mathcal O( ext{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 << 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;

}