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

· · 算法·理论

开篇声明一下，本文是对 _fewq \mathcal O(n(\log\log n)^2) Dirichlet 卷积算法的详细解读。

集合 \text{SF},\text{PN} 分别指所有 square-free number, powerful number 构成的集合。

定义函数 \omega(i)=\sum_{p|i}1。

Square free 处的 Dirichlet 卷积计算

对于数列 F_{i,j}(i\in\text{SF},0\le j\le \omega(i))，定义两种变换 \mathsf P,\mathsf Q 为：

(\mathsf P F)_{n,k}=\sum_{d|n}F_{d,k} (\mathsf Q F)_{n,k}=\sum_{d|n}F_{d,k+\omega(n/d)}

这两种变换都是可逆的，我们这里只关心 \mathsf Q 的逆变换，它的表达式为：

(\mathsf Q^{-1}F)_{n,k}=\sum_{d|n}F_{d,k+\omega(n/d)}\mu(n/d)

形象地说，\mathsf P,\mathsf Q^{-1} 可以分别视为 DFT 和 IDFT，所以我们现在还差一步「点值相乘」的操作。因此，我们定义数列 F_{i,j},G_{i,j} 的乘法为：

(F\otimes G)_{n,k}=\sum_{i+j=\omega(n)+k}F_{n,i}G_{n,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}

如果令 F_{n,k}=f(n)[k=\omega(n)],G_{n,k}=g(n)[k=\omega(n)]，就有：

\mathsf Q^{-1}(\mathsf P F\otimes \mathsf P G)_{n,0}=\sum_{st=n}f(s)g(t)

这恰是 sqaure free 上的 Dirichlet 卷积。

而变换 \mathsf P,\mathsf Q^{-1} 和运算 \otimes 的时间复杂度均为 \mathcal O(n(\log \log n)^2)，故我们能在 \mathcal O(n(\log \log n)^2) 时间内计算 square free 处的 Dirichlet 卷积。

补全其他部分的 Dirichlet 卷积

对于正整数 n，我们熟知其有唯一的 PN-SF 分解。即存在唯一的数对 (d,s)，使得 n=ds,d\perp s,d\in\text{PN},s\in\text{SF}。我们上面已经展示了当 d=1 时，如何计算 (f*g)(n)，下面考虑 d>1 的情况。

记 h=f*g，则当 ij=n 时，f(i)g(j) 对 h(n) 有贡献。但是因为 n 不是 square free 的，我们没法用上面的方法。不过，设 v=\gcd(i,d)，则 (i/v)(j/(d/v))=n/d，此时 n/d 是 square free 的，就可以设

f_*(i)=f(iv)[i\perp d] g_*(i)=g(id/v)[i\perp d]

然后去计算 f_*,g_* 的 square free 处的 Dirichlet 卷积，并把 s 处的值累加到 h(n) 上。

总的来说，我们只需枚举所有可能的 (d,v)（PN 及其因子），按上面的方法计算出所有 n/d 为 square free，v=\gcd(i,d) 的贡献。总时间复杂度为：

\sum_{d\in\text{PN}}\sigma_0(d)(n/d)[\log\log(n/d)]^2=\mathcal O(n(\log\log n)^2)

一个代码实现：

#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 的地方可以线性预处理，\mathcal O(1) 查询。因为若记 r(n)=\prod_{p|n}p，则 PN 的 r 值小于等于 \sqrt n，我们只需预处理出所有 \gcd(i,j),i,j\le \sqrt n 即可。这样比原先 _fewq 的实现更省空间。

我们得到了什么

结合我之前的 DGF 牛顿迭代理论，我们可以得到以下问题的 \mathcal O(n(\log \log n)^2) 解法。

DGF 或将焕发第二春！