投稿
  1. 高性能架构探索首页
  2. C/C++

浅谈莫比乌斯反演

admin • • C/C++ • 阅读 123

狄利克雷卷积

定义两个数论函数的狄利克雷卷积为：

\[f*g=\sum\limits_{d|n}f(d)g(\dfrac{n}{d})
\]

其中单位元\(e(n)=[n=1]\)

可以证明狄利克雷卷积满足交换律及结合律

莫比乌斯函数

对于\(n=\prod\limits_{i=1}^{k}p_i^{c_i}(p_i\in \text{Prime})\)，定义莫比乌斯函数\(\mu(n)\)：

\[\mu(n)=
\begin{cases}
1\quad&(n=1)\\
0&(\exists\, c_i\ge 2)\\
(-1)^k&(\forall i\in[1,k],c_i<2)
\end{cases}
\]

对于莫比乌斯函数，存在一个很重要的性质：

\[\sum\limits_{d|n}\mu(d)=e(n)=[n=1]
\]

证明：

我们考虑枚举\(n\)的所有因子，显然某个质因子的指数\(\ge 2\)时该函数值为\(0\)，对结果不产生贡献，故当我们对于\(n=\prod\limits_{i=1}^{k}p_i^{c_i}\)定义\(x=\prod\limits_{i=1}^kp_i\)时，有：

\[\sum\limits_{d|n}\mu(d)=\sum\limits_{d|x}\mu(d)
\]

考虑\(x\)的所有质因子\(p_1,p_2,p_3\dots p_k\)，枚举因数等价于从中选择任意个数的质因子进行组合，所以显然有：

\[\sum\limits_{d|x}\mu(d)=\sum\limits_{i=0}^{k}\begin{pmatrix}k\\i\end{pmatrix}(-1)^i
\]

反向利用二项式定理：

\[\begin{equation}
\begin{aligned}
\sum\limits_{i=0}^{k}\begin{pmatrix}k\\i\end{pmatrix}(-1)^i&=\sum\limits_{i=0}^{k}\begin{pmatrix}k\\i\end{pmatrix}(-1)^i\times 1^{k-i}\\
&=(-1+1)^k\\
&=0^k
\end{aligned}
\end{equation}
\]

显然原式当且仅当\(k=0\)时等于\(1\)，否则等于\(0\)

当\(k=0\)时有\(n=1\)，故我们证明了：

\[\sum\limits_{d|n}\mu(d)=[n=1]=e(n)
\]

同时该结论也可以写成狄利克雷卷积的形式：

\[e=\mu * 1
\]

注意到莫比乌斯函数是积性函数，所以很好求，线性筛就可以了

inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
        }
    }
}

莫比乌斯反演

对于数论函数\(f(n),g(n)\)，如果有：

\[f(n)=\sum\limits_{d|n}g(d)
\]

则可以推出：

\[g(n)=\sum\limits_{d|n}f(d)\mu(\dfrac{n}{d})
\]

证明：

将等式右边写成狄利克雷卷积的形式：

\[\sum\limits_{d|n}f(d)\mu(\dfrac{n}{d})=f*\mu
\]

同样将条件写成狄利克雷卷积的形式：

\[f(n)=\sum\limits_{d|n}g(d)=g*1
\]

所以：

\[\begin{equation}
\begin{aligned}
\sum\limits_{d|n}f(d)\mu(\dfrac{n}{d})&=f*\mu\\
&=g*1*\mu\\
&=g*(1*\mu)\\
&=g*e\\
&=g
\end{aligned}
\end{equation}
\]

即：

\[g(n)=\sum\limits_{d|n}f(d)\mu(\dfrac{n}{d})
\]

于是我们证明了莫比乌斯反演定理

易证上式也有个等价形式：

\[g(n)=\sum\limits_{d|n}\mu(d)f(\dfrac{n}{d})
\]

同时我们还有第二种反演形式：

若

\[g(n)=\sum\limits_{n|d}f(d)
\]

则同样可以推出

\[f(n)=\sum\limits_{n|d}\mu(\dfrac{d}{n})g(d)
\]

证明：

设\(d=kn\)：

