莫比乌斯反演
莫比乌斯反演
引例：[HAOI2011]Problem b
- 题目描述：
对于给出的T个询问，每次求有多少个数对（x，y）满足 a ≤ x ≤ b ,c ≤y ≤ d，且gcd（x,y）=k？
-
输入格式：
第一行为一个整数T，接下来T行每行五个整数，分别表示a，b，c，d，k。
-
输出格式：
共T行，每行一个整数表示满足要求的数对（x，y）的个数
-
样例输入 样例输出
2
2 5 1 5 1 14
1 5 1 5 2 3
-
数据范围：
20%的数据满足：1 ≤T,k ≤100，1 ≤ a ≤ b ≤ 100，1 ≤ c ≤ d ≤ 100。
50%的数据满足：1 ≤T,k ≤1000，1 ≤ a ≤ b ≤ 1000，1 ≤ c ≤ d ≤ 1000。
70%的数据满足：1 ≤T ≤1000，1 ≤ a ≤ b ≤ 10000，1 ≤ c ≤ d ≤ 10000，1 ≤k ≤10000。
100%的数据满足：1 ≤T,k ≤50000，1 ≤ a ≤ b ≤ 50000，1 ≤ c ≤ d ≤ 50000。
-
算法一
-
枚举[a,b]中的每一个数x，枚举[c,d]中的每一个数y，判断gcd(x,y)是否等于k，单次询问的时间复杂度是O(Tnm\(lnn\))。
当然，我们可以对上述暴力的算法进行优化，令x=k×x1，y=k×y1，则\(\lceil\frac{a}{k}\rceil\)≤x1≤\(\lceil\frac{b}{k}\rceil\)，\(\lceil\frac{c}{k}\rceil\)≤x1≤\(\lceil\frac{d}{k}\rceil\)，枚举x1,y1，判断gcd(x1,y1)是否等于1即可，时间复杂度为O(T\(\frac{nmln\frac{n}{k}}{k^2}\))，预计得分20。
-
算法二
设f(n,m)表示有多少个数对（x,y）满足1≤x≤n，1≤y≤m，且gcd(x,y)=k，则原问题的解=f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1)，计算f(n,m)的时候，我们先枚举1到\(\lfloor\frac{n}{k}\rfloor\)的每一个数x，再枚举1到\(\lfloor\frac{m}{k}\rfloor\)的每一个数y，判断gcd(x,y）=1，公式为f(n,m) = \(\sum_{x=1}^{\frac{n}{k}}\sum_{y=1}^{\frac{m}{k}}{gcd(x,y)=1}\)。
设g(x)表示1到\(\lfloor\frac{m}{k}\rfloor\)中与x互质的数的个数，则f（n,m）=\(\sum_{x=1}^{\frac{n}{k}}{g(x)}\)，如果能求出来g(x)，那么就会省去第二层的枚举。如何求g(x)呢？
根据算数唯一分解定理，设x = \(p_1^{a_1}p_2^{a_2}p_3^{a_3}……p_t^{a_t}\)，假设A集合表示1到\(\lfloor\frac{m}{k}\rfloor\)中含\(p_i\)的集合，利用容斥原理：
g(x)=\(\lfloor\frac{m}{k}\rfloor\)—|\(\cup_{i=1}^{t}{A_i}\)|=\(\lfloor\frac{m}{k}\rfloor\)—\(\sum_{i=1}^{t}{|A_i|}\) + \(\sum_{i,j=1}^{t}{|A_i\cap{A_j}|}\)— \(\sum_{i,j,k=1}^{t}{|A_i\cap{A_j}\cap{A_k}|}\)+……
假设x=30，\(\lfloor\frac{m}{k}\rfloor\)=100，\(|A_i|\)={2,3,5}
g(x)=100—\(\lfloor\frac{100}{2}\rfloor\)—\(\lfloor\frac{100}{3}\rfloor\)—\(\lfloor\frac{100}{5}\rfloor\)+\(\lfloor\frac{100}{2*2}\rfloor\)+\(\lfloor\frac{100}{2*5}\rfloor\)+\(\lfloor\frac{100}{3*5}\rfloor\)—\(\lfloor\frac{100}{2*3*5}\rfloor\)=26，这样的话，计算g(x)的时间复杂度为O(\(2^{t(x)}\))（枚举子集的容斥的时间复杂度，如果是本质不同的数集，则时间复杂度为\(2^n\)，n为不同的数字的个数），t(x)为x的质因子的个数，对于50%的数字，n<=1000，按照最小的质因子乘起来，\(2*3*5*7*11\)=2310 > 1000,所里这里的t(x)最多为4，所以整体的时间复杂度为O(T\(\lfloor\frac{n}{k}\rfloor*\sum_{}2^{t(x)}\))，只能过50%的数据
-
算法三
设f(k)表示有多少个数对（x,y）满足1≤x≤n，1≤y≤m，且gcd(x,y)=k，在这里，我们设定n≤m
设g(k)表示1≤x≤n，1≤y≤m，且k|gcd(x,y)的数对的个数，这里gcd(x,y)是k的倍数，例如gcd(x,y)=k,2k,3k,4k……\({\lfloor\frac{n}{k}\rfloor} *k\)，所以g(k)=\(\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}f(k*d)\) ,g(k)的设定非常重要
g(k)的求值简单，找到1—n中k的倍数x，找到1—m中k的倍数y，gcd(x,y)一定是k的倍数，那么1—n中k的倍数有多少个呢？\({\lfloor\frac{n}{k}\rfloor}\)个，同理，1—m中k的倍数有\({\frac{m}{k}}\)个，根据乘法原理，g(k)=\({\lfloor\frac{n}{k}\rfloor}\)*\(\lfloor{\frac{m}{k}}\rfloor\)。
举例说明：k=2，n=4，m=6
n里面2的倍数有2个{2，4}，m里面2的倍数有3个{2, 4, 6}
{2,2}、{2,4}、{4,2}、{4,4}、
-
对于70%和100%的数据，无法解决
引入：莫比乌斯函数
-
定义f(n)和g(n)是定义再正整数集合上的两个函数，满足以下的关系
f(n)=\(\sum_{d|n}{g{(d)}}\)
根据定义，我们可以知道： 根据f(n)推导出g(n)：
f(1)=g(1) g(1)=f(1)
f(2)=g(1)+g(2) g(2)=f(2) - f(1)
f(3)=g(1)+g(3) g(3)=f(3) - f(1)
f(4)=g(1)+g(2)+g(4) g(4)=f(4) - g(2) -g(1)=f(4) -f(2)
f(5)=g(1)+g(5) g(5)=f(5) - f(1)
f(6)=g(1)+g(2)+g(3)+g(6) g(6)=f(6) -f(2) -f(3) +f(1)
f(7)=g(1)+g(7) g(7)=f(7) - f(1)
f(8)=g(1)+g(2)+g(4)+g(8) g(8)=f(8) - f(4)
f(9)=g(1)+g(3)+g(9) g(9)=f(9) - f(3)
f(10)=g(1)+g(2)+g(5)+g(10) g(10)=f(10) - f(5) - f(2) +f(1)
f(11)=g(1)+g(11) g(11)=f(11) - f(1)
f(12)=g(1)+g(2)+g(3)+g(4)+g(6)+g(12) g(12)=f(12) - f(6) -f(4) +f(2)
我们观察发现:
g(n)=\(\sum_{d|n}{f(d)*?}\)
问号( ? )部分与d无关，与\(\frac{n}{d}\)有关，我们设这部分的值为\(\mu\)(\(\frac{n}{d}\))，我们来猜测一下\(\mu\)(\(\frac{n}{d}\))的值：
观察g(2)的值，我们猜测
\(\mu(1)\)= 1 , \(\mu(2)\)= -1, \(\mu(3)\)= -1，\(\mu(4)\)= 0，\(\mu(5)\)= -1，\(\mu(6)\)= 1，
\(\mu(7)\)= -1，\(\mu(8)\)= 0，\(\mu(9)\)= 0，\(\mu(10)\)= 1，\(\mu(11)\)= -1，\(\mu(12)\)= 0
\(\mu(d)\)为莫比乌斯函数，定义如下：
若d=1，\(\mu(1)\)= 1；
若d= \(p_1^{a_1}p_2^{a_2}p_3^{a_3}……p_t^{a_t}\)，只要有一个质数的指数大于1，则\(\mu(d)\)=0；
若d= \(p_1p_2p_3……p_t\)，则\(\mu(d)\)=\((-1)^t\)
莫比乌斯函数的性质
- 性质一：对于任意正整数n，有\(\sum_{d|n}{\mu{(d)}}=\begin{cases}1, \quad if\quad n=1\\0,\quad if\quad n\neq0\\ \end{cases}\)
证明方法一：
- n=1 时，\(\sum_{d|n}{\mu{(d)}}\)=1；
- n>1时，设n=\(p_1^{x_1}p_2^{x_2}p_3^{x_3}……p_k^{x_k}\)，设约数d=\(p_1^{y_1}p_2^{y_2}p_3^{y_3}……p_k^{y_k}\)，如果\(y_i\)的值大于等于2，则\(\mu\)的值为0，所以我们只考虑\(y_i\)为0或者为1的情况，假设约数d中有r个互异的质因子，质因子的指数为1，则对应的\(\mu(d)\)=\((-1)^r\)，这样的约数d有\(C_k^r\)个，所以\(\sum_{d|n}{\mu{(d)}}=\sum_{r=0}^{k}C_k^r(-1)^r*1^{k-r}\)=\(（1-1）^k\)=0。
上面式子可以进行重要替换，将n替换成gcd(x,y)，
- 性质二： 莫比乌斯函数\(\mu{(n)}\)是积性函数
证明：
-
n=1时，\(\mu{(1)}\)=1；
-
如果a和b互质，我们很容易证明\(\mu{(ab)}=\mu{(a)}*\mu{(b)}\)。分以下几种情况讨论：
假设a或者b中有一个的指数大于等于2，则\(\mu{(a)}\)=0或者\(\mu{(b)}\)=0，\(\mu{(ab)}\)=0；
假设a或者b中的质因子都等于1，因为a和b互质，则易证\(\mu{(ab)}=\mu{(a)}*\mu{(b)}\)；
- 如何求莫比乌斯函数？
由于莫比乌斯函数是积性函数，所以我们可以在线性筛中求函数值
void seive(int n){
mu[1]=1;
for(int i=2;i<=n;i++){
if(vis[i]==0){
pri[++cnt]=i;
mu[i]=-1;//说明这是一个质数
}
for(int j=1;j<=cnt;j++){
int t=i*pri[j];//即将被筛掉的数字
if(t>n) break;//线性筛的判断条件
vis[t]=1;
if(i%pri[j]==0) break;//说明pri[j]的指数大于等于2
else mu[t]=-mu[i];//增加了一个质因子pri[j]，所以变成-1
}
}
}
常用的数论函数
我们不必熟悉整个数论函数的体系，我们只需要了解一些常见数论函数即可:
- \(\phi(n)\) 欧拉函数 表示1～n中与n互质的数的个数
- \(μ(n)\) 莫比乌斯函数
- \(d(n)\) 约数个数 表示n的正因子个数
- \(σ(n)\) 约数和函数 表示n的所有正因子之和
- \(ϵ(n)\) 元函数 定义\(ϵ(n)=[n==1]\)
- \(1(n)\) 恒等函数 定义：\(1(n)=1\)
- \(Id(n)\) 单位函数 定义：\(Id(n)=n\)
- \(I_k(n)\) 单位函数 定义：\(I_k(n)=n^k\)
如果一个数论函数 f ，对于任意两个的正互质的正整数a,b，如果满足\(f(a×b)=f(a)×f(b)\)，我们把函数 f 称作积性函数，注意积性函数与积性函数的乘积仍为积性函数。上述数论函数都是积性函数；
若函数f对于任意一个正整数对(a,b)都满足\(f(a×b)=f(a)×f(b)\)，那么我们把函数f称作完全积性函数。上面的5,6,7,8都是完全积性函数。
狄利克雷卷积
狄利克雷卷积：定义一种运算∗，对于函数f，g满足
，我们把这种运算，称作狄利克雷卷积。狄利克雷卷积满足交换律，结合律以及分配律。元函数\(ϵ(n)\)是狄利克雷卷积的乘法单位，元函数可以在狄利克雷卷积中充当单位元，即\(f∗ϵ=f\)(f函数卷上单位函数，等于f函数)。
交换律：\(f*g=g*f\) 。
结合律：\(f*(g*h)=(f*g)*h\)。
分配律：\(f*(g+h)=f*g+f*h\)。
常用的数论函数卷积
-
\(\mu*1=\epsilon\)
\[(\mu*1)(n)=\sum_{d|n}{\mu(d)}*1(\frac{n}{d})=\sum_{d|n}{\mu(d)}=\begin{cases}1, \quad if\quad n=1\\0,\quad if\quad n\neq0\\ \end{cases} \quad=\quad[n==1]=\epsilon \] -
\(\phi*1=I_k\)
\[n=\prod_{i=1}^{k}{p_i^{x_i}}，因为\phi函数为积性函数，只需要证明n^，=p_i^{x_i}时\phi*1=id即可\\ \phi*1=\sum_{d|n^，}{\phi(d)}=\sum_{i=0}^{k}{\phi(p_i^{x_i})}\\ =1+(p_i-1)+p_i*(p_i-1)+p_i^2*(p_i-1)+……+p_i^{x_i-1}*(p_i-1)\\ =p_i^{x_i}=I_k \]
莫比乌斯反演
莫比乌斯反演的公式如下：
证明方法一：在这里先将f函数替换成g函数
接下来考虑以g(i)为主体的等价形式，如何变化呢？
第一步为什么可以变到第二步？
两步之间的转换必须具备两个条件：条件1是把n分成两部分，是一部分是x，一部分是\(\frac{n}{x}\)部分，条件2是这两部分分别属于两个数论函数，且两个函数都需求和，所以无论形式如何，只要满足这两个条件即可。
举个例子，当n=6时
接下来我们考虑\(\sum_{d|\frac{n}{i}}{\mu{(d)}}\)(主要利用\(\mu{}\)函数的性质一)：
- 如果i=n，\(\sum_{d|\frac{n}{i}}{\mu{(d)}}=\mu{(1)}=1\)，所以\(\sum_{}{g(i)}\sum_{d|\frac{n}{i}}{\mu{(d)}}=g(n)\)
- 如果i取小于n的约数，\(\sum_{d|\frac{n}{i}}{\mu{(d)}}=0\)，\(\sum_{}{g(i)}\sum_{d|\frac{n}{i}}{\mu{(d)}}\)=0
综上，\(\sum_{d|n}f(\frac{n}{d})*\mu{(d)}=g(n)\)，得证！
证明方法二：
使用狄利克雷卷积证明：
莫比乌斯反演的变形
- 变形一
根据f(i)的定义，进行下列式子的变化
上式中我们需要枚举k，枚举\(k_1\)。接下来的变化很关键，设T=\(k_1\)*k，那么T的范围是多少呢？考虑两个边界
- 当k=1，i的值为n，\(k_1\)为1，T的值为1
- 当k=n，i的值为1，\(k_1\)为1，T的值为n
所以T的取值范围为[1，n]。
则k是T的约数，上面式子可以进行如下变化：
根据莫比乌斯函数的性质一，可得：
应用
回想引例中的问题，我们令g(i)表示 x\(\in\)[1, n]，y\(\in\)[1，m]，gcd(x,y)=i的数对的个数
令f(i)表示 x\(\in\)[1, n]，y\(\in\)[1，m]，且i|gcd(x,y)的数对的个数，那么f(i)和g(i)正好满足变形一中的式子关系。
所以变形一适用于g(i)很难求，但是\(\sum_{k=1}^{\lfloor\frac{n}{i}\rfloor}g(k*i)\)很好求，可以利用反演来计算g(i)
- 变形二
变形二和变形一的意思其实是一样的。
莫比乌斯反演的性质
性质：
f(n)和g(n)都是定义在正整数集合上的两个函数，并且满足\(f(n)=\sum_{d|n}{g(d)}\)，
则g(n)是积性函数的充分必要条件是f(n)是积性函数
充分性：\(f(n)\)是积性函数\(\longrightarrow\ g(n)\)是积性函数
设gcd(a,b)=1，\(a=\prod{p_i}^{x_i},b=\prod{q_i}^{y_i}\)，根据莫比乌斯反演可得，\(g(n)=\sum_{d|n}{f(\frac{n}{d})}{\mu(d)}\)
\(g(a)=\sum_{d_1|a}{f(\frac{a}{d_1})}{\mu(d_1)},g(b)=\sum_{d_2|b}{f(\frac{b}{d_2})}{\mu(d_2)}，g(a)g(b)=\sum_{d_1|a}{f(\frac{a}{d_1})}{\mu(d_1)}\sum_{d_2|b}{f(\frac{b}{d_2})}{\mu(d_2)}\),
其中gcd(\(d_1,d_2\))=1，gcd(\(\frac{a}{d_1},\frac{b}{d_2}\))=1，同时\(d_1|a, d_2|b, d_1d_2|ab\)，
因为f和\(\mu\)函数都是积性函数，所以
必要性：\(g(n)\)是积性函数\(\longrightarrow\ f(n)\)是积性函数
设gcd(a,b)=1，\(d_1|a,d_2|b\), 那么
引例的算法四
70%的数据满足：1 ≤T ≤1000，1 ≤ a ≤ b ≤ 10000，1 ≤ c ≤ d ≤ 10000，1 ≤k ≤10000。
100%的数据满足：1 ≤T,k ≤50000，1 ≤ a ≤ b ≤ 50000，1 ≤ c ≤ d ≤ 50000。
设f(k)表示有多少个数对（x,y）满足1≤x≤n，1≤y≤m，且gcd(x,y)=k，在这里，我们设定n≤m
设g(k)表示1≤x≤n，1≤y≤m，且k|gcd(x,y)的数对的个数，这里gcd(x,y)是k的倍数，例如gcd(x,y)=k,2k,3k,4k……\({\lfloor\frac{n}{k}\rfloor} *k\)，所以g(k)=\(\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}f(k*d)\) ,g(k)的设定非常重要
根据g(k)和f(k)的关系，完全符合莫比乌斯变形一的式子，所以
\(\mu(d)\)的值，我们可以在线性筛中O(n)中预处理出来，这样单词询问的时间复杂度为\(O(\frac{n}{k})\)，总的时间复杂度为\(O(T*\frac{n}{k})\)，可以得到70%的分数
引例的算法五
观察\(f(k)\)的式子，以n=32，m=40，k=2为例，观察\(\lfloor\frac{n}{kd}\rfloor*\lfloor\frac{m}{kd}\rfloor\)的变化：
如上图所示，蓝色部分为\(\lfloor\frac{n}{kd}\rfloor\)相同的部分，绿色部分为\(\lfloor\frac{m}{kd}\rfloor\)相同的部分，这部分的内容为整除分块的知识点
整除分块
对于\(f(d)=\lfloor\frac{n}{d}\rfloor\)这样的式子，我们需要计算f(d)的值，暴力枚举d的话，时间复杂度为\(O(n)\)
但正如上图中看到的，f(d)的值呈现中一段一段的连续形态，如果我们能找到每一段的左端点和右端点的话，只需要枚举每一段的左右端点即可，不需要全部枚举，那么对于每一段来说，左右端点如何寻找呢？假设左右端点为\(l，r\)
显然，\(l=r+1\)，那么r是多少呢？\(r=\lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor\)，证明如下
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);//找到右边界
……
}
时间复杂度是多少呢？\(O(\sqrt{n})\)，证明如下：
- 当\(1 \leq d \leq \sqrt{n}时，\lfloor\frac{n}{d}\rfloor最多有\sqrt{n}个不同的取值\)
- 当\(\sqrt{n}+1 \leq d \leq n时，由于\lfloor\frac{n}{d}\rfloor \leq \lfloor\sqrt{n}\rfloor,\lfloor\frac{n}{d}\rfloor最多有\sqrt{n}个不同的取值\)
引例中的算法五继续
观察\(f(k)\)的式子，以n=32，m=40，k=2为例，观察\(\lfloor\frac{n}{kd}\rfloor*\lfloor\frac{m}{kd}\rfloor\)的变化：
如上图所示，蓝色部分为\(\lfloor\frac{n}{kd}\rfloor\)相同的部分，绿色部分为\(\lfloor\frac{m}{kd}\rfloor\)相同的部分，红色部分为\(\lfloor\frac{n}{kd}\rfloor和\lfloor\frac{m}{kd}\rfloor\)相同的部分，按照整除分块的知识来处理这个这些式子(其中k是定值)，相同的部分的值是定值，但是\(\mu\)值不一样，我们可以在线性筛中预处理出\(\mu\)函数的前缀和，接下来在\(O(\sqrt{\frac{n}{k}})\)内处理单次询问。100分可以实现。
代码如下：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int t,k,a,b,c,d,cnt;
ll pri[50005],vis[50005],mu[50005],sum[50005];
//solve其实求的是min([1,x],[1,y])之间的数字
void seive(int n){
mu[1]=1;
for(int i=2;i<=n;i++){
if(vis[i]==0){
pri[++cnt]=i;
vis[i]=i;
mu[i]=-1;//说明这是一个质数
}
for(int j=1;j<=cnt;j++){
int t=i*pri[j];//即将被筛掉的数字
if(t>n) break;//线性筛的判断条件
vis[t]=pri[j];//能够筛掉t的最小质因子是pri[j];
if(i%pri[j]==0) break;
else mu[t]=-mu[i];//增加了一个质因子pri[j]，所以变成-1
}
//不能在这里求前缀和，因为缺少mu[1]
}
for(int i=1;i<=n;i++) sum[i]=sum[i-1]+mu[i];//求莫比乌斯函数的前缀和
}
ll solve(int x,int y,int k){
ll ans=0,r1,r2,r;
int rj=min(x,y);
//所以，l到r，其实枚举的是d的不同赋值
for(int l=1;l<=rj;l=r+1){
r1=x/(x/l);//找到第一个右边界
r2=y/(y/l);;//找到第二个右边界
r=min(r1,r2);//是k*d使得他们呈现出分块的性质，k是一个定值，所以其实是d使得它呈现出分块的性质
ans+=(x*1ll/(1ll*l*k))*(y*1ll/(1ll*l*k))*(sum[r]-sum[l-1]);
}
return ans;
}
int main(){
scanf("%d",&t);
seive(50005);
for(int i=1;i<=t;i++){
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
ll ans=solve(b,d,k)-solve(a-1,d,k)-solve(b,c-1,k)+solve(a-1,c-1,k);
printf("%lld\n",ans);
}
return 0;
}
例2：YY的GCD
题目描述：神犇 YY 虐完数论后给傻× kAc 出了一题，给定 N, M求 \(1 \leq x \leq N,1 \leq y\leq M\)且 $\gcd(x,y) $为质数的 (x, y)有多少对？
输入格式：第一行一个整数 T表述数据组数，接下来 T行，每行两个正整数N, M。
输出格式：T行，每行一个整数表示第 i 组数据的结果。
输入： 输出：
2
10 10 30
100 100 2791
数据：\(T=10^4，N, M \leq 10^7\)。
算法一：
算法二：
#include<bits/stdc++.h>
using namespace std;
const int maxn=10000005;
typedef long long ll;
int t,mu[maxn],pri[1000005],sum[maxn],n,m,vis[maxn],cnt,f[maxn];
void seive(int x){
mu[1]=1;
for(int i=2;i<=x;i++){
if(!vis[i]){
pri[++cnt]=i;
vis[i]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>x||vis[i]<pri[j]) break;
vis[tt]=pri[j];
if(i%pri[j]==0) mu[tt]=0;
else mu[tt]=-mu[i];
}
}
for(int j=1;j<=cnt;j++)
//这里是枚举倍数
for(int k=1;k*pri[j]<=x;k++){
f[k*pri[j]]+=mu[k];//
}
for(int i=1;i<=x;i++) sum[i]=sum[i-1]+f[i];
}
ll solve(int x,int y){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=((1ll*x)/l)*((1ll*y)/l)*(1ll*(sum[r]-sum[l-1]));
}
return ans;
}
int main(){
scanf("%d",&t);
seive(10000000);
for(int i=1;i<=t;i++){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=solve(n,m);
printf("%lld\n",ans);
}
}
算法三：
#include<bits/stdc++.h>
using namespace std;
const int maxn=10000005;
typedef long long ll;
int t,pri[1000005],sum[maxn],n,m,vis[maxn],cnt,f[maxn];
int zs[maxn],tc[maxn];
bool flag[maxn];
//其中zs表示i里面不同的质因子的个数,tc表示i里面含有指数为2的质因子的个数
int js(int x){
int ans=1;
for(int i=1;i<=zs[x];i++){
ans=ans*(-1);
}
return ans;
}
int js1(int x){
int ans=0,k=1;
for(int i=1;i<=zs[x];i++){
ans++;
}
for(int i=1;i<=(zs[x]-1);i++){
k=k*(-1);
}
return ans*k;
}
void seive(int x){
for(int i=2;i<=x;i++){
if(!vis[i]){
pri[++cnt]=i;
vis[i]=i;
zs[i]=1;
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>x||vis[i]<pri[j]) break;
vis[tt]=pri[j];
if(flag[i]==true) {
flag[tt]=true;continue;
}
zs[tt]=zs[i];//开始i和tt的质因子个数是一样的
tc[tt]=tc[i];//质因子指数为2的质因子个数也是一样的
if(i%pri[j]) zs[tt]++;
else{
if(i%(pri[j]*pri[j])==0) flag[tt]=true;//说明含有3个pri[j];
else {
tc[tt]++;
if(tc[tt]>=2) flag[tt]=true;
}
}
}
}
for(int i=2;i<=x;i++)
{
if(flag[i]==true) f[i]=0;
else{
if(tc[i]==1) f[i]=js(i);
else f[i]=js1(i);
}
}
for(int i=1;i<=x;i++) sum[i]=sum[i-1]+f[i];
}
ll solve(int x,int y){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=((1ll*x)/l)*((1ll*y)/l)*(1ll*(sum[r]-sum[l-1]));
}
return ans;
}
int main(){
scanf("%d",&t);
seive(10000000);
for(int i=1;i<=t;i++){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=solve(n,m);
printf("%lld\n",ans);
}
}
算法四：
#include<bits/stdc++.h>
using namespace std;
const int maxn=10000005;
typedef long long ll;
int t,pri[1000005],sum[maxn],n,m,vis[maxn],cnt,f[maxn],mu[maxn];
void seive(int x){
mu[1]=1;
for(int i=2;i<=x;i++){
if(!vis[i]){
pri[++cnt]=i;
vis[i]=i;
mu[i]=-1;
f[i]=1;
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>x||vis[i]<pri[j]) break;
vis[tt]=pri[j];
if(i%pri[j]==0) mu[tt]=0,f[tt]=mu[i];
else mu[tt]=-mu[i],f[tt]=mu[i]-f[i];
}
}
for(int i=1;i<=x;i++) sum[i]=sum[i-1]+f[i];
}
ll solve(int x,int y){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=((1ll*x)/l)*((1ll*y)/l)*(1ll*(sum[r]-sum[l-1]));
}
return ans;
}
int main(){
scanf("%d",&t);
seive(10000000);
for(int i=1;i<=t;i++){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=solve(n,m);
printf("%lld\n",ans);
}
}
例3：于神之怒
题目描述：给定$ n,m,k$，计算
对 \(10^9+7\) 取模的结果。
输入格式：本题单测试点内有多组测试数据。
第一行有两个整数，分别表示数据组数T和给定的 k。
接下来 T行，每行两个整数，表示一组数据的n和 m。
输出格式：对于每组数据，输出一行一个整数表示答案。
样例输入 样例输出
1 2
3 3 20
数据规模与约定:
对于全部的测试点，保证$ 1 \leq T \leq 2 \times 10^3$，\(1 \leq n, m, k \leq 5 \times 10^6\)。
算法一
算法二
#include<bits/stdc++.h>
using namespace std;
const int maxn=5000005;
const int mod=1000000007;
typedef long long ll;
int n,m,k,t,g[maxn],pri[maxn],vis[maxn],cnt,f[maxn],sum[maxn];
int mul(int x,int s){
int ans=1;
while(s){
if(s&1) ans=(1ll*ans)*(1ll*x)%mod;
x=(1ll*x)*(1ll*x)%mod;
s=s>>1;
}
return ans;
}
void seive(int x){
g[1]=1;
for(int i=2;i<=x;i++){
if(!vis[i]){
pri[++cnt]=i;
vis[i]=i;
f[cnt]=mul(i,k);
g[i]=(f[cnt]-1+mod)%mod;//f是一个模完了之后的数字，所以很可能等于0或者1
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>x||vis[i]<pri[j]) break;
vis[tt]=pri[j];
//如果出现一次重复的质因子，乘以f函数，一共乘以xi-1次
if(i%pri[j]==0) g[tt]=(1ll*g[i])*(1ll*f[j])%mod;
//新增一个质因子的时候，乘以g函数
else g[tt]=(1ll*g[i])*(1ll*g[pri[j]])%mod;
}
}
for(int i=1;i<=x;i++) sum[i]=(sum[i-1]+g[i])%mod;
}
ll solve(int x,int y){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=(1ll*x/l)*(1ll*y/l)*(1ll*(sum[r]-sum[l-1]+mod))%mod;
}
return ans%mod;
}
int main(){
scanf("%d%d",&t,&k);
seive(5000000);
for(int i=1;i<=t;i++){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=solve(n,m);
printf("%lld\n",ans);
}
}
算法三
例4：Crash的数字表格 / JZPTAB
今天的数学课上，Crash 小朋友学习了最小公倍数（Least Common Multiple）。对于两个正整数 \(a 和 b\)，\(lcm(a,b)\) 表示能同时整除 $a和 b \(的最小正整数。例如，\)lcm(6,8)$=24。
回到家后，Crash 还在想着课上学的东西，为了研究最小公倍数，他画了一张$ n×m\(的表格。每个格子里写了一个数字，其中第 i行第 j列的那个格子里写着数为\)lcm(i,j)$。
看着这个表格，Crash 想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当 \(n 和 m\)很大时，Crash 就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash 只想知道表格里所有数的和 mod 20101009 的值。
输入格式：输入包含一行两个整数，分别表示 n和 m。
输出格式：输出一个正整数，表示表格中所有数的和 mod 20101009的值。
**输入 ** 输出
4 5 122
数据规模与约定
- 对于 30% 的数据，保证 \(n, m \le 10^3\)。
- 对于70% 的数据，保证 \(n, m \le 10^5\)。
- 对于 100% 的数据，保证\(1\le n,m \le 10^7\)。
算法一：
算法二
例5：【SDOI2014】数表
有一张 n×m 的数表，其第 i行第 j列（\(1\le i\le n，1\le j\le m\)）的数值为能同时整除 i和 j的所有自然数之和。给定 a，计算数表中不大于 a 的数之和。
输入格式:输入包含多组数据。输入的第一行一个整数 Q表示测试点内的数据组数
接下来 Q行，每行三个整数 n，m，a（\(|a|\le 10^9\)）描述一组数据。
输出格式:对每组数据，输出一行一个整数，表示答案模 \(2^{31}\)的值。
**输入 ** 输出
2
4 4 3 20
10 10 5 148
数据范围：
\(1\le n,m\le 10^5，1\le Q\le 2\times 10^4\)
算法一
/**
代码小技巧：如何对2^31次方取模
答案要求对2的31次方取模，为什么我们需要将int定义成unsigned int？
我们知道unsigned int的取值范围是0--2^32-1,这样的取值范围正好是对2^32取模
有一个前置知识，如果我们对x取模，可以先对kx取模，最后再对x取模也行
所以，这里的k是2，x是2^31,kx正好是2^32
相当于先落在0—kx-1,0-kx-1会分成k个0-x-1的区间，再对x取模后，就落在0-x-1区间
**/
#include<bits/stdc++.h>
using namespace std;
#define int unsigned int //unsigned int的取值范围为2^32
const int mod=pow(2,31);
const int maxn=100005;
int pri[maxn],mu[maxn],vis[maxn],n,m,q,cnt,f[maxn],mx=100000;
struct Ys{
int id,sum;
}ys[maxn];
struct Qu{
int id,n,m,a;
}qu[20005];
int ans1[20005];
bool cmp(const Qu x,const Qu y){
return x.a<y.a;
}
bool cmp1(const Ys x,const Ys y){
return x.sum<y.sum;
}
void seive(){
mu[1]=1;
for(int i=2;i<=mx;i++){
if(!vis[i]){
pri[++cnt]=i;
mu[i]=-1;
vis[i]=i;
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>mx||vis[i]<pri[j]) break;
vis[tt]=pri[j];
if(i%pri[j]==0) mu[tt]=0;
else mu[tt]=-mu[i];
}
}
for(int i=1;i<=mx;i++){
ys[i].id=i;
for(int j=1;j<=mx/i;j++)
ys[i*j].sum+=i;
}
}
signed lowbit(int x){
return x&(-x);
}
void insert(int x,int val){
for(;x<=mx;x+=lowbit(x))f[x]+=val;
}
int ask(int x){
if(x==0) return 0;
int sum=0;
for(;x>0;x=x-lowbit(x)) sum+=f[x];
return sum;
}
void solve(int x,int y,int id){
int ans=0;
int bj=min(x,y);
for(int l=1,r;l<=bj;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=(x/l)*(y/l)*(ask(r)-ask(l-1));
}
ans1[id]=ans;
}
signed main(){//为什么这里必须这么写？
seive();
scanf("%ud",&q);
for(int i=1;i<=q;i++){
scanf("%u",&qu[i].n);
scanf("%u",&qu[i].m);
scanf("%u",&qu[i].a);
qu[i].id=i;
}
sort(qu+1,qu+1+q,cmp);
sort(ys+1,ys+1+mx,cmp1);
//for(int i=1;i<=20;i++) cout<<mu[i]<<endl;
int st=1;
for(int i=1;i<=q;i++){
/**如果st的约数和小于当前的a，那么在树状数组上更新她的值，但是她的倍数也受影响
所以他的倍数也应该在线段树上进行更新
但是他的倍数不一定被更新完成了**/
while(ys[st].sum<=qu[i].a&&st<=mx){
int k=ys[st].id;
for(int j=k;j<=mx;j+=k){
insert(j,ys[st].sum*mu[j/k]);
}
st++;
}
solve(qu[i].n,qu[i].m,qu[i].id);
}
for(int i=1;i<=q;i++)
printf("%u\n",ans1[i]%mod);
return 0;
}
例6：[SDOI2015]约数个数和
题目描述
设$ d(x)表示x \(的约数个数，给定 n,m，求\)\sum_{i=1}n\sum_{j=1}md(ij)$
输入格式:输入文件包含多组测试数据。第一行，一个整数 T，表示测试数据的组数。接下来的 T行，每行两个整数 n,m。
输出格式:T行，每行一个整数，表示你所求的答案。
输入 输出
2
7 4 110
5 6 121
数据范围
对于 100%的数据，\(1\le T,n,m \le 50000。\)
#include<bits/stdc++.h>
using namespace std;
const int maxn=50005;
typedef long long ll;
int n,m,mu[maxn],sum[maxn],pri[maxn],vis[maxn],cnt,t;
ll g[maxn];
ll f(int x){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=x/(x/l);
ans+=1ll*(r-l+1)*(x/l);
}
return ans;
}
void seive(int x){
mu[1]=1;
for(int i=2;i<=x;i++){
if(!vis[i]){
vis[i]=i;
mu[i]=-1;
pri[++cnt]=i;
}
for(int j=1;j<=cnt;j++){
int tt=i*pri[j];
if(tt>x||vis[i]<pri[j]) break;
vis[tt]=pri[j];
if(i%pri[j]==0) mu[tt]=0;
else mu[tt]=-mu[i];
}
}
for(int i=1;i<=x;i++) sum[i]=sum[i-1]+mu[i];
for(int l=1,r=0;l<=x;l++){
g[l]=f(l);
}
}
ll solve(int x,int y){
ll ans=0;
for(int l=1,r=0;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=1ll*(sum[r]-sum[l-1])*g[x/l]*g[y/l];
}
return ans;
}
int main(){
seive(50000);
scanf("%d",&t);
for(int i=1;i<=t;i++){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=solve(n,m);
printf("%lld\n",ans);
}
}
总结
莫比乌斯反演应用的四大要点：
- 公式推导：利用莫比乌斯反演及其变形
- 等价变化：变量替换，内层变量外移，转换枚举顺序
- 合理利用线性筛预处理各种信息
- 整除分块处理
浙公网安备 33010602011771号