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

算法原理

本文参考了 zzq's blog 。

$\text{powerful number}$ 的定义是每个质因子次数都 $\ge 2$ 的数，有个结论是 $\le 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 << 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\n", (ans + Mod) % Mod);

	return 0;

}

__EOF__

本文作者：zjp_shadow
本文链接：https://www.cnblogs.com/zjp-shadow/p/11093242.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！