跳转至

狄利克雷卷积 & 莫比乌斯反演

定义两个数论函数 \(f(x),g(x)\) 的狄利克雷卷积 \((f*g)(x)\) 为

\[ (f*g)(x)=\sum_{i|x}f(i)g(\frac xi) \]

性质 1：狄利克雷卷积具有交换律，结合律；
性质 2：其单位元 \(\epsilon(x)\) 为 \(\epsilon(x)=[x=1]\)；
性质 3：两个积性函数的狄利克雷卷积还是积性函数；
性质 4：若一个数论函数 \(f\) 满足 \(f(1)\ne 0\)，那么其存在逆元；

现在介绍几个重要的数论函数：

\(\operatorname I(x)=1\)：常函数；
\(\operatorname{Id}(x)=x\)：线性函数（我瞎叫的）；
\(\operatorname{Id}^k(x)=x^k\)：幂函数；
\(\varphi(x)\)：不用多说了吧；
\(\mu(x)\)：莫比乌斯函数；
\(d(x)\)：因数个数函数；
\(\sigma^k(x)\)：因数（\(k\) 次方和）和函数；

这里说一下莫比乌斯函数，对于正整数 \(n=\prod{p_i^{\alpha_i}}\)（质因数分解），\(\mu(x)\) 定义为：

\[ \mu(x)= \begin{cases} 0,&\exist\alpha_i\ge 2\\ (-1)^{\sum[\alpha_i=1]},&\forall\alpha_i\in\{0,1\} \end{cases} \]

不难证明 \(\mu(x)\) 是积性函数，因此可以使用线性筛。

推论 1：\(\mu*\operatorname{I}=\epsilon\)（莫比乌斯反演）；
推论 2：\(\varphi*\operatorname{I}=\operatorname{Id}\)（欧拉反演）；
推论 3：\(\operatorname{Id}^k*\operatorname{I}=\sigma^k\)（废话）；

证明：

（推论 1）由于 \(\mu*\operatorname{I}\) 还是积性函数，因此只需证明 \((\mu*\operatorname{I})(p^k)=[k=0]\)：

\[ \sum_{i|p^k}\mu(i)=\sum_{i=0}^{k}\mu(p^i)= \begin{cases} 1+(-1)+0+0+\cdots,&(k\ge 1)\\ 1,&(k=0) \end{cases} \]

（推论 2）考虑拆贡献，把 \(1\sim n\) 的每一个数字 \(x\) 放到 \(\gcd(x,n)\) 处统计：

\[ \sum_{i=1}^{n}1=\sum_{i|n}\sum_{i|j\wedge j\le n}[\frac ji\bot\frac ni]=\sum_{i|n}\varphi(\frac ni) \]

将数字 \(n\) 质因数分解为 \(\prod p_i^{\alpha_i}\)，然后就可以写成一个无穷维向量 \((\alpha_1,\alpha_2,\alpha_3,\alpha_4,\cdots)\) 分别表示每个质因子的次数。那么数论函数就是一个向量到数值的映射，卷 \(\operatorname{I}\) 对应着高维前缀和，卷 \(\mu\) 对应高维差分。

狄利克雷前缀和 & 差分

给定数论函数 \(f(x)\)，如何快速计算 \(f*\operatorname{I}\) 和 \(f*\mu\) 呢？暴力卷积显然是 \(O(n\ln n)\) 的。注意到高维前缀和可以逐维度考虑，依次进行前缀和。因此枚举质数 \(p\)，然后执行 \(f(i)\to f(pi)\)（箭头表示贡献）。这样只有 \(\lfloor\frac np\rfloor\) 个 \(i\) 会被处理，因此时间复杂度同埃氏筛，为 \(O(n\ln\ln n)\)。

P2522 [HAOI2011] Problem b

给定 \(n,m\)，求

\[ \sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=1] \]

\(n,m\le 5\times 10^4\)，多组询问，询问次数 \(\le 5\times 10^4\)

考虑到 \([x=1]=\sum\limits_{i|x}\mu(i)\)，因此上式可以写成

