专栏文章

DGF 不存在了：Dirichlet 卷积的 n(loglogn)^2 做法

飞飞雨烟雁
算法·理论参与者 47已保存评论 51

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
51 条
当前快照
1 份
快照标识符
@mipt2fdw
此快照首次捕获于
2025/12/03 17:29
3 个月前
此快照最后确认于
2025/12/03 17:29
3 个月前
查看原文
开篇声明一下，本文是对 _fewq O(n(log⁡log⁡n)2)\mathcal O(n(\log\log n)^2) Dirichlet 卷积算法的详细解读。
集合 SF,PN\text{SF},\text{PN} 分别指所有 square-free number, powerful number 构成的集合。
定义函数 ω(i)=∑p∣i1\omega(i)=\sum_{p|i}1。

Square free 处的 Dirichlet 卷积计算

对于数列 Fi,j(i∈SF,0≤j≤ω(i))F_{i,j}(i\in\text{SF},0\le j\le \omega(i))，定义两种变换 P,Q\mathsf P,\mathsf Q 为：
(PF)n,k=∑d∣nFd,k(\mathsf P F)_{n,k}=\sum_{d|n}F_{d,k}
(QF)n,k=∑d∣nFd,k+ω(n/d)(\mathsf Q F)_{n,k}=\sum_{d|n}F_{d,k+\omega(n/d)}
这两种变换都是可逆的，我们这里只关心 Q\mathsf Q 的逆变换，它的表达式为：
(Q−1F)n,k=∑d∣nFd,k+ω(n/d)μ(n/d)(\mathsf Q^{-1}F)_{n,k}=\sum_{d|n}F_{d,k+\omega(n/d)}\mu(n/d)
形象地说，P,Q−1\mathsf P,\mathsf Q^{-1} 可以分别视为 DFT 和 IDFT，所以我们现在还差一步「点值相乘」的操作。因此，我们定义数列 Fi,j,Gi,jF_{i,j},G_{i,j} 的乘法为：
(F⊗G)n,k=∑i+j=ω(n)+kFn,iGn,j(F\otimes G)_{n,k}=\sum_{i+j=\omega(n)+k}F_{n,i}G_{n,j}
那么我们就有如下定理：
Q−1(PF⊗PG)n,k=∑i+j=ω(n)+k∑lcm(s,t)=nFs,iGt,j\mathsf Q^{-1}(\mathsf P F\otimes \mathsf P G)_{n,k}=\sum_{i+j=\omega(n)+k}\sum_{\text{lcm}(s,t)=n}F_{s,i}G_{t,j}
如果令 Fn,k=f(n)[k=ω(n)],Gn,k=g(n)[k=ω(n)]F_{n,k}=f(n)[k=\omega(n)],G_{n,k}=g(n)[k=\omega(n)]，就有：
Q−1(PF⊗PG)n,0=∑st=nf(s)g(t)\mathsf Q^{-1}(\mathsf P F\otimes \mathsf P G)_{n,0}=\sum_{st=n}f(s)g(t)
这恰是 sqaure free 上的 Dirichlet 卷积。
而变换 P,Q−1\mathsf P,\mathsf Q^{-1} 和运算 ⊗\otimes 的时间复杂度均为 O(n(log⁡log⁡n)2)\mathcal O(n(\log \log n)^2)，故我们能在 O(n(log⁡log⁡n)2)\mathcal O(n(\log \log n)^2) 时间内计算 square free 处的 Dirichlet 卷积。

补全其他部分的 Dirichlet 卷积

对于正整数 nn，我们熟知其有唯一的 PN-SF 分解。即存在唯一的数对 (d,s)(d,s)，使得 n=ds,d⊥s,d∈PN,s∈SFn=ds,d\perp s,d\in\text{PN},s\in\text{SF}。我们上面已经展示了当 d=1d=1 时，如何计算 (f∗g)(n)(f*g)(n)，下面考虑 d>1d>1 的情况。
记 h=f∗gh=f*g，则当 ij=nij=n 时，f(i)g(j)f(i)g(j) 对 h(n)h(n) 有贡献。但是因为 nn 不是 square free 的，我们没法用上面的方法。不过，设 v=gcd⁡(i,d)v=\gcd(i,d)，则 (i/v)(j/(d/v))=n/d(i/v)(j/(d/v))=n/d，此时 n/dn/d 是 square free 的，就可以设
f∗(i)=f(iv)[i⊥d]f_*(i)=f(iv)[i\perp d] g∗(i)=g(id/v)[i⊥d]g_*(i)=g(id/v)[i\perp d]
然后去计算 f∗,g∗f_*,g_* 的 square free 处的 Dirichlet 卷积，并把 ss 处的值累加到 h(n)h(n) 上。
总的来说，我们只需枚举所有可能的 (d,v)(d,v)（PN 及其因子），按上面的方法计算出所有 n/dn/d 为 square free，v=gcd⁡(i,d)v=\gcd(i,d) 的贡献。总时间复杂度为：
∑d∈PNσ0(d)(n/d)[log⁡log⁡(n/d)]2=O(n(log⁡log⁡n)2)\sum_{d\in\text{PN}}\sigma_0(d)(n/d)[\log\log(n/d)]^2=\mathcal O(n(\log\log n)^2)
一个代码实现：
CPP
#include <iostream>
#include <cstdio>
#include <cmath>
#define ll long long
#define int unsigned int
using namespace std;
const int Mx = 1e6 + 10, Nx = log(Mx) / log(2);
bool vis[Mx], squarefree[Mx];
int prime[Mx], tot;
int omega[Mx]; // for square free  
void sieve(){
	squarefree[1] = 1;
	for(int i = 2; i < Mx; ++i){
		if(!vis[i]) prime[++tot] = i, squarefree[i] = omega[i] = 1;
		for(int j = 1; j <= tot && prime[j] * i < Mx; ++j){
			vis[i * prime[j]] = true;
			if(i % prime[j] == 0) break;
			omega[i * prime[j]] = omega[i] + 1;
			squarefree[i * prime[j]] = squarefree[i];
		}
	}
}

