浅谈莫比乌斯反演（c++）

前置知识

一些约定

•[A]：当且仅当表达式 A 为真时 [A] = 1，否则 [A] = 0.
•gcd(x, y) 或 (x, y)：表示 x 和 y 的最大公因数.
•lcm(x, y)：表示 x 和 y 的最小公倍数.
•d | n：表示 d 整除 n.
•d ∤ n：表示 d 不整除 n.
•$\left[\begin{array}{c}n\\ d\end{array}\right]$：下取整符号.
• 正体、粗体或希腊字母，如 1, F, σ, d 等：表示数论函数.
• 其他字母如 i, j, d, g, T 等：表示数字

数论分块

问题 1：给定正整数 n，求$n\sum _{i=1}^n\left[\begin{array}{c}n\\ i\end{array}\right]$?
不难发现 $\left[\begin{array}{c}n\\ i\end{array}\right]$的值呈块状分布，每一块内$\left[\begin{array}{c}n\\ i\end{array}\right]$的值相同，且每
个块最后位置均为$\frac{n}{\left[\begin{array}{c}n\\ i\end{array}\right]}$
，所以可以计算，时间复杂度 O(√n)。
求 1 到 n 的 µ 值，可以使用线性筛做到 O(n)。

mu[1] = 1;
for(int i = 2; i <= n; i++) {
if(!not_prime[i]) prime[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i*prime[j] <= n; j++) {
not_prime[i*prime[j]] = true;
if(i%prime[j] == 0) break;
mu[i*prime[j]] = -mu[i];
}
}

也可以打包带走连着 φ 一起筛：

mu[1] = phi[1] = 1;
for(int i = 2; i <= n; i++) {
if(!not_prime[i]) prime[++tot] = i,
mu[i] = -1, phi[i] = i-1;
for(int j = 1; j <= tot && i*prime[j] <= n; j++) {
not_prime[i*prime[j]] = true;
if(i%prime[j] == 0) {
phi[i*prime[j]] = phi[i] * prime[j];
break;
}
mu[i*prime[j]] = -mu[i];
phi[i*prime[j]] = phi[i] * phi[prime[j]];
}
}

狄利克雷卷积

卷积、旋积或褶积(英语：Convolution)是通过两个函数f和g生成第三个函数的一种数学运算，其本质是一种特殊的积分变换，表征函数f与g经过翻转和平移的重叠部分函数值乘积对重叠长度的积分。（源于：百度百科）
上面是百度百科给的定义，第一次接触的时候也是一头雾水，就看懂了一句：卷积是一种数学运算，单着一句就够用了
卷积就是一种数学运算

定义

若数论函数 F, G, H 满足：$F(n)=\sum _{d|n}G(d)H(\frac{n}{d})$
则称 F 是 G 和 H 的狄利克雷卷积，记作 F = G ∗ H

一些常见的狄利克雷卷积

•φ ∗ 1 = id：即 $\sum _{d|n}\phi (d)=n.$
• 对于任意数论函数 F，都有 F ∗ ϵ = F.
•$i{d}^k\ast 1={\sigma }_k.$
•µ ∗ 1 = ϵ
也就是说，µ 实际上是 1 的逆元。

莫比乌斯反演

莫比乌斯函数的性质/莫比乌斯变换

因为 µ 是 1 的逆元，对于任意数论函数 F 有：
F ∗ 1 = G ⇐⇒ G ∗ µ = F
也就是说：
•id ∗ µ = φ.
•${\sigma }_k\ast \mu =i{d}^k.$
•µ ∗ 1 = ϵ.
然后……就结束了
莫比乌斯反演的精髓：数学推导

例题讲解

公约数的和

题目背景

有一天，TIBBAR 和 LXL 比赛谁先算出 $1\sim n$ 这 $n$ 个数中每任意两个不同的数的最大公约数的和。LXL 还在敲一个复杂而冗长的程序，争取能在 $100s$ 内出解。而 TIBBAR 则直接想 $1s$ 秒过而获得完胜，请你帮他完成这个任务。

题目描述

给定 $n$，求
$$\sum _{i=1}^n\sum _{j=i+1}^n\gcd (i,j)$$

其中 $\gcd (i,j)$ 表示 $i$ 和 $j$ 的最大公约数。

输入格式