\[ \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|\gcd(i,j)}\mu(d)=\sum_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\frac{n}{d}\right\rfloor\left\lfloor\frac{m}{d}\right\rfloor \]

筛出 \(\mu(x)\) 后整除分块即可。

P5518 [MtOI2019] 幽灵乐团 / 莫比乌斯反演基础练习题

对于 \(f(i,j,k)=1,~ijk,~\gcd(i,j,k)\) 分别求出

\[ \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}\biggl(\frac{[i,j]}{(i,k)}\bigg)^{f(i,j,k)} \]

多测，\(a,b,c\le 10^5,~T\le 70\)；对质数 \(p\) 取模，\(10^7\le p\le 1.05\times 10^9\)。

\[ \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}\biggl(\frac{ij}{(i,j)(i,k)}\bigg)^{f(i,j,k)} \]

分子分母分别求，\(i,j\) 和 \((i,j),(i,k)\) 也可以分离，问题转化为

\[ \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}i^{f(i,j,k)}\\ \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}(i,j)^{f(i,j,k)} \]

总共 \(6\) 个数分别记为 \(A,B,C,D,E,F,G\)。先看第一行，\(f(i,j,k)=1\) 和 \(ijk\) 时显然平凡（\(A,B\)）。考虑

\[ \begin{align*} C=\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}i^{(i,j,k)}&=\prod_{i=1}^{a}i^{\sum_{j=1}^{B}\sum_{k=1}^{C}(i,j,k)}\\ &=\prod_{i=1}^{a}i^{\sum_{d|i}\varphi(d)\lfloor B/d\rfloor\lfloor C/d\rfloor}\\ &=\prod_{d=1}^{a}(\prod_{d|i}i)^{\varphi(d)\lfloor B/d\rfloor\lfloor C/d\rfloor}\\ &=\prod_{d=1}^{a}\Bigl(\lfloor\frac Ad\rfloor!d^{\lfloor\frac Ad\rfloor}\Big)^{\varphi(d)\lfloor B/d\rfloor\lfloor C/d\rfloor} \end{align*} \]

直接对 \(d\) 整除分块即可。

再看第二行（\(D,E,F\)），由于外面是 \(\prod\)，里面的 \((i,j)\) 不能用反演拆开，因此枚举 \((i,j)\)：

\[ \prod_{d=1}d^{\sum_{i,j,k}[(i,j)=d]f(i,j,k)} \]

把指数拿出来，莫反：

\[ \sum_{t}\mu(t)\sum_{i=1}^{\lfloor\frac{A}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{B}{dt}\rfloor}\sum_{k=1}^{C}f(dti,dtj,k) \]

带回原式，再改成枚举 \(dt=x\)：

\[ \begin{align*} \prod_{d=1}d^{\sum_{t}\mu(t)\sum_{i=1}^{\lfloor\frac{A}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{B}{dt}\rfloor}\sum_{k=1}^{C}f(dti,dtj,k)}&=\prod_{x}\biggl(\prod_{d}d^{\mu(x/d)}\bigg)^{\sum_{i=1}^{\lfloor\frac{A}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{B}{x}\rfloor}\sum_{k=1}^{C}f(xi,xj,k)}\\ &=\prod_{x}g(x)^{\sum_{i=1}^{\lfloor\frac{A}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{B}{x}\rfloor}\sum_{k=1}^{C}f(xi,xj,k)} \end{align*} \]

\(g(x)\) 可以 \(O(n\ln n)\) 预处理。\(f(i,j,k)=1,~ijk\) 时直接整除分块（\(D,E\)）。\(f(i,j,k)=\gcd(i,j,k)\) 时，指数上那一坨不好用欧拉反演，因此考虑枚举 \((i,j,k)\)，再表示出 \((i,j)\)：

\[ \begin{align*} F=\prod_{d}\prod_{d|i}\prod_{d|j}(di,dj)^{d\sum_{k}^{C/d}[(i,j,k)=1]}&=\prod_{d}\prod_{i=1}^{A/d}\prod_{j=1}^{B/d}\bigl(d\times\gcd(i,j)\big)^{d\sum_{k}^{C/d}[(i,j,k)=1]}\\ &=\prod_d\prod_t\prod_{i=1}^{A/dt}\prod_{j=1}^{B/dt}\bigl(dt\times\gcd(i,j)\big)^{\mu(t)d\lfloor C/dt\rfloor} \end{align*} \]

