数论 题解|狄利克雷前缀和

先撇开这题不谈,通过这题我倒是发现了推了十几道莫反后连狄利克雷前缀和都没学扎实。

狄利克雷前缀和定义:f1=g,我们称 gf 的狄利克雷前缀和。

给定 f1,f2...fn,考虑如何在 O(nloglogn) 内求出 g 函数。

沿用高维前缀和的思路,我们可以把每一个质因子看作一维,然后套用高维前缀和的求法。

对于二位前缀和,我们是先累加其中一维,然后再累加另一维,形象地来说,如果把二位前缀和理解为一个二维矩阵,就是先对每行求出前缀和。

那么我们可以沿用一样的思路求出狄利克雷前缀和:

for (int i = 1; i <= cnt; ++ i)
  for (int j = 1; j * Prime[i] <= n; ++ j) f[j * Prime[i]] += f[j];

狄利克雷后缀和:

gi=i|ddnfd,我们称 gf 的狄利克雷后缀和。

求法:与狄利克雷前缀和不能说是非常相像,只能说是一模一样

for (int i = 1; i <= cnt; ++ i)
  for (int j = n / Prime[i]; j; -- j) f[j] += f[j * Prime[i]];

倒推狄利克雷前缀和:

如果我们已知 f1=g,如何快速求出 f 函数呢?

解法:前缀和的差分数组就是原数组,我们还是拿二维前缀和举例。首先对行做差分,得到的数组是对应列的前缀和。再做一遍差分即可。

但注意差分要从大到小做。

for (int i = 1; i <= cnt; ++ i)
  for (int j = n / Prime[i]; j; -- j) f[j * Prime[i]] -= f[j * Prime[i]];

然后我们回到这道题:求

i=1nj=1mgcdn(i,j)k=1ij[ik][jk]k

1n,m107

将原式变为

i=1nj=1mgcdn(i,j)k=1ij[ijk]k

然后我们发现这个 k=1ij[ijk]k 很不好动,但如果 ijk,则 ijijk,所以这样的 k 总是成对出现,因此 k=1ij[ijk]k=φ(ij)ij2ij=1 例外,需要特判。
我怎么会告诉你我是打表发现的呢

于是原式化为(我们忽略那个 ÷2ij=1 的特例,因为这无关紧要)

i=1nj=1mgcdn(i,j)φ(ij)ij

仍然不好动,因为 ij 很大。这个时候就只有利用 φ(ij)=φ(i)φ(j)gcd(i,j)φ(gcd(i,j))

然后变成

i=1nj=1mgcdn+1(i,j)φ(gcd(i,j))φ(i)φ(j)ij=d=1ndn+1φ(d)k=1μ(k)i=1ndkφ(idk)idkj=1mdkφ(jdk)jdk

f1(i)=j=1niφ(ij)ij,f2(i)=j=1miφ(ij)ij,这俩玩意儿就是个狄利克雷后缀和,可以快速预处理。

枚举 dk,化为

i=1nf1(i)f2(i)d|idn+1φ(d)μ(id)

f(i)=in+1φ(i),g(i)=d|idn+1φ(d)μ(id),显然 g=fμ

狄利克雷卷积要 O(nlogn) 求,无法通过本题。但 g=fμ 意味着 f=g1,也就是倒推狄利克雷前缀和。综上,总时间复杂度 O(nloglogn)

然后别忘了最后有个 +1÷2

#include <cstdio>

const int mod = 998244353, N = 1e7 + 5;
int qpow(int a, int b) {
	int ret = 1ll;
	while (b) {
		if (b & 1) ret = 1ll * ret * a % mod;
		a = 1ll * a * a % mod, b >>= 1;
	}
	return ret;
}
int phi[N], inv[N], Prime[N / 10], cnt, f1[N], f2[N], g[N], n, m;
bool mark[10000005];
inline void add(int &x, const int y) {
	if ((x += y) >= mod) x -= mod;
}
inline void mns(int &x, const int y) {
	if ((x -= y) < 0) x += mod;
}

void init(int N) {
	g[1] = inv[1] = phi[1] = 1;
	for (int i = 2; i <= N; ++ i) {
		if (!mark[i]) Prime[++ cnt] = i, phi[i] = i - 1, g[i] = qpow(i, n + 1), inv[i] = qpow(i, mod - 2);
		for (int j = 1; i * Prime[j] <= N; ++ j) {
			mark[i * Prime[j]] = true, g[i * Prime[j]] = 1ll * g[i] * g[Prime[j]] % mod;
			inv[i * Prime[j]] = 1ll * inv[i] * inv[Prime[j]] % mod;
			if (i % Prime[j] == 0) {phi[i * Prime[j]] = phi[i] * Prime[j]; break;}
			phi[i * Prime[j]] = phi[i] * (Prime[j] - 1);
		}
	}
	for (int i = 1; i <= n; ++ i) g[i] = 1ll * g[i] * inv[phi[i]] % mod;
	for (int i = 1; i <= n; ++ i) f1[i] = 1ll * phi[i] * i % mod;
	for (int i = 1; i <= m; ++ i) f2[i] = 1ll * phi[i] * i % mod;
	for (register int i = cnt; i; -- i) 
		for (register int j = n / Prime[i]; j; -- j) mns(g[j * Prime[i]], g[j]);
	for (register int i = 1; i <= cnt && Prime[i] <= n; ++ i)
		for (register int j = n / Prime[i]; j; -- j) add(f1[j], f1[j * Prime[i]]);
	for (register int i = 1; i <= cnt && Prime[i] <= m; ++ i)
		for (register int j = m / Prime[i]; j; -- j) add(f2[j], f2[j * Prime[i]]);
}

int main() {
	scanf("%d%d", &n, &m);
	init(n > m ? n : m);
	int ans = 0;
	for (int i = 1; i <= n; ++ i) add(ans, 1ll * f1[i] * f2[i] % mod * g[i] % mod);
	printf("%d", ((int)((1ll * ans + 1) * qpow(2, mod - 2) % mod) + mod) % mod);
	return 0;
}

本文作者:zqs2020

本文链接:https://www.cnblogs.com/stinger/p/15778463.html

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   zqs2020  阅读(830)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起
  1. 1 404 not found REOL
404 not found - REOL
00:00 / 00:00
An audio error has occurred.