\[\begin{equation}
\begin{aligned}
\sum\limits_{n|d}\mu(\dfrac{d}{n})g(d)&=\sum\limits_{k\ge 1}\mu(k)g(nk)\\
&=\sum\limits_{k\ge 1}\mu(k)\sum\limits_{nk|d}f(d)\\
&=\sum\limits_{n|d}f(d)\sum\limits_{k|\frac{d}{n}}\mu(k)\\
&=\sum\limits_{n|d}f(d)[\frac{d}{n}=1]\\
&=f(n)
\end{aligned}
\end{equation}
\]

整除分块

求\(
\sum\limits_{i=1}^{n}\lfloor\dfrac{n}{i}\rfloor
\)

显然暴力计算复杂度为\(O(n)\)，但我们注意到很多值其实是一样的，所以我们可以把值相同的统一计算，分块处理：

for (int l = 1, r; l <= n; l = r + 1) {
    r = (n / (n / l));
    ans += (r - l + 1) * (n / l);
}

复杂度\(O(\sqrt n)\)

例题

[POI2007]ZAP-Queries

\(T\)组询问，给出\(a,b,d\)，求\(\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=d]\)
\(1\le T\le 5\times 10^4\)
\(1\le d\le a,b\le 5\times 10^4\)

Solution 1

我们设：

\[f(k)=\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=k]\\
g(n)=\sum\limits_{n|k}f(k)
\]

可以发现\(f(d)\)就是我们要求的答案

显然根据莫比乌斯反演定理可以得到：

\[f(n)=\sum\limits_{n|d}g(d)\mu(\dfrac{d}{n})
\]

注意到：

\[\begin{equation}
\begin{aligned}
g(n)&=\sum\limits_{n|k}f(k)\\
&=\sum\limits_{n|k}\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=k]\\
&=\lfloor\dfrac{a}{n}\rfloor\lfloor\dfrac{b}{n}\rfloor
\end{aligned}
\end{equation}
\]

所以：

\[\begin{equation}
\begin{aligned}
f(n)&=\sum\limits_{n|d}g(d)\mu(\dfrac{d}{n})\\
&=\sum\limits_{n|d}\mu(\dfrac{d}{n})\lfloor\dfrac{a}{d}\rfloor\lfloor\dfrac{b}{d}\rfloor
\end{aligned}
\end{equation}
\]

设\(\dfrac{d}{n}=t\)，枚举\(t\)：

\[f(n)=\sum\limits_{t=1}^{min(\lfloor\frac{a}{n}\rfloor,\lfloor\frac{b}{n}\rfloor)}\mu(t)\lfloor\dfrac{a}{nt}\rfloor\lfloor\dfrac{b}{nt}\rfloor
\]

前面的\(\mu(t)\)可以前缀和预处理，后面用整除分块即可

总复杂度\(O(T\sqrt{min(a,b)})\)

#include <bits/stdc++.h>
using namespace std;
#define MAXN 50005

int mu[MAXN], sum[MAXN];
bool isprime[MAXN];
vector<int> prime;
inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
        }
    }

    for (i = 1; i <= n; ++i) sum[i] = sum[i - 1] + mu[i];
}

int T, a, b, d;
signed main() {
    init(MAXN);
    scanf("%d", &T);
    while (T--) {
        scanf("%d%d%d", &a, &b, &d);
        int n = min(a / d, b / d);
        int l, r;

        int ans = 0;
        for (l = 1; l <= n; l = r + 1) {
            r = min((a / d) / ((a / d) / l), (b / d) / ((b / d) / l));
            ans += ((a / d) / l) * ((b / d) / l) * (sum[r] - sum[l - 1]);
        }
        printf("%d\n", ans);
    }

    return 0;
}

Solution 2

可以发现上述利用莫比乌斯反演定理构造函数进行求解的过程不太好想，其实我们还可以利用莫比乌斯函数的性质顺向思考求解

我们要求：

\[\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=x]
\]

（这里用\(x\)替换\(d\)避免变量名重复）

约去\(x\)：

\[\sum\limits_{i=1}^{\lfloor\frac{a}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{x}\rfloor}[gcd(i,j)=1]
\]

莫比乌斯函数有个性质：

\[\sum\limits_{d|n}\mu(d)=[n=1]
\]