把底数拆开：

\[ \biggl(\prod_d\prod_t(dt)^{\mu(t)d\lfloor A/dt\rfloor\lfloor B/dt\rfloor\lfloor C/dt\rfloor}\bigg)\biggl(\prod_d\prod_t\prod_{i=1}^{A/dt}\prod_{j=1}^{B/dt}\gcd(i,j)^{\mu(t)d\lfloor C/dt\rfloor}\bigg) \]

左边那个是

\[ F_1=\prod_{x}x^{(\sum_{dt=x}\mu(t)d)\lfloor A/x\rfloor\lfloor B/x\rfloor\lfloor C/x\rfloor}=\prod_{x}x^{\varphi(x)\lfloor A/x\rfloor\lfloor B/x\rfloor\lfloor C/x\rfloor} \]

直接整除分块。右边那个是

\[ \begin{align*} F_2=\prod_x\prod_{i=1}^{A/x}\prod_{j=1}^{B/x}\gcd(i,j)^{\varphi(x)\lfloor C/x\rfloor}&=\prod_x\prod_d\prod_t\prod_{i=1}^{A/dtx}\prod_{j=1}^{B/dtx}d^{\mu(t)\varphi(x)\lfloor C/x\rfloor}\\ &=\prod_x\prod_yg(y)^{\varphi(x)\lfloor A/xy\rfloor\lfloor B/xy\rfloor\lfloor C/x\rfloor} \end{align*} \]

时间复杂度 \(O(Tn^{3/4}\log n)\)，瓶颈在 \(F_2\)，其它的都是 \(O(T\log n)\) 或者 \(O(T\sqrt n\log n)\)。把能预处理的都预处理了，尽可能把瓶颈处的运算都挪到外面，还是能卡过的。

代码 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include<iostream>
#define int long long
using namespace std;
const int N = 1e5 + 10;

int T, MOD, MOD2;
int a, b, c;

inline int qpow(int a, int b) {
    int res = 1;
    while(b) { (b & 1) && (res = res * a % MOD); a = a * a % MOD; b >>= 1; }
    return res;
}

inline int get_inv(int a) { return qpow(a, MOD - 2); }

int npri[N], pri[N], phi[N], mu[N], pcnt;
int fact[N], ifact[N], inv[N], s[N];

int sphi[N];
int f[N], sf[N], sif[N], si2f[N], sii2f[N], g[N], ig[N], h[N];

inline int A(int a, int b, int c) { return qpow(fact[a], b * c % MOD2); }
inline int B(int a, int b, int c) { return qpow(h[a], s[b] * s[c] % MOD2); }
inline int C(int a, int b, int c) {
    int res = 1;
    for(int l = 1, r; l <= a && l <= b && l <= c; l = r + 1) {
        r = min(a / (a / l), min(b / (b / l), c / (c / l)));
        int x = g[r] * ig[l - 1] % MOD;
        x = qpow(x, a / l);
        int y = qpow(fact[a / l], (sphi[r] - sphi[l - 1]) % MOD2);
        (res *= qpow(x * y % MOD, (b / l) * (c / l) % MOD2)) %= MOD;
    }
    return res;
}
inline int D(int a, int b, int c) {
    int res = 1;
    for(int l = 1, r; l <= a && l <= b; l = r + 1) {
        r = min(a / (a / l), b / (b / l));
        (res *= qpow(sf[r] * sif[l - 1] % MOD, (a / l) * (b / l) * c % MOD2)) %= MOD;
    }
    return res;
}
inline int E(int a, int b, int c) {
    int res = 1;
    for(int l = 1, r; l <= a && l <= b; l = r + 1) {
        r = min(a / (a / l), b / (b / l));
        (res *= qpow(si2f[r] * sii2f[l - 1] % MOD, s[a / l] * s[b / l] % MOD2 * s[c] % MOD2)) %= MOD;
    }
    return res;
}
inline int F1(int a, int b, int c) {
    int res = 1;
    for(int l = 1, r; l <= a && l <= b && l <= c; l = r + 1) {
        r = min(a / (a / l), min(b / (b / l), c / (c / l)));
        int x = g[r] * ig[l - 1] % MOD;
        (res *= qpow(x, (a / l) * (b / l) * (c / l) % MOD2)) %= MOD;
    }
    return res;
}
inline int F2(int a, int b, int c) {
    int res = 1;
    for(int l = 1, r; l <= a && l <= b && l <= c; l = r + 1) {
        r = min(a / (a / l), min(b / (b / l), c / (c / l)));
        int tmp = 1, a1 = a / l, b1 = b / l, c1 = c / l;
        for(int l = 1, r; l <= a1 && l <= b1; l = r + 1) {
            r = min(a1 / (a1 / l), b1 / (b1 / l));
            (tmp *= qpow(sf[r] * sif[l - 1] % MOD, (a1 / l) * (b1 / l) % MOD2)) %= MOD;
        }
        (res *= qpow(tmp, (sphi[r] - sphi[l - 1]) * c1 % MOD2)) %= MOD;
    }
    return res;
}
inline int F(int a, int b, int c) { return F1(a, b, c) * F2(a, b, c) % MOD; }