输入只有一行一个整数，表示 $n$。

输出格式

输出一行一个整数表示答案。

样例 #1

样例输入 #1

10

样例输出 #1

67

提示

数据规模与约定

• 对于 $40\%$ 的数据，保证 $n\le 2\times 10^3$。
• 对于 $100\%$ 的数据，保证 $2\le n\le 2\times 10^6$。

思路

AC代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2e6+5;
ll p[N],n;ll phi[N];
bool notp[N];
void seive(ll n){
    phi[1] = 1;
    for(ll i=2;i<=n;++i) {
        if(!notp[i]) p[++p[0]] = i, phi[i] = i-1;
        for(ll j=1;j<=p[0] && i*p[j]<=n;++j) {
            notp[i*p[j]] = 1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]] = phi[i]*(p[j]-1);
        }
    }
    for(ll i=1;i<=n;++i) phi[i] += phi[i-1];
}
inline ll cal(ll n,ll m){
    ll ans = 0,r;
    for(ll i=1;i<=n;i=r+1) {
        r = min(n/(n/i), m/(m/i));
        ans += (phi[r]-phi[i-1]) * (n/i) * (m/i);
    }
    return ans;
}
int main() {
    cin>>n;
    seive(n);
    cout<<(cal(n,n)-n*(n+1)/2)/2<<endl;
}

[HAOI2011] Problem b

题目描述

对于给出的 $n$ 个询问，每次求有多少个数对 $(x,y)$，满足 $a\le x\le b$，$c\le y\le d$，且 $\gcd (x,y)=k$，$\gcd (x,y)$ 函数为 $x$ 和 $y$ 的最大公约数。

输入格式

第一行一个整数 $n$，接下来 $n$ 行每行五个整数，分别表示 $a,b,c,d,k$。

输出格式

共 $n$ 行，每行一个整数表示满足要求的数对 $(x,y)$ 的个数。

样例 #1

样例输入 #1

2
2 5 1 5 1
1 5 1 5 2

样例输出 #1

14
3

提示

对于 $100\%$ 的数据满足：$1\le n,k\le 5\times 10^4$，$1\le a\le b\le 5\times 10^4$，$1\le c\le d\le 5\times 10^4$。

思路

AC代码

#include<bits/stdc++.h>
#define N 60010
using namespace std;
bool vis[N];
int prim[N],mu[N],sum[N],cnt,k;
inline void get_mu(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]){mu[i]=-1;prim[++cnt]=i;}
        for(int j=1;j<=cnt&&i*prim[j]<=n;j++)
        {
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)break;
            else mu[i*prim[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++)sum[i]=sum[i-1]+mu[i];
}
inline long long calc(int a,int b)
{
    static int max_rep;
    static long long ans;
    max_rep=min(a,b);ans=0;
    for(int l=1,r;l<=max_rep;l=r+1)
    {
        r=min(a/(a/l),b/(b/l));
        ans+=(1ll*a/(1ll*l*k))*(1ll*b/(1ll*l*k))*(sum[r]-sum[l-1]);
    }
    return ans;
}
int main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
    int t;
    cin>>t;
    get_mu(50000);
    while(t--)
    {
        int a,b,c,d;
        cin>>a>>b>>c>>d>>k;
        cout<<calc(b,d)-calc(b,c-1)-calc(a-1,d)+calc(a-1,c-1)<<endl;
    }
    return 0;
}

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

题目描述

今天的数学课上，Crash 小朋友学习了最小公倍数（Least Common Multiple）。对于两个正整数 $a$ 和 $b$，$\text{lcm}(a,b)$ 表示能同时被 $a$ 和 $b$ 整除的最小正整数。例如，$\text{lcm}(6,8)=24$。

回到家后，Crash 还在想着课上学的东西，为了研究最小公倍数，他画了一张 $ n \times m$ 的表格。每个格子里写了一个数字，其中第 $i$ 行第 $j$ 列的那个格子里写着数为 $\text{lcm}(i,j)$。

看着这个表格，Crash 想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当 $n$ 和 $m$ 很大时，Crash 就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash 只想知道表格里所有数的和对 $20101009$ 取模后的值。

输入格式

输入包含一行两个整数，分别表示 $n$ 和 $m$。

输出格式