所以：

\[[gcd(i,j)=1]=\sum\limits_{d|gcd(i,j)}\mu(d)
\]

带回原式：

\[\sum\limits_{i=1}^{\lfloor\frac{a}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{x}\rfloor}\sum\limits_{d|gcd(i,j)}\mu(d)
\]

枚举\(d\)：

\[\sum\limits_{d=1}^{min(\lfloor\frac{a}{x}\rfloor,\lfloor\frac{b}{x}\rfloor)}\sum\limits_{i=1}^{\lfloor\frac{a}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{x}\rfloor}[d|gcd(i,j)]
\]

可以发现当且仅当\(i,j\)均为\(d\)的倍数时\(d|gcd(i,j)\)，同时在\(1-n\)中\(d\)的倍数有\(\lfloor\dfrac{n}{d}\rfloor\)个，所以原式等价于：

\[\sum\limits_{d=1}^{min(\lfloor\frac{a}{x}\rfloor,\lfloor\frac{b}{x}\rfloor)}\mu(d)\lfloor\dfrac{a}{dx}\rfloor\lfloor\dfrac{b}{dx}\rfloor
\]

结果与\(\text{Solution 1}\)一致

[HAOI2011]Problem b

\(T\)组询问，给出\(a,b,c,d,k\)，求\(\sum\limits_{i=a}^b\sum\limits_{j=c}^d[gcd(i,j)=k]\)
\(1\le T\le 5\times 10^4\)
\(1\le k\le a,c\le b,d\le 5\times 10^4\)

做法和上一题几乎完全一致，容斥一下即可

令：

\[solve(a,b,c,d)=\sum\limits_{i=a}^b\sum\limits_{j=c}^d[gcd(i,j)=k]
\]

显然有：

\[ans=solve(1,b,1,d)-solve(1,b,1,c-1)-solve(1,a-1,1,d)+solve(1,a-1,1,c-1)
\]

#include <bits/stdc++.h>
using namespace std;
#define MAXN 50005

int mu[MAXN], sum[MAXN];
bool isprime[MAXN];
vector<int> prime;
inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
        }
    }

    for (i = 1; i <= n; ++i) sum[i] = sum[i - 1] + mu[i];
}

inline int solve(int a, int b, int d) {
    int n = min(a / d, b / d);
    int l, r;

    int ans = 0;
    for (l = 1; l <= n; l = r + 1) {
        r = min((a / d) / ((a / d) / l), (b / d) / ((b / d) / l));
        ans += ((a / d) / l) * ((b / d) / l) * (sum[r] - sum[l - 1]);
    }

    return ans;
}

int T, a, b, c, d, k;
signed main() {
    init(MAXN);
    scanf("%d", &T);
    while (T--) {
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        printf("%d\n", solve(b, d, k) - solve(b, c - 1, k) -
                           solve(a - 1, d, k) + solve(a - 1, c - 1, k));
    }

    return 0;
}

P2257 YY的GCD

给定 \(N, M\)，求 \(1 \le x \le N,1 \leq y \leq M\)且 \(\gcd(x, y)\) 为质数的 \((x,y)\) 有多少对。

简化题意就是要求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m [\gcd(i,j)\in \text{Prime}]
\]

Solution 1

枚举质数：

\[\sum\limits_{p\in \text{Prime} }\sum\limits_{i=1}^n\sum\limits_{j=1}^m [\gcd(i,j)=p]
\]

提出\(p\)：

\[\sum\limits_{p\in \text{Prime} }\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor
} [\gcd(i,j)=1]
\]

替换\([\gcd(i,j)=1]\)：

\[\sum\limits_{p\in \text{Prime} }\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor
} \sum\limits_{d|\gcd(i,j)}\mu(d)
\]

枚举\(d\)：