void FDT(int F[][Nx], int n){ // fast dirichlet transform
	for(int i = 1; prime[i] <= n && i <= tot; ++i){
		int p = prime[i];
		for(int j = n / p; j; --j) if(squarefree[j * p]){
			for(int r = 0; r <= omega[j]; ++r) F[j * p][r] += F[j][r];
		}
	}
}

void IFDT(int F[][Nx], int n){ // inverse fast dirichlet transform
	for(int i = 1; prime[i] <= n && i <= tot; ++i){
		int p = prime[i];
		for(int j = 1; j <= n / p; ++j) if(squarefree[j * p])
			for(int r = 1; r <= omega[j]; ++r) F[j * p][r - 1] -= F[j][r];
	}
}

int F[Mx][Nx], G[Mx][Nx], H[Nx];
void mul(int f[], int g[], int h[], int n){ // h = f * g for square free
	for(int i = 1; i <= n; ++i) if(squarefree[i]){
		F[i][omega[i]] = f[i], G[i][omega[i]] = g[i];
		for(int r = 0; r < omega[i]; ++r) F[i][r] = G[i][r] = 0;
	}
	FDT(F, n), FDT(G, n);
	for(int i = 1; i <= n; ++i) if(squarefree[i]){
		for(int r = 0; r <= omega[i]; ++r) H[r] = 0;
		for(int r = 0; r <= omega[i]; ++r) for(int s = omega[i] - r; s <= omega[i]; ++s)
			H[r + s - omega[i]] += F[i][r] * G[i][s];
		for(int r = 0; r <= omega[i]; ++r) F[i][r] = H[r];
	}
	IFDT(F, n);
	for(int i = 1; i <= n; ++i) h[i] = F[i][0];
}

int gcd(int x, int y){ return y ? gcd(y, x % y) : x;}

int f[Mx], g[Mx], h[Mx];
int fd[Mx], gd[Mx], hd[Mx];
void dfs(int A, int B, int pos, int n){ // A * B is PN
	if((ll)prime[pos] * prime[pos] > n / A / B || pos > tot){
		int m = n / A / B, pn = A * B;
		// 下面用 gcd 来写其实会使复杂度变劣
		// 严格的写法需要按 fewq 那种写法，每次 dfs 到下一层时就分离掉不互质的部分 
		// 或者使用快速 gcd，毕竟这样不需要额外的空间，对常数也许有帮助
		for(int i = 1; i <= m; ++i){
			if(gcd(i, pn) == 1) fd[i] = f[i * A], gd[i] = g[i * B];
			else fd[i] = gd[i] = 0;
		}
		mul(fd, gd, hd, m);
		for(int i = 1; i <= m; ++i) h[i * pn] += hd[i];
		return;
	}
	dfs(A, B, pos + 1, n);
	int m = n / A / B / prime[pos] / prime[pos];
	int prod = prime[pos] * prime[pos];
	for(int i = 2; m; ++i, m /= prime[pos], prod *= prime[pos])
		for(int j = 0, k = 1; j <= i; ++j, k *= prime[pos])
			dfs(A * k, B * (prod / k), pos + 1, n);
}

signed main(){
	sieve();
	int n;
	cin >> n;
	for(int i = 1; i <= n; ++i) cin >> f[i];
	for(int i = 1; i <= n; ++i) cin >> g[i];
	dfs(1, 1, 1, n);
	for(int i = 1; i <= n; ++i) cout << h[i] << ' ';
	return 0;
}
注：里面调用 gcd 的地方可以线性预处理，O(1)\mathcal O(1) 查询。因为若记 r(n)=∏p∣npr(n)=\prod_{p|n}p，则 PN 的 rr 值小于等于 n\sqrt n，我们只需预处理出所有 gcd⁡(i,j),i,j≤n\gcd(i,j),i,j\le \sqrt n 即可。这样比原先 _fewq 的实现更省空间。

我们得到了什么

结合我之前的 DGF 牛顿迭代理论，我们可以得到以下问题的 O(n(log⁡log⁡n)2)\mathcal O(n(\log \log n)^2) 解法。
  • DGF 乘法逆
  • DGF ln / exp
  • DGF 快速幂 / 开根
DGF 或将焕发第二春！

评论

共 51 条评论，欢迎与作者交流。

正在加载评论...