莫比乌斯反演
本文符号说明
([A]=1) 代表表达式 (A) 是真，([A]=0) 代表其为假。
(gcd(x,y)) 代表最大公约数，(lcm(x,y)) 代表最小公倍数。
(Omega(x)) 代表 (x) 的质因子个数。
(p) 如果没有特殊说明均指一个素数。
前置芝士
包括积性函数性质（线性筛），狄利克雷卷积以及数论分块，杜教筛，min_25筛的内容。如果已经会了的可以自行跳过。
除了这些前置芝士你可能还需要会的知识：约数相关、素数相关、组合数学、二项式定理。
积性函数
若一个数论函数（定义域和值域均为整数（这个定义可能不严谨，但不妨碍我们的推到））满足以下性质
则称 (f(x)) 为数论函数。
如果不用 (gcd(x,y)=1) 则称完全积性函数。
一些常见的的积性函数
欧拉函数
(varphi (x)) 代表 ([1,x]) 中与 (x) 互质的数的个数，即 (varphi(x)=sumlimits_{i=1}^x [gcd(i,x)=1])。
通过类容斥原理我们可以的到欧拉函数的快速计算公式：
若 (x=p_1^{c_1} imes p_2^{c_2} imes cdots imes p_m^{c_m})，则 (varphi(x)=x imes frac{p_1-1}{p_1} imes frac{p_2-1}{p_2} imes cdots imes frac{p_n-1}{p_n})。
证明：
设 (p,q) 均为 (x) 的质因子，则每 (p) 个数中只有一个与 (x) 不互质，同理每 (q) 个数中也只有一个与 (x) 不互质。
而 ((p,q)=1)，所以互质的数一定不冲突，就有 (x imes frac{p-1}{p} imes frac{q-1}{q}) 个互质的数。显然这个可以推到 (m) 项。
通过上面的例子我们还可以发现欧拉函数的一个性质：(varphi(p) = p - 1)，其中 (p) 是质数。
例题：求 (sumlimits_{i=1}^n i imes [gcd(n,i)=1])。
即求 (n) 以内所有与 (n) 互质的数的和。根据辗转相减法有 (gcd(a,b) =gcd(a,a-b))，所以互质的数一定是成对出现且和为 (n) 的。所以题目所求即为 (frac{varphi(n) imes n}{2})
根据其求和式，显然欧拉函数为积性函数。根据积性函数的性质我们可以推出以下结论：
对于一个质数 (p)，若 ((x,p)=1)，则 (varphi(xp)=varphi(x)varphi(p)=(p-1)varphi(x))。否则 (p|x)，根据欧拉函数的计算公式 (varphi(xp)=pvarphi(x))。
根据上述两个性质我们又可以得到线性筛欧拉函数的方法，类似欧拉筛质数，还是用最小的质因子更新：
int N;
int prime[MAXN], phi[MAXN], cnt;
bool mark[MAXN];
void sieve() {
phi[1] = 1;
for (int i = 2; i <= N; i++) {
if (!mark[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= cnt && prime[j] * i <= N; j++) {
mark[i * prime[j]] = true;
if (i % prime[j] == 0) {
phi[i * prime[j]] = prime[j] * phi[i];
break;
} else {
phi[i * prime[j]] = (prime[j] - 1) * phi[i];
}
}
}
}
总结一下，如果是积性函数且素数情况易得，倍数情况易推导的函数即可线性筛。
莫比乌斯函数
通过欧拉函数我们已经简单了解了积性函数，现在我们请出我们今天的主角——莫比乌斯函数。
有些教材会把容斥原理与莫比乌斯函数一起讲，原因是莫比乌斯函数其实就是容斥系数，我们先来看莫比乌斯函数的定义。
莫比乌斯函数记作 (mu(x))，若 (x=1) 则 (mu(x)=1)；若 (x) 为若干素数乘积则 (mu(x)=(-1)^{Omega(x)})；其它情况为 (0)。
先来证明其为积性函数：
对于两个数 ((x,y)=1)，若 (mu(x)=0) 或 (mu(y)=0) 则 (mu(xy)) 一定也为 (0)，得证。若 (x=1) 或 (y=1) 也显然成立。
我们主要考虑 (mu(x)=(-1)^{Omega(x)},mu(y)=(-1)^{Omega(y)}) 这种情况。((x,y)=1) 有 (Omega(xy)=Omega(x)Omega(y))，所以显然 (mu(xy)=mu(x)mu(y))。
为了线性筛出莫比乌斯函数所以我们还需找到 (mu(p))，显然 (mu(p)=-1)。
根据积性函数的性质若 ((p,x)=1)，则 (mu(px)=-mu(x))，而当其为 (p) 的倍数时直接就是 (0)。
int N;
int prime[MAXN], mu[MAXN], cnt;
bool mark[MAXN];
void sieve() {
mu[1] = 1;
for (int i = 2; i <= N; i++) {
if (!mark[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i <= N; j++) {
mark[i * prime[j]] = true;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
} else {
mu[i * prime[j]] = -mu[i];
}
}
}
}
约数个数与约数和函数
题目中有可能会出现这里就简单介绍一下求和式，不过相信大家小学奥数应该也都学过。
若 (n=p_1^{c_1} imes p_2^{c_2} imes cdots imes p_m^{c_m})。
约数个数函数：(d(n)=(c_1+1) imes (c_2+1) imes cdots imes (c_m+1))。
约数和函数：(sigma(n)=(sum_{i=0}^{c_1} p_1^i) imes cdots imes(sum_{i=0}^{c_m}p_m^i))。
性质的证明与线性筛的实现读者可以自己思考。
一道练习题（数据还没造完）
这里再提一嘴埃氏筛，这玩意虽然 (O(nlog{n}))，但有的时候是即好想又好用，如果时间允许写这个比想破脑袋写线性筛好多了。其实就是用每个数去筛其倍数。
上面那道题没卡埃氏筛。
一些常见的完全积性函数
常函数 (I(n)=1)。
恒等函数 (id(n)=n)。
单位函数 (varepsilon(n)=[n=1])。
狄利克雷卷积
对于三个数论函数 (F(x),f(x),g(x))，若 (F(x)=sum_{d|x}f(x)g(frac{x}{d}))，则称 (F) 是 (f) 和 (g) 的狄利克雷卷简记 (F=f*g)。显然 (f*g=g*f)。
性质1：(f * varepsilon=f)。
证明：由狄利克雷卷积的定义就可以证明。
性质2：(mu * I=varepsilon)。
证明：稍微复杂，展开其实是 (sum_{d|n}mu(d))，当 (n=1) 时 (mu(n)=1) 得证；而 (n e 1) 时，因为对于不是若干质数乘积的数 (mu) 为 (0)，所以我们只考虑是若干质数乘积的 (d)，所以我们只用考虑选或不选某个质因子，从而转变为在 (Omega(n)) 个数中选 (0,1,2dotsOmega(n)) 个数的方案。
所以有：
得证。这条性质很重要，可能在做题时比一般的反演用得还多。反演本质。
性质3：(varphi * I = id)。
证明：两种证明方法，这里用逆推法证，另外一种用积性函数证的可以参考李煜东的蓝书。
即证：(sum_{d|n}varphi(d)=n)。我们把 ([1,n]) 的数与和 (n) 的 (gcd) 分类，即 (n=sum_{i|n}sum_{j=1}^n [gcd(j,n)=i])，再根据 (gcd) 的性质，有 (sum_{i|n}sum_{j=1}^n [gcd(j,n)=i]=sum_{i|n}sum_{j=1}^{frac{n}{i}} [gcd(j,frac{n}{i})=1]=sum_{i|n} varphi(frac{n}{i})=sum_{i|n}varphi(i))。
得证。这条结论还被称为欧拉反演。
另一种证明方法：
设 (f = varphi * I)，易得 (f) 是积性函数。(f(p^m)=sum_{i=1}^{m}varphi(p^i)=p^m)，所以 (f(n)=n)，(f=id)。
性质4：(mu*id=varphi)。
这其实已经是一个莫比乌斯反演的形式了。通过性质2和性质3可得：(mu*I*id=varepsilon*id=id=varphi*I)，然后两边同时搞掉 (I) 就是性质4的式子了。
性质5、6：
做题可能会用到
(I*I=d)。
(I*id=sigma)。
数论分块
求 (sumlimits_{i=1}^nlfloorfrac{n}{i} floor)。([i,n/(n/i)]) 的数都是相等的。最多只会有 (2sqrt{n}) 块，所以可以 (O(sqrt{n})) 解决。
for (int i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ans += (j - i + 1) * (n / i);
}
莫比乌斯反演
终于到正题了qwq。
莫比乌斯反演的约数形式：(F=f*I)，则(f=F*mu)。
证明和性质4很像：(F*mu=f*mu*I=f*epsilon=f)。
倍数形式：(F(n)=sum_{n|d}f(d))， (f(n)=sum_{n|d}F(d)mu(d/n))。
好了，没了，就这，就这！！！
但是毒瘤出题人可以把莫比乌斯反演出到你怀疑人生。
杜教筛
我们在数论分块的时候经常要用到前缀和，而有的时候毒瘤出题人喜欢把 (n) 出到 (10^{10}) 或更高，现在线性筛也没用了，我们需要一个更快的求出积性函数前缀和的方法。
假设我们现在想求 (f) 的前缀和，设为 (S)，然后我们找到一个积性函数 (g) 使得 (f*g) 的前缀和很好求，(g) 的前缀和也很好求。假设我们现在已经找到了这么一个 (g)。
可以把 (d) 提到前面来转而枚举 (d)：
我们现在想求 (g(1)S(n))，我们发现直接一减就是：
前面的很好求（这是要求），后面的可以数论分块然后递归去做。
时间复杂度的分析：
如果通过线性筛预处理前面的一部分，假设是 (k)。
则为 (Theta(frac{n}{sqrt{k}}))，显然 (k=n^{2/3}) 最优，复杂度为 (n^{2/3})。
看到 (mu) 想到 (mu*I=varepsilon)，然后 (varepsilon) 和 (I) 的前缀和都很好求，就做完了，注意记忆化。
(varphi) 的自己思考。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 2000000;
int t;
ll prime[MAXN], cnt;
ll sum_mu[MAXN], sum_phi[MAXN];
bool mark[MAXN];
map<ll, ll> mp_mu, mp_phi;
ll getsum_mu(ll n) {
if (n < MAXN) return sum_mu[n];
if (mp_mu.find(n) != mp_mu.end()) return mp_mu[n];
ll ret = 1;
for (ll i = 2, j; i <= n; i = j + 1) {
j = n / (n / i);
ret -= 1ll * (j - i + 1) * getsum_mu(n / i);
}
return mp_mu[n] = ret;
}
ll getsum_phi(ll n) {
if (n < MAXN) return sum_phi[n];
if (mp_phi.find(n) != mp_phi.end()) return mp_phi[n];
ll ret = n * (n + 1) >> 1;
for (ll i = 2, j; i <= n; i = j + 1) {
j = n / (n / i);
ret -= 1ll * (j - i + 1) * getsum_phi(n / i);
}
return mp_phi[n] = ret;
}
void sieve() {
sum_mu[1] = sum_phi[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!mark[i]) {
prime[++cnt] = i;
sum_mu[i] = -1;
sum_phi[i] = i -1 ;
}
for (int j = 1; j <= cnt && prime[j] * i < MAXN; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
sum_mu[i * prime[j]] = 0;
sum_phi[i * prime[j]] = sum_phi[i] * prime[j];
break;
} else {
sum_mu[i * prime[j]] = -sum_mu[i];
sum_phi[i * prime[j]] = sum_phi[i] * (prime[j] - 1);
}
}
}
for (int i = 2; i < MAXN; i++) sum_mu[i] += sum_mu[i - 1], sum_phi[i] += sum_phi[i - 1];
}
int main() {
sieve();
read(t);
while (t--) {
ll n;
read(n);
printf("%lld %lld
", getsum_phi(n), getsum_mu(n));
}
return 0;
}
例题
[POI2007]ZAP-Queries
题目即求 (sum_{i=1}^asum_{j=1}^b[gcd(i,j)=d])。
这里介绍几种不同的算法，其实不同只在推到过程，最后的式子都是相同的。
这里不妨设 (a ge b)。
算法1
只用到了上述性质2。根据 (gcd) 的性质，有(gcd(i,j)=d Rightarrow gcd(i/d,j/d)=1)。
所以这题所求可以转换为
继续推为
可以把 (mu(k)) 提到前面来：
显然这个东西可以数论分块来搞，时间复杂度 (O(sqrt{min(a,b)}))。
算法2
看完算法1我们发现这和莫比乌斯反演莫得关系啊。那来个有关系的？
设 (f(d)) 为所求。我们需要找到一个合适的 (F(d))，这里我们发现若选取 (F(n)=sum_{n|d}f(d)=[frac{a}{n}] cdot[frac{b}{n}])，比较好搞。
根据莫比乌斯反演的倍数公式 ：
于是我们得到了....和上面一样的东西！
惊不惊喜，意不意外！
但我觉得还是算法1更无脑一点，毕竟算法2需要选择合适的 (F)。不知道有没有神仙教教我怎么快速找 (F) 啊。
不过你说算法1叫莫比乌斯反演也可，因为其用到了莫比乌斯反演最基础的公式（反演本质），而且最后一步也推到了类似莫比乌斯反演公式的东西。算法1最重要的核心是找到 (varepsilon)，也就是 ([...=1]) 这种东西。
对于看算法2爽的同学，对于这种 (gcd) 的题 (F) 一般都是这么选（倍数形式）的。
这道题的代码实现：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 50010;
ll mu[MAXN], prime[MAXN], cnt, sum[MAXN];
bool mark[MAXN];
void sieve() {
mu[1] = sum[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!mark[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i < MAXN; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
} else {
mu[i * prime[j]] = -mu[i];
}
}
}
for (int i = 2; i < MAXN ; i++) sum[i] = mu[i] + sum[i - 1];
}
ll solve(int a, int b, int d) {
a /= d, b /= d;
if (a < b) swap(a, b);
ll ret = 0;
for (int i = 1, j; i <= b; i = j + 1) {
j = min(a / (a / i), b / (b / i));
ret += (sum[j] - sum[i - 1]) * (long long)(a / i) * (long long)(b / i);
}
return ret;
}
int main() {
sieve();
int n;
scanf("%d", &n);
while (n--) {
int a, b, d;
scanf("%d%d%d", &a, &b, &d);
printf("%lld
", solve(a, b, d));
}
return 0;
}
YY的GCD
这题和上题类似，不过怕还有同学（比如说我）不会推式子就再推一遍。
求：(sum_{pin Prime}sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=p])。
老套路，把 (k) 搞到前面来：
这里看 (pk) 不爽，再把它搞到前面去，其实也就是把 (p) 搞到里面去：
后面那玩意可以用埃氏筛求出来，时间复杂度约为：(O(frac{n}{ln{n}} imes ln{{n}})=O(n))。当然还可用标准的线性筛（那也是个积性函数），不过没有必要（上面有提到）。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 10000010;
int T;
int N, M;
int mu[MAXN], prime[MAXN], cnt;
ll sum[MAXN];
bool mark[MAXN];
void sieve() {
mu[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!mark[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i < MAXN; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
} else {
mu[i * prime[j]] = -mu[i];
}
}
}
for (int i = 1; i <= cnt; i++)
for (int j = 1; j * prime[i] < MAXN; j++)
sum[j * prime[i]] += mu[j];
for (int i = 1; i < MAXN; i++) sum[i] += sum[i - 1];
}
ll solve(int n, int m) {
ll ret = 0;
if (n > m) swap(n, m);
for (int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
ret += (sum[j] - sum[i - 1]) * (n / i) * (m / i);
}
return ret;
}
int main() {
sieve();
scanf("%d", &T);
while (T--) {
scanf("%d%d", &N, &M);
printf("%lld
", solve(N, M));
}
return 0;
}
类似的题：[HAOI2011]Problem b（套上简单容斥）、[SDOI2014]数表（树状数组维护前缀和）
[SDOI2015]约数个数和
这题虽然也和上面类似，但是有一个引理还是很有意思的。
引理：(d(ij)=sum_{a|i}sum_{b|j}[gcd(i,j)=1])。
证明：考虑一个质数 (p)，设 (i=i' imes p^{k_i},j=j' imes p^{k_j})。由上面的约数个数公式可得 (p) 的贡献为 (k_i+k_j+1)。考虑等式右边，设 (a=a' imes p^{k_a},b=b' imes p^{k_b})，若有 (gcd(a,b)=gcd(a' imes p^{k_a},b' imes p^{k_b})=1)，必有 (k_a=0) 或 (k_b=0)。而前者和后者分别还有 (k_j+1) 和 (k_i+1) 种，减去都为 (0) 的一种，共 (k_i+k_j+1) 种，刚好为左边。
然后我们处理的式子即为：
然后就是套路：
然后 (a,b) 作为约数没用，提到前面去：
后面那俩玩意可以预处理，然后就可以数论分块 (O(sqrt{n})) 做了。
[国家集训队]Crash的数字表格
根据小学数学 (lcm(a,b)=frac{a imes b}{gcd(a,b)})，然后套路地枚举 (gcd)，得到小面的式子：
发现是两个数论分块套一起，时间复杂度 (O(n))。注意一下二维数论分块时一定要判 (n,m) 大小关系。
当然这题还有 (O(sqrt{n})) 的解法，具体来说枚举第二个分块的除数为 (T)，然后后面搞出来的是个积性函数，可以线性筛。当你想降一下时间复杂度时，枚举第二个分块中的某一项再进行处理可能是一个好选择。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 10000010;
const ll mod = 20101009;
int T;
int prime[MAXN], mu[MAXN], cnt;
ll sum[MAXN];
bool mark[MAXN];
void sieve() {
mu[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!mark[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i < MAXN; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
} else {
mu[i * prime[j]] = -mu[i];
}
}
}
for (ll i = 1; i < MAXN; i++) sum[i] = sum[i - 1] + mu[i] * i * i, sum[i] %= mod;
}
ll calc(ll a, ll b) {
return ((a + b) * (b - a + 1) >> 1) % mod;
}
ll count(int n, int m) {
ll ret = 0;
if (n > m) swap(n, m);
for (int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
ret += (sum[j] - sum[i - 1] + mod) * calc(1, n / i) % mod * calc(1, m / i) % mod;
ret %= mod;
}
return ret % mod;
}
ll solve(int n, int m) {
ll ret = 0;
if (n > m) swap(n, m);
for (int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
ret += calc(i, j) * count(n / i, m / i) % mod;
ret %= mod;
}
return ret % mod;
}
int main() {
sieve();
int n, m;
scanf("%d%d", &n, &m);
printf("%lld
", solve(n, m));
return 0;
}
简单的数学题
欧拉反演×min_25筛或杜教筛
发现这玩意一顿乱推需要我们求一个这样的东西的前缀和：(id^2varphi)。于是我们要找到一个合适的 (g)，看到 (varphi) 我们想到 (I)，但是还有个 (id^2) 所以我们也要再搞一个 (id^2)。于是我们的 (g) 就是 (id^2I)，(((id^2varphi) * (id^2I))(n)=n^2(varphi*I)(n)=n^3)，前缀和就是个立方和公式，相信大家小学奥数都学过。然后 (g) 的前缀和就是个平方和公式，也可以求，于是就可以杜教筛了。杜教筛的时间复杂度不变。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 5000000;
ll mod, inv2, inv6;
int T;
ll prime[MAXN + 1], phi[MAXN + 1], cnt;
ll sum[MAXN + 1];
bool mark[MAXN + 1];
map<ll, ll> mp;
ll ksm(ll a, ll b) {
ll ret = 1;
while (b) {
if (b & 1) ret = (ret * a) % mod;
a = (a * a) % mod;
b >>= 1;
}
return ret;
}
void sieve() {
phi[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!mark[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= cnt && prime[j] * i < MAXN; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
} else {
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
}
for (ll i = 1; i < MAXN; i++) sum[i] = (sum[i - 1] + phi[i] * i % mod * i) % mod;
}
ll calc(ll n) {
n %= mod;
return n * (n + 1) % mod * inv2 % mod;
}
ll cal(ll n) {
n %= mod;
return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
ll get_sum(ll n) {
if (n < MAXN) return sum[n];
if (mp.find(n) != mp.end()) return mp[n];
ll ret = calc(n) * calc(n) % mod;
for (ll i = 2, j; i <= n; i = j + 1) {
j = n / (n / i);
ret -= (cal(j) - cal(i - 1) + mod) * get_sum(n / i) % mod;
ret = (ret + mod) % mod;
}
ret = (ret + mod) % mod;
return mp[n] = ret;
}
ll solve(int n) {
ll ret = 0;
for (int i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ret += (get_sum(j) - get_sum(i - 1) + mod) * calc(n / i) % mod * calc(n / i) % mod;
ret %= mod;
}
return ret;
}
int main() {
ll p, n;
scanf("%lld%lld", &p, &n);
mod = p;
sieve();
inv2 = ksm(2, mod - 2); inv6 = ksm(6, mod - 2);
printf("%lld
", solve(n));
return 0;
}
最小公倍数之和
这时不一定是 (1) 到 (n) 了，可是值域很小，所以我们不妨看成这样：
设 (c_i) 代表 (i) 出现的次数，然后求的即为 (sum_{i=1}^{50000}sum_{j=1}^{50000}c_ic_jlcm(i,j))。然后莫比乌斯反演套路就可以处理了。
类似题目：[AGC038C] LCMs
[SDOI2017]数字表格
最后一题了。
我们发现它求的是积，不过没有关系，还是一样的套路。
预处理 (Pi_{d|T}f_d^{mu(frac{T}{d})}) 的时间是 (O(nlog{n})) 的。然后可以数论分块，单次 (O(sqrt{n}))。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9 + 7;
const int MAXN = 1000010;
template <typename T> void read(T &x) {
T f = 1;
char ch = getchar();
for (; !isdigit(ch); ch = getchar()) if (ch == '-') f = -1;
for (x = 0; isdigit(ch); ch = getchar()) x = x * 10 + ch - '0';
x *= f;
}
int T;
ll prime[MAXN], mu[MAXN], cnt;
ll sum[MAXN], F[MAXN], invF[MAXN];
bool mark[MAXN];
ll ksm(ll a, ll b) {
ll ret = 1;
while (b) {
if (b & 1) ret = (ret * a) % mod;
a = (a * a) % mod;
b >>= 1;
}
return ret;
}
void sieve() {
for (int i = 1; i <= 1000000; i++) sum[i] = 1;
F[0] = 0;
F[1] = 1;
invF[1] = 1;
for (int i = 2; i <= 1000000; i++) F[i] = (F[i - 1] + F[i - 2]) % mod, invF[i] = ksm(F[i], mod - 2);
mu[1] = 1;
for (int i = 2; i <= 1000000; i++) {
if (!mark[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i <= 1000000; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
mu[i * prime[j]] = -mu[i];
}
}
sum[0] = 1;
for (int i = 1; i <= 1000000; i++) {
for (int j = i; j <= 1000000; j += i) {
if (mu[j / i] > -1) sum[j] = (sum[j] * ksm(F[i], mu[j / i])) % mod;
else sum[j] = (sum[j] * invF[i]) % mod;
}
}
for (int i = 2; i <= 1000000; i++) {
sum[i] = (sum[i - 1] * sum[i]) % mod;
}
}
ll solve(ll n, ll m) {
ll ret = 1;
if (n > m) swap(n, m);
for (int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
ret = (ret * ksm(sum[j] * ksm(sum[i - 1], mod - 2) % mod, (n / i) * (m / i))) % mod;
}
return ret;
}
int main() {
sieve();
read(T);
while (T--) {
ll n, m;
read(n); read(m);
printf("%lld
", solve(n, m));
}
return 0;
}
处理技巧
- (sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=k]=sum_{i=1}^{n/k}sum_{j=1}^{m/k}[gcd(i,j)=1])。
- 枚举约数不好搞，不如把约数提前枚举然后枚举倍数，因为倍数很容易算。(n) 以内 (x) 的倍数有 (lfloorfrac{n}{x} floor) 个。
- 对于一个整体 (xy)，我们可以把它提到前面去，转而枚举其约数然后形成类卷积形式（一般都是积性函数）埃氏筛或线性筛都可搞。
- 求的是一个数列可以转换为出现次数
- 乘积转换为幂
参考资料
莫比乌斯反演+常见数论函数的性质+狄利克雷卷积+数论分块+杜教筛学习笔记 from birchtree
算法竞赛进阶指南 from 李煜东
题解 P3455 【[POI2007]ZAP-Queries】 from pengym
莫比乌斯反演-让我们从基础开始 from An_Account
完结撒花❀