\[\sum\limits_{p\in \text{Prime}}\sum\limits_{d=1}^{min(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor
}[d|\gcd(i,j)]
\]

显然由于之前的经验上式等于：

\[\sum\limits_{p\in \text{Prime}}\sum\limits_{d=1}^{min(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor)}\mu(d)\lfloor\dfrac{n}{pd}\rfloor\lfloor\dfrac{m}{pd}\rfloor
\]

令\(k=pd\)，枚举\(k\)：

\[\sum\limits_{k=1}^{min(n,m)}\sum\limits_{p\in \text{Prime},p|k}\mu(\dfrac{k}{p})\lfloor\dfrac{n}{k}\rfloor\lfloor\dfrac{m}{k}\rfloor
\]

令\(f(x)=\sum\limits_{p\in \text{Prime},p|x}\mu(\dfrac{x}{p})\)，则原式等于：

\[\sum\limits_{k=1}^{min(n,m)}f(k)\lfloor\dfrac{n}{k}\rfloor\lfloor\dfrac{m}{k}\rfloor
\]

注意到如果我们求出了\(f\)的前缀和，即可像之前一样通过整除分块快速求解

考虑\(f(x)\)的取值：

令\(x=i\times y\)，其中\(y\)是\(x\)的最小质因子

\(1.\)若\(x\in \text{Prime}\)，显然\(f(x)=1\)

否则，考虑\(\mu(\dfrac{x}{p})=\mu(\dfrac{i\times y}{p})=\mu(i)\times \mu(\dfrac{y}{p})=\mu(\dfrac{i}{p})\times \mu(y)\)

\(2.\)若\(i\mod y=0\)

\(\quad 2.1\)若\(\mu(i)\not =0\)，则当切仅当枚举的\(p=y\)时\(f(x)\not =0\)，所以有\(f(x)=\mu(i)\)

\(\quad 2.2\)若\(\mu(i) =0\)，则无论何时\(f(x) =0\)，所以有\(f(x)=\mu(i)\)

\(3.\)若\(i\mod y\not =0\)，则\(\mu(\dfrac{x}{p})=\mu(\dfrac{i}{p})\times \mu(y)=-\mu(\dfrac{i}{p})\)，即\(f(x)\)中包含\(-f(i)\)，同时和\(f(i)\)相比\(f(x)\)多出的一项就是当\(p=y\)时，\(f(x)=\mu(\dfrac{i\times y}{y})=\mu(i)\)，所以\(f(x)=-f(i)+\mu(i)\)

于是我们得到了\(f(x)\)的方程：

\[f(x)=\begin{cases}1& x\in \text{Prime}\\
\mu(i)&i\mod y=0\\
-f(i)+\mu(i)&i\mod y\not =0
\end{cases}
\]

可以通过线性筛求出\(f(x)\)的前缀和

#include <bits/stdc++.h>
using namespace std;
#define MAXN 10000005
#define int long long

int mu[MAXN], sum[MAXN], f[MAXN];
bool isprime[MAXN];
vector<int> prime;
inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1, f[i] = 1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                f[i * prime[j]] = mu[i];
                mu[i * prime[j]] = 0;
                break;
            } else {
                f[i * prime[j]] = -f[i] + mu[i];
                mu[i * prime[j]] = mu[i] * mu[prime[j]];
            }
        }
        f[i] += f[i - 1];
    }
}

int T, a, b, d;
signed main() {
    init(MAXN);
    scanf("%lld", &T);
    while (T--) {
        scanf("%lld%lld", &a, &b);
        int n = min(a, b);
        int l, r;

        int ans = 0;
        for (l = 1; l <= n; l = r + 1) {
            r = min(a / (a / l), b / (b / l));
            ans += (a / l) * (b / l) * (f[r] - f[l - 1]);
        }
        printf("%lld\n", ans);
    }

    return 0;
}

在你谷因为用了vector所以上述代码吸氧才能过

Solution 2

我们以前做过这道题：

\[\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[gcd(i,j)=1]
\]

当时我们利用了\(\sum\limits_{d|n}\mu(d)=[n=1]\)这个性质

现在我们要求：

\[\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[gcd(i,j)\in \text{Prime}]
\]

那我们也可以自己构造一个函数满足\(\sum\limits_{d|n}f(d)=[n\in \text{Prime}]\)

我们先给\(f\)赋初值：

\[f(x)=[x\in \text{Prime}]
\]

考虑如何让\(\sum\limits_{d|n}f(d)=[n\in \text{Prime}]\)：