signed main() {

    cin >> T >> MOD; MOD2 = MOD - 1;

    fact[0] = ifact[0] = inv[1] = 1;
    for(int i = 2; i < N; i++) inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
    for(int i = 1; i < N; i++) fact[i] = fact[i - 1] * i % MOD;
    for(int i = 1; i < N; i++) ifact[i] = ifact[i - 1] * inv[i] % MOD;
    phi[1] = mu[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!npri[i]) { pri[++pcnt] = i; phi[i] = i - 1; mu[i] = -1; }
        for(int j = 1; j <= pcnt; j++) {
            if(i * pri[j] >= N) break;
            npri[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                mu[i * pri[j]] = 0;
                break;
            } else {
                phi[i * pri[j]] = phi[i] * (pri[j] - 1);
                mu[i * pri[j]] = -mu[i];
            }
        }
    }
    for(int i = 1; i < N; i++) sphi[i] = sphi[i - 1] + phi[i];

    for(int i = 1; i < N; i++) s[i] = i * (i + 1) / 2 % MOD2;
    for(int i = 1; i < N; i++) f[i] = 1;
    for(int i = 1; i < N; i++) {
        for(int j = 1; i * j < N; j++) {
            if(mu[j] == 1) (f[i * j] *= i) %= MOD;
            else if(mu[j] == -1) (f[i * j] *= inv[i]) %= MOD;
        }
    }
    sf[0] = sif[0] = si2f[0] = sii2f[0] = 1;
    for(int i = 1; i < N; i++) sf[i] = sf[i - 1] * f[i] % MOD;
    for(int i = 1; i < N; i++) sif[i] = get_inv(sf[i]);
    for(int i = 1; i < N; i++) si2f[i] = si2f[i - 1] * qpow(f[i], i * i % MOD2) % MOD;
    for(int i = 1; i < N; i++) sii2f[i] = get_inv(si2f[i]);

    g[0] = ig[0] = 1;
    for(int i = 1; i < N; i++) g[i] = g[i - 1] * qpow(i, phi[i]) % MOD;
    for(int i = 1; i < N; i++) ig[i] = get_inv(g[i]);

    h[0] = 1;
    for(int i = 1; i < N; i++) h[i] = h[i - 1] * qpow(i, i) % MOD;

    while(T--) {
        cin >> a >> b >> c;
        cout << (A(a, b, c) * A(b, a, c) % MOD * get_inv(D(a, b, c) * D(a, c, b) % MOD) % MOD) << ' ';
        cout << (B(a, b, c) * B(b, a, c) % MOD * get_inv(E(a, b, c) * E(a, c, b) % MOD) % MOD) << ' ';
        cout << (C(a, b, c) * C(b, a, c) % MOD * get_inv(F(a, b, c) * F(a, c, b) % MOD) % MOD) << '\n';
    }

    return 0;
}