输出一个正整数，表示表格中所有数的和对 $20101009$ 取模后的值。

样例 #1

样例输入 #1

4 5

样例输出 #1

122

提示

样例输入输出 1 解释

该表格为：

$1$	$2$	$3$	$4$	$5$
$2$	$2$	$6$	$4$	$10$
$3$	$6$	$3$	$12$	$15$
$4$	$4$	$12$	$4$	$20$

数据规模与约定

• 对于 $30\%$ 的数据，保证 $n,m\le 10^3$。
• 对于 $70\%$ 的数据，保证 $n,m\le 10^5$。
• 对于 $100\%$ 的数据，保证 $1\le n,m\le 10^7$。

思路

因此前半部分可以预处理出前缀和，这样 f(n, m) 就能数论分块
求了. 因为原问题也可以数论分块，所以求解的复杂度为 $O(n^{\frac{3}{4}})$，证明过程略.
但是由于需要预处理 $\mu (d)d^2$ 的前缀和，瓶颈在线性筛和预处理，
所以总复杂度是线性的 O(n)

AC代码

#include <cstdio>
#include <algorithm>
using namespace std;
const int N=1e7+5,mod=20101009;
int n,m,mu[N],p[N/10],sum[N];
bool flg[N];
inline void init() {
    mu[1]=1;
    int tot=0,k=min(n,m);
    for(int i=2;i<=k;++i) {
        if(!flg[i]) p[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*p[j]<=k;++j) {
            flg[i*p[j]]=1;
            if(i%p[j]==0) {mu[i*p[j]]=0;break;}
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=k;++i) sum[i]=(sum[i-1]+1LL*i*i%mod*(mu[i]+mod))%mod;
}
inline int Sum(int x,int y) {
    return (1LL*x*(x+1)/2%mod)*(1LL*y*(y+1)/2%mod)%mod;
}
inline int func(int x,int y) {
    int res=0;
    for(int i=1,j;i<=min(x,y);i=j+1) {
        j=min(x/(x/i),y/(y/i));
        res=(res+1LL*(sum[j]-sum[i-1]+mod)*Sum(x/i,y/i)%mod)%mod;
    }
    return res;
}
inline int solve(int x,int y) {
    int res=0;
    for(int i=1,j;i<=min(x,y);i=j+1) {
        j=min(x/(x/i),y/(y/i));
        res=(res+1LL*(j-i+1)*(i+j)/2%mod*func(x/i,y/i)%mod)%mod;
    }
    return res;
}
int main() {
    scanf("%d%d",&n,&m);
    init();
    printf("%d\n",solve(n,m));
}

[SDOI2015] 约数个数和

题目描述

设 $d(x)$ 为 $x$ 的约数个数，给定 $n,m$，求
$$\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

输入格式

输入文件包含多组测试数据。
第一行，一个整数 $T$，表示测试数据的组数。
接下来的 $T$ 行，每行两个整数 $n,m$。

输出格式

$T$ 行，每行一个整数，表示你所求的答案。

样例 #1

样例输入 #1

2
7 4
5 6

样例输出 #1

110
121

提示

【数据范围】
对于 $100\%$ 的数据，$1\le T,n,m\le 50000$。

思路

AC代码

#include <bits/stdc++.h>
using namespace std;
const long long N=5e4+5;
long long tot,mu[N],p[N];
long long s[N];
bool flg[N];
void init() {
    mu[1]=1;
    for(long long i=2;i<=5e4;++i) {
        if(!flg[i]) p[++tot]=i,mu[i]=-1;
        for(long long j=1;j<=tot&&i*p[j]<=5e4;++j) {
            flg[i*p[j]]=1;
            if(i%p[j]==0) {
                mu[i*p[j]]=0;
                break;
            } else {
                mu[i*p[j]]=-mu[i];
            }
        }
    }
    for(long long i=1;i<=5e4;++i) mu[i]+=mu[i-1];
    for(long long x=1;x<=5e4;++x) {
        long long res=0;
        for(long long i=1,j;i<=x;i=j+1) j=x/(x/i),res+=(j-i+1)*(x/i);
        s[x]=res;
    }
    return ;
}
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	init();
	long long T;
	cin>>T;
	while(T--){
		long long n,m;
		cin>>n>>m;
		if(n>m)swap(n,m);
		long long ans=0;
		for(long long i=1,j;i<=n;i=j+1){
			j=min(n/(n/i),m/(m/i));
			ans+=(mu[j]-mu[i-1])*s[n/i]*s[m/i];
		}
		cout<<ans<<endl;
	}
}