for (i = 2; i <= MAXN; ++i) {
        for (j = i << 1; j <= MAXN; j += i) {
            f[j] -= f[i];
        }
    }

不难发现如此处理即可得到合法的\(f\)

[SDOI2015]约数个数和

设\(d(x)\)为\(x\)的约数个数，\(T\)组询问，求\(\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}d(xy)\)
\(1\le T,n,m\le 5\times 10^4\)

首先有一个结论：

\[d(xy)=\sum\limits_{i|x}\sum\limits_{j|y}[\gcd(i,j)=1]
\]

证明：

设\(x=\prod\limits_{i=1}^{k}p_i^{a_i},y=\prod\limits_{i=1}^{k}p_i^{b_i}\)（不存在该质因子则指数为\(0\)），显然有：

\[d(xy)=\prod\limits_{i=1}^{k}(a_i+b_i+1)
\]

对于某个质因子\(p\)，考虑等式右边如果要满足\(\gcd(i,j)=1\)，显然不能同时拥有\(p\)因子。令\(x\)中含有\(p^a\)，\(y\)中含有\(p^b\)：

\(1.\)如果\(i\)中不含\(p\)因子，则\(j\)中可选质因子\(p\)次数为\([0,b]\)

\(2.\)如果\(j\)中不含\(p\)因子，则\(i\)中可选质因子\(p\)次数为\([0,a]\)

二者相加再减去\(i,j\)都不含\(p\)的一种情况，结果数为\(a+b+1\)

考虑乘法原理，对于每个\(p_i\)结果数都是\(a_i+b_i+1\)，所以：

\[\sum\limits_{i|x}\sum\limits_{j|y}[\gcd(i,j)=1]=\prod\limits_{i=1}^{k}(a_i+b_i+1)
\]

证毕

有了上述结论，我们就可以很轻松的化简原式：

\[\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}\sum\limits_{i|x}\sum\limits_{j|y}[\gcd(i,j)=1]
\]

枚举\(i,j\)：

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor[\gcd(i,j)=1]
\]

设：

\[f(x)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor[\gcd(i,j)=x]\\
g(x)=\sum\limits_{x|d}f(d)
\]

莫比乌斯反演：

\[f(n)=\sum\limits_{n|d}\mu(\dfrac{d}{n})g(d)
\]

注意到：

\[g(x)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor[x|\gcd(i,j)]
\]

约去\(x\)：

\[g(x)=\sum\limits_{i=1}^{\frac{n}{x}}\sum\limits_{j=1}^{\frac{m}{x}}\lfloor\dfrac{n}{ix}\rfloor\lfloor\dfrac{m}{jx}\rfloor
\]

显然答案为：

\[f(1)=\sum\limits_{1|d}\mu(\dfrac{d}{1})g(d)=\sum\limits_{i=1}^n\mu(i)g(i)
\]

可以整除分块预处理出\(s(n)=\sum\limits_{i=1}^n\lfloor\dfrac{n}{i}\rfloor\)，从而快速求出\(g(i)\)

#include <bits/stdc++.h>
using namespace std;
#define MAXN 50005
#define int long long

int mu[MAXN], sum[MAXN], s[MAXN];
bool isprime[MAXN];
vector<int> prime;
inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
        }
    }

    for (i = 1; i <= n; ++i) sum[i] = sum[i - 1] + mu[i];
    for (i = 1; i <= n; ++i) {
        int ans = 0;
        int l, r;
        for (l = 1; l <= i; l = r + 1) {
            r = i / (i / l);
            ans += (i / l) * (r - l + 1);
        }
        s[i] = ans;
    }
}

int T, a, b;
signed main() {
    init(MAXN);
    scanf("%lld", &T);
    while (T--) {
        scanf("%lld%lld", &a, &b);
        int n = min(a, b);
        int l, r;

        int ans = 0;
        for (l = 1; l <= n; l = r + 1) {
            r = min(a / (a / l), b / (b / l));
            ans += s[a / l] * s[b / l] * (sum[r] - sum[l - 1]);
        }
        printf("%lld\n", ans);
    }

    return 0;
}

[国家集训队]Crash的数字表格

给出\(n,m\)，求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)\)
\(1\le n,m\le 10^7\)
答案对\(20101009\)取模

\[\begin{equation}
\begin{aligned}
\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)&=\sum\limits_{i=1}^n\sum\limits_{j=1}^m \dfrac{ij}{\gcd(i,j)}\\
&=\sum\limits_{x=1}^{min(n,m)}\sum\limits_{i=1}^n\sum\limits_{j=1}^m [\gcd(i,j)=x]\dfrac{ij}{x}\\
&=\sum\limits_{x=1}^{min(n,m)}x\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor} [\gcd(i,j)=1]ij\\
&=\sum\limits_{x=1}^{min(n,m)}x\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor} \sum\limits_{d|gcd(i,j)}\mu(d)ij\\
&=\sum\limits_{x=1}^{min(n,m)}x\sum\limits_{d=1}^{min(\lfloor\frac{n}{x}\rfloor,\lfloor\frac{m}{x}\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor}ij[d|\gcd(i,j)]
\end{aligned}
\end{equation}
\]

令\(i=ud,j=vd\)，则显然有\(d|\gcd(i,j)\)

所以：

\[\begin{equation}
\begin{aligned}
\text{原式}&=\sum\limits_{x=1}^{min(n,m)}x\sum\limits_{d=1}^{min(\lfloor\frac{n}{x}\rfloor,\lfloor\frac{m}{x}\rfloor)}\mu(d)\sum\limits_{ud=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{vd=1}^{\lfloor\frac{m}{x}\rfloor}d^2uv\\
&=\sum\limits_{x=1}^{min(n,m)}x\sum\limits_{d=1}^{min(\lfloor\frac{n}{x}\rfloor,\lfloor\frac{m}{x}\rfloor)}d^2\mu(d)\sum\limits_{u=1}^{\lfloor\frac{n}{dx}\rfloor}u\sum\limits_{v=1}^{\lfloor\frac{m}{dx}\rfloor}v
\end{aligned}
\end{equation}
\]

最后两项显然是等差数列求和，可以用整除分块优化

#include <bits/stdc++.h>
using namespace std;
#define MAXN 10000005
#define mod 20101009
#define int long long

int mu[MAXN], sum[MAXN], s[MAXN];
bool isprime[MAXN];
vector<int> prime;
inline void init(int n) {
    int i, j;
    mu[1] = 1;
    for (i = 1; i <= n; ++i) isprime[i] = 1;
    for (i = 2; i <= n; ++i) {
        if (isprime[i]) prime.push_back(i), mu[i] = -1;
        for (j = 0; j < prime.size() && i * prime[j] <= n; ++j) {
            isprime[i * prime[j]] = 0;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
        }
    }

    for (i = 1; i <= n; ++i)
        sum[i] = (sum[i - 1] + mu[i] * 1ll * i % mod * 1ll * i % mod) % mod;
}

int a, b;
signed main() {
    scanf("%lld%lld", &a, &b);
    init(min(a, b));
    int n = min(a, b);
    int l, r, i, j;

    int inv = (mod + 1ll) / 2ll;

    int ans = 0, tsum = 0;
    for (i = 1; i <= n; ++i) {
        int x = a / i, y = b / i;
        int nn = min(x, y);
        tsum = 0;
        for (l = 1; l <= nn; l = r + 1) {
            r = min(x / (x / l), y / (y / l));
            tsum =
                (tsum +
                 (sum[r] - sum[l - 1]) % mod * 1ll *
                     (((1 + (x / l)) % mod * 1ll * (x / l) % mod) * inv % mod) %
                     mod * 1ll *
                     (((1 + (y / l)) % mod * 1ll * (y / l) % mod) * inv % mod) %
                     mod) %
                mod;
        }
        ans = (ans + tsum * i % mod) % mod;
    }
    printf("%lld\n", (ans % mod + mod) % mod);

    return 0;
}

原文链接: https://www.cnblogs.com/luisvacson/p/15806700.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍

    浅谈莫比乌斯反演

原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/185850

非原创文章文中已经注明原地址，如有侵权，联系删除

关注公众号【高性能架构探索】，第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处！

赞 (0)
adminadmin
0 0
上一篇 2023年2月12日 上午11:01
下一篇 2023年2月12日 上午11:01

相关推荐