[SDOI2008] 仪仗队

题目描述

作为体育委员，C 君负责这次运动会仪仗队的训练。仪仗队是由学生组成的 $N\times N$ 的方阵，为了保证队伍在行进中整齐划一，C 君会跟在仪仗队的左后方，根据其视线所及的学生人数来判断队伍是否整齐（如下图）。

现在，C 君希望你告诉他队伍整齐时能看到的学生人数。

输入格式

一行，一个正整数 $N$。

输出格式

输出一行一个数，即 C 君应看到的学生人数。

样例 #1

样例输入 #1

4

样例输出 #1

9

提示

对于 $100\%$ 的数据，$1\le N\le 40000$。

思路

AC代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<queue>
#include<algorithm>
using namespace std;
typedef long long ll;
  
const int maxn=50000;
  
int vis[maxn];
int prime[maxn];
int phi[maxn];
int sum[maxn];
  
int main()
{
    phi[1]=1;
    sum[1]=1;
    int k=-1;
    for(int i=2;i<=40000;i++)
    {
        if(!vis[i])
        {
            prime[++k]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<=k&&i*prime[j]<=40000;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
        sum[i]=sum[i-1]+phi[i];
    }
    int n;
    scanf("%d",&n);
    if(n==1)
        printf("0\n");
    else
        printf("%d\n",2*sum[n-1]+1);
    return 0;
}

Product

题目背景

$\mathrm{CYJian}$：“听说$gcd$和$\sum$套起来比较好玩？？那我就…”

题目描述

$\mathrm{CYJian}$最近闲的玩起了$gcd$。。他想到了一个非常简单而有意思的式子：

$$\prod _{i=1}^N\prod _{j=1}^N\frac{lcm(i,j)}{gcd(i,j)}(\mathrm{mod}104857601)$$

$\mathrm{CYJian}$已经算出来这个式子的值了。现在请你帮他验算一下吧。$\mathrm{CYJian}$只给你$0.2s$的时间哦。

2024.5.11 upd: 放宽时空限制。

输入格式

一行一个正整数$N$。

输出格式

一行一个正整数，表示答案模$104857601$的值。

样例 #1

样例输入 #1

5

样例输出 #1

585494

提示

样例解释：

$\frac{lcm}{gcd}$	1	2	3	4	5
1	1	2	3	4	5
2	2	1	6	2	10
3	3	6	1	12	15
4	4	2	12	1	20
5	5	10	15	20	1

对于$30\%$的数据：$1\le N\le 5000$

对于$100\%$的数据：$1\le N\le 1000000$

思路

AC代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=1000005,mod=104857601;
ll n,cnt,ans1=1,prim[N],ans2=1,pai[N];
bool vis[N];
inline ll power( ll a,ll b){
	 ll r=1;
    while(b)
    {
        if(b&1ll) r=1ll*r*a%mod;
        b>>=1ll;
        a=1ll*a*a%mod;
    }
    return r;
}
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	cin>>n;
    pai[1]=1;
    for(ll i=2;i<=n;++i)
    {
    	ans1=ans1*i%mod;
        if(!vis[i]) prim[++cnt]=i,pai[i]=i-1;
        for(ll j=1;j<=cnt;++j)
        {
            if(prim[j]*i>n) break;
            vis[prim[j]*i]=1;
            if(i%prim[j]==0) {pai[i*prim[j]]=pai[i]*prim[j];break;}
            pai[i*prim[j]]=pai[prim[j]]*pai[i];
        }
    }
    for(ll i=1;i<=n;++i) pai[i]=pai[i]*2+pai[i-1]%(mod-1);
    ans1=power(ans1,2*n);
    for(ll i=2;i<=n;++i) ans2=ans2*power(i,pai[n/i]-1)%mod;
    cout<<(ans1*power(1ll*ans2*ans2%mod,mod-2))%mod<<endl;
    return 0;
}

这是我的第十五篇文章，如有纰漏也请各位大佬指正

辛苦创作不易，还望看官点赞收藏打赏，后续还会更新新的内容。