莫比乌斯反演学习笔记
前置知识
目录
一、莫比乌斯函数
- 定义
- 性质
- 证明
二、莫比乌斯反演
- 反演公式
- 形式 \(1\)
- 证明
- 形式 \(2\)
三、例题
- POI2007 zap-Queries
- YY 的 GCD
- SDOI2014 数表
- BZOJ3309 DZY Loves Math
- SDOI2015 约数个数和
- SDOI2017 数字表格
- BZOJ4407 于神之怒加强版
- 国家集训队 Crash 的数字表格、JZPTAB
莫比乌斯函数
定义
\[\mu(n)=
\begin{cases}
1&n=1\\
0&n\text{ 含有平方因子}\\
(-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}\\
\end{cases}
\]
设 \(n=\prod^{k}_{i=1}p_i^{c_i}\)，其中 \(p_i\) 为 \(n\) 的质因子，\(c_i \geq 1\)。
性质
莫比乌斯函数是积性函数，有以下性质
\[\sum_{d \mid n}\mu(d)=\left[n=1 \right]
\]
常用于莫比乌斯反演中。
证明
- 当 \(n=1\) 时显然成立；
- 当 \(n \neq 1\) 时，设
\[n=\prod^{k}_{i=1}p_i^{c_i}
\]
那么有
\[\sum_{d \mid n}\mu(d)=\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i
\]
考虑我们选 \(i\) 个本质不同的质因子，有 \(\binom{k}{i}\) 种方法，其中每个数的莫比乌斯函数值为 \((-1)^i\)，所以这样计算出原式的贡献。
那么我们需要证明
\[\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i=0
\]
考虑二项式定理
\[\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i=\sum^{n}_{i=0}\dbinom{n}{i} \cdot (-1)^i \cdot 1^{n-i}=(-1+1)^n=0
\]
证毕。
莫比乌斯反演
反演公式
设 \(f(n),g(n)\) 是定义在正整数集上的两个函数。
形式 \(1\)
若有
\[f(n)=\sum_{d \mid n}g(d)=\sum
_{d \mid n}g(\dfrac{n}{d})\]
则有
\[g(n)=\sum_{d \mid n}f(\dfrac{n}{d})\mu(d)=\sum_{d \mid n}f(d)\mu(\dfrac{n}{d})
\]
证明
考虑对下式进行变换
\[\sum_{d \mid n}\mu(d)f(\dfrac{n}{d})
\]
代入定义，有
\[\sum_{d \mid n}\mu(d)f(\dfrac{n}{d})=\sum_{d\mid n}\mu(d)\sum_{t\mid \frac{n}{d}}g(t)
\]
可以交换求和号
\[\sum_{t\mid n}g(t)\sum_{d\mid \frac{n}{t}}\mu(d)
\]
详细解释
对于上面这一步变换，有以下考虑。由于 \(d \mid n\)，存在 \(t\) 满足 \(t \mid \! \frac{n}{d}\)。那么对于每一个符合条件的 \(t\) 都有与之对应的 \(d\)，那么就可以转化为
\[\sum_{dt \mid n}\mu(d)g(t)
\]
考虑原式，相当于先枚举符合条件的 \(d\)，再枚举在此时 \(d\) 基础上符合条件的 \(t\)，我们发现，换成现在的式子后，对答案的统计仍然不重不漏。
还可以换一种方式考虑，对于原式，我们可以转化为
\[\sum_d \mu(d) \left[ d \mid n\right]\sum_t g(t)\left[dt \mid n\right]
\]
显然，对于两个限制条件，我们可以直接消去第一个条件，再将艾弗森括号转换进和式，得到
\[\sum_{dt \mid n}\mu(d)g(t)
\]
得到这个式子，我们就可以直接交换 \(\mu(d)\) 和 \(g(t)\) 的位置了，并可以将和式中的一个条件拆开成为两个和式
\[\sum_{t \mid n}g(t)\sum_{d \mid \frac{n}{t}}\mu(d)
\]
又由 \(\sum_{d\mid n}\mu(d)=\left[ n=1 \right]\) 可以得到
\[=\sum_{t \mid n}g(t)\left[ \frac{n}{t}=1 \right]
\]
把艾弗森括号里的部分转化一下
\[=\sum_{t \mid n}g(t)\left[ n=t \right]
\]
只有当 \(n=t\) 时，该式子才有值，即
\[=g(n)
\]
证毕。
形式 \(2\)
若有
\[f(n)=\sum_{n \mid d}g(d)=\sum_{n \mid d}g( \frac{n}{d})
\]
则有
\[g(n)=\sum_{n \mid d}\mu( \frac{d}{n})f(d)=\sum_{n \mid d}\mu(d)f( \frac{d}{n})
\]
证明与形式 \(1\) 的证明类似，不再展开说了。
例题
POI2007 ZAP-Queries
题目大意
给定 \(n\) 和 \(m\)，求
\[\sum_{i = 1}^n\sum_{j = 1}^m \left[ \gcd(i,j) = d\right]
\]
证明
为了方便处理，我们钦定 \(n \leq m\)。
首先考虑 \(d\) 同时是 \(i\) 和 \(j\) 的最大公约数，所以满足 \(d \ | \ i\) 且 \(d \ | \ j\)。
于是我们可以把 \(d\) 提出来。
\[\sum_{i = 1}^{\left\lfloor\ \frac{n}{d} \ \right\rfloor}\sum_{j = 1}^{\left\lfloor\ \frac{m}{d} \ \right\rfloor} \left[ \gcd(i,j) = 1\right]
\]
然后我们可以利用莫比乌斯函数的结论，变为
\[\sum_{i = 1}^{\left\lfloor\ \frac{n}{d} \ \right\rfloor}\sum_{j = 1}^{\left\lfloor\ \frac{m}{d} \ \right\rfloor}\sum_{t |\gcd(i,j)} \mu(t)
\]
此时，实际上 \(t\) 是随意枚举的，只要保证 \(t\ | \gcd(i,j)\) 就可以。既然 \(t\) 不受到前面两个和式的限制，就可以向前提，先枚举 \(t\)，又考虑到满足 \(t\ | \gcd(i,j)\) 的最大的值也就是 \(\min(\left\lfloor\ \frac{n}{d} \ \right\rfloor,\left\lfloor\ \frac{m}{d} \ \right\rfloor)\)，所以 \(t\) 就枚举到这就行。
又考虑对 \(\mu(t)\) 产生限制的只有 \(t\)，那么就可以把她也向前提前面去，放在 \(t\) 的和式后面。我们也可以顺便把其他和式的限制条件放在她后面。
\[\sum_{t = 1}^{\left\lfloor\ \frac{n}{d} \ \right\rfloor}\mu(t) \sum_{i = 1}^{\left\lfloor \frac{n}{d} \ \right\rfloor}\left[t \ | \ i \right]\sum_{j=1}^{\left\lfloor\ \frac{m}{d} \ \right\rfloor} \left[t \ | \ j \right]
\]
然后我们可以发现，观察后面两个和式。例如第二个和式，满足 \(\left[t \ | \ i \right]\) 的 \(i\) 最多有 \(\left\lfloor\ \frac{n}{dt} \ \right\rfloor\) 个。
那么我们这个式子又可以改变，变成下面这个样子。
\[\sum_{t = 1}^{\left\lfloor\ \frac{n}{d} \ \right\rfloor}\mu(t) \left\lfloor\ \frac{n}{dt} \ \right\rfloor \left\lfloor\ \frac{m}{dt} \ \right\rfloor
\]
对于后面两个向下取整的部分，我们可以整除分块求，前面的部分我们可以预处理求出莫比乌斯函数。
Code
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
#define int long long
const int N = 50050;
int T,n,m,d;
int pri[N],p[N],cnt;
int miu[N],mu[N];
void Ready() {
p[1] = 1;
mu[1] = 1;
for(int i = 2;i < N; i++) {
if(p[i] == 0) {
pri[++cnt] = i;
mu[i] = -1;
}
for(int j = 1;j <= cnt && i * pri[j] < N; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0)
break;
mu[i * pri[j]] = -mu[i];
}
}
for(int i = 1;i < N; i++)
miu[i] = mu[i] + miu[i - 1];
return ;
}
int ans = 0;
signed main() {
cin >> T;
Ready();
while(T--) {
cin >> n >> m >> d;
if(n > m)
swap(n,m);
n /= d;
m /= d;
int l = 1,r = 0;
while(l <= n) {
r = min(n / (n / l),m / (m / l));
ans += (n / l) * (m / l) * (miu[r] - miu[l - 1]);
l = r + 1;
}
cout << ans << "\n";
ans = 0;
}
return 0;
}
YY 的 GCD
题目大意
给定 \(n\) 和 \(m\)，求
\[\sum_{i=1}^n\sum^{m}_{j=1}\left[\gcd(i,j) \in \mathbb{P} \right]
\]
即求最大公约数为质数的数对个数。
思路
为了方便处理，我们还是钦定 \(n \leq m\)。
我们可以枚举一个 \(d\)，把原式转化为
\[\sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]\sum_{i=1}^n\sum^{m}_{j=1}\left[\gcd(i,j) = d \right]
\]
利用莫比乌斯函数性质，有
\[\sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]\sum_{i=1}^{\left\lfloor \frac{n}{d}\right\rfloor} \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\sum^{}_{t \mid \gcd(i,j)}\mu(t)
\]
交换求和号，令 \(i= \frac{i}{t},j= \frac{j}{t}\)，有
\[\begin{aligned}
&= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1}\mu(t)
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}\left[i \mid t\right]
\sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\left[j \mid t\right] \\
&= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1}\mu(t)
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}\left[i \mid t\right]
\sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\left[j \mid t\right] \\
&= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1}
\mu(t)
{\left\lfloor \frac{n}{dt}\right\rfloor}
{\left\lfloor \frac{m}{dt}\right\rfloor} \\
\end{aligned}
\]
式子推到这里，如果我们去枚举 \(d\)，肯定会 TLE，我们令 \(T = dt\)，有
\[\sum^{n}_{T=1}{\left\lfloor \frac{n}{T} \right\rfloor}{\left\lfloor \frac{m}{T} \right\rfloor}\sum^{}_{d \mid T}\left[ d \in \mathbb{P} \right]\mu(\frac{T}{d})
\]
我们发现，后面那部分可以预处理。
考虑每个质数 \(d\)，对于它的倍数 \(T\)，其值加上 \(\mu( \frac{T}{d})\)。
Code
#include <bits/stdc++.h>
using namespace std;
const int N = 10050000;
int T,n,m;
bool p[N];
int mu[N],sum[N],f[N];
int prime[N],cnt;
void Seive() {
mu[1] = 1;
for(int i = 2;i <= 10000000; i++) {
if(!p[i]) {
cnt ++;
prime[cnt] = i;
mu[i] = -1;
}
for(int j = 1;i * prime[j] <= 10000000 && j <= cnt; j++) {
p[i * prime[j]] = 1;
if(i % prime[j] == 0)
break;
mu[i * prime[j]] = - mu[i];
}
}
for(int i = 1;i <= cnt; i++)
for(int j = 1;prime[i] * j <= 10000000; j++)
f[prime[i] * j] += mu[j];
for(int i = 1;i <= 10000000; i++)
sum[i] = sum[i - 1] + f[i];
}
long long H(int a,int b) {
long long res = 0;
if(a == 0)
return 0;
int l = 1,r = 0;
while(l <= n) {
r = min(a / (a / l),b / (b / l));
res += (sum[r] - sum[l - 1]) * (long long)(a / l) * (long long)(b / l);
l = r + 1;
}
return res;
}
int main() {
ios::sync_with_stdio(false);
cin >> T;
Seive();
while(T --> 0) {
cin >> n >> m;
if(n > m)
swap(n,m);
cout << H(n,m) << "\n";
}
return 0;
}
SDOI2014 数表
题目大意
有一张 \(n\times m\) 的数表，其第 \(i\) 行第 \(j\) 列（\(1\le i\le n\)，\(1\le j\le m\)）的数值为能同时整除 \(i\) 和 \(j\) 的所有自然数之和。给定 \(a\)，计算数表中不大于 \(a\) 的数之和。
（\(T \leq 2 \times 10^4,1 \leq n,m \leq 10^5,a\leq 10^9\)）。
思路
形式化题意：
给定 \(n,m,a\)，求
\[\sum^{n}_{i=1}
\sum^{m}_{j=1}
\sigma_1(\gcd\left(i, j\right))
\left[ \sigma_1(\gcd\left(i, j\right)) \le a\right]
\]
为了方便处理，我们还是钦定 \(n \leq m\)。
先不考虑 \(a\) 的限制，那么我们要求的就是
\[\sum^{n}_{i=1}
\sum^{m}_{j=1}
\sigma(\gcd\left(i, j\right))
\]
可以转化为
\[\sum^{n}_{d=1}
\sum^{n}_{i=1}
\sum^{m}_{j=1}
\sigma(d)
\left[ \gcd\left(i, j\right)=d \right]
\]
即
\[\sum^{n}_{d=1}
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}
\sigma(d)
\left[ \gcd\left(i, j\right)=1 \right]
\]
利用常规的变换技巧
\[\sum^{n}_{d=1}
\sigma(d)
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}
\sum_{x \mid \gcd(i,j)}\mu(x)
\]
换成整除分块形式
\[\sum^{n}_{d=1}
\sigma(d)
\sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{x=1}
\mu(x)
{\left\lfloor \frac{n}{dx}\right\rfloor}
{\left\lfloor \frac{m}{dx}\right\rfloor}
\]
和上一道题类似，我们令 \(T=dx\)，然后更改枚举顺序，有
\[\sum^{n}_{T=1}
{\left\lfloor \frac{n}{T}\right\rfloor}
{\left\lfloor \frac{m}{T}\right\rfloor}
\sum_{d \mid T}
\sigma(d)
\mu( \frac{T}{d})
\]
这道题没有 \(a\) 的限制的部分就做完了，但还要考虑我们 \(10^9\) 的 \(a\)。
我们设 \(g(T)=\sum_{d \mid T}\sigma(d)\mu( \frac{T}{d})\)，我们发现，只有在 \(\sigma(d) \leq a\) 时，才会对 \(g(T)\) 产生贡献。
那我们可以离线，把询问按照 \(a\) 从小到大排序，在从小到大逐个枚举询问的时候，\(a\) 不断变大，使得一些 \(\sigma(d)\) 对 \(g(T)\) 产生贡献。
我们枚举倍数来找到所有的 \(T\)，需要动态修改 \(g(T)\) 且区间查询，考虑树状数组。
Code
#include <bits/stdc++.h>
using namespace std;
const int M = 1e5;
const int N = 100500;
const long long Mod = 1 << 31;
int T;
bool p[N];
int prime[N],cnt;
int num[N];
int sigma[N],mu[N];
unsigned long long ans[N];
struct Query{
int id;
int n,m,a;
bool operator < (const Query &t) const {
return a < t.a;
}
}q[N];
struct Num{
int x;
bool operator < (const Num &t) const {
return sigma[x] < sigma[t.x];
}
}a[N];
void Seive() {
sigma[1] = mu[1] = 1;
for(int i = 2;i <= M; i++) {
if(!p[i]) {
cnt ++;
prime[cnt] = i;
mu[i] = -1;
sigma[i] = num[i] = i + 1;
}
for(int j = 1;j <= cnt && i * prime[j] <= M; j++) {
p[i * prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i * prime[j]] = 0;
num[i * prime[j]] = num[i] * prime[j] + 1;
sigma[i * prime[j]] = sigma[i] / num[i] * num[i * prime[j]];
break;
}
mu[i * prime[j]] = - mu[i];
num[i * prime[j]] = prime[j] + 1;
sigma[i * prime[j]] = sigma[i] * (prime[j] + 1);
}
}
for(int i = 1;i <= M; i++)
a[i].x = i;
sort(a + 1,a + M + 1);
return ;
}
int bit[N];
class BIT{
private:
int lowbit(int x) {
return x & (-x);
}
public:
void Add(int x,int y) {
for(int i = x;i <= M; i += lowbit(i))
bit[i] = (1ll * bit[i] + y) % Mod;
return ;
}
int Query(int x) {
int ans = 0;
for(int i = x;i >= 1; i -= lowbit(i))
ans = (1ll * ans + bit[i]) % Mod;
return ans;
}
}tree;
void Insert(int x) {
for(int k = 1;x * k <= M; k++)
tree.Add(x * k,1ll * mu[k] * sigma[x] % Mod);
// 预处理
}
int query(int n,int m) {
long long res = 0,last = 0,cur;
for(int l = 1,r = 0;l <= n;l = r + 1,last = cur) {
r = min(n / (n / l),m / (m / l));
cur = tree.Query(r);
res = (res + 1ll * (1ll * cur - last) * (n / l) % Mod * (m / l) % Mod) % Mod;
}
return res % Mod;
}
signed main() {
ios::sync_with_stdio(false);
cin >> T;
Seive();
for(int i = 1;i <= T; i++) {
cin >> q[i].n >> q[i].m >> q[i].a;
q[i].id = i;
if(q[i].n > q[i].m)
swap(q[i].n,q[i].m);
}
sort(q + 1,q + T + 1);
int j = 1;
for(int i = 1;i <= T; i++) {
for( ;j <= M && sigma[a[j].x] <= q[i].a; j++) {
Insert(a[j].x);
}
ans[q[i].id] = query(q[i].n,q[i].m);
}
for(int i = 1;i <= T; i++)
cout << (ans[i] % Mod + Mod) % Mod << "\n";
return 0;
}
BZOJ3309 DZY Loves Math
题目大意
对于正整数 \(n\)，定义 \(f(n)\) 为 \(n\) 所包含质因子的最大幂指数。例如 \(f(1960)=f(2^3 \times 5^1 \times 7^2)=3\)，\(f(10007)=1\)，\(f(1)=0\)。
给定正整数 \(a,b\)，求下式的值：
\[\sum^{a}_{i=1}\sum^{b}_{j=1}f(\gcd(i,j))
\]
推导
首先记 \(n=\min(a,b),m=\max(a,b)\)。
\[\begin{aligned}
& \ \ \ \ \ \ \sum^{a}_{i=1}\sum^{b}_{j=1}f(\gcd(i,j)) \\
&= \sum^{n}_{d=1}f(d)\sum^{n}_{i=1}\sum^{m}_{j=1}\left[ \gcd(i,j)=d \right] \\
&= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\left[ \gcd(i,j)=1 \right] \\
&= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\sum_{t \mid \gcd(i,j)}\mu(t) \\
&= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\sum_{t \mid i \land t \mid j}\mu(t) \\
&= \sum^{n}_{d=1}f(d)
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t =1}\mu(t)
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1} \left[ t \mid i \right]
\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1} \left[ t \mid j \right] \\
&= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}\mu(t) \left\lfloor \frac{n}{dt} \right\rfloor \left\lfloor \frac{m}{dt} \right\rfloor \\
\end{aligned}
\]
设 \(T=dt\)（十分常用的技巧），那么有
\[= \sum^{n}_{T=1} \left\lfloor \frac{n}{T} \right\rfloor \left\lfloor \frac{m}{T} \right\rfloor \sum_{d \mid T}\mu(\frac{T}{d})f(d)
\]
记 \(h=f*\mu\)，那么有
\[= \sum^{n}_{T=1} \left\lfloor \frac{n}{T} \right\rfloor \left\lfloor \frac{m}{T} \right\rfloor h(T)
\]
那么，现在的问题是如何求 \(h\)？考虑求 \(h(n)\)。
我们可以把 \(n\) 进行质因数分解：
\[n=\prod^{m}_{i=1}p_i^{c_i}
\]
考虑 \(h(T)\) 的原本写法：
\[h(T)=\sum_{d \mid T}\mu(\frac{T}{d})f(d)
\]
记 \(T=\prod^{m}_{i=1}p_i^{a_i}\)，\(d=\prod^{m}_{i=1}p_i^{b_i}\)。
考虑 \(\mu(n)\) 的定义，当 \(n\) 中含有平方质因子时 \(\mu(n)=0\)，即不会在 \(h\) 中产生贡献，可以不考虑。
所以 \(\dfrac{T}{d}\) 的质因数要满足最大幂指数小于 \(2\)，即 \(a_i-b_i \in \left\{ 0,1 \right\}\)。
所以将 \(h(T)\) 改写为如下形式时
\[\displaystyle h(T) = \sum_{ab = T} f(a) \mu(b)
\]
有
\[b=\prod^{m}_{i=1}p_i^{d_i}(d_i \in \left\{ 0,1 \right\})
\]
记 \(l=\max^m\limits_{i=1}c_i,k=\sum^{m}\limits_{i=1} \left[c_i=l \right]\)。即 \(l\) 为 \(n\) 所包含质因子的最大幂指数，\(k\) 为 \(n\) 所包含质因子幂指数中最大幂指数的个数。
我们发现 \(f(a)\) 的取值只有 \(l\) 和 \(l-1\) 两种可能（\(b\) 可能把最大幂指数的都取走，导致 \(f(a)\) 的取值少了 \(1\)）。
我们先按 \(k=m\) 和 \(k \neq m\) 两种情况分类讨论，在每一项讨论中，我们再分 \(f(a)=l\) 和 \(f(a)=l-1\) 两种子情况。
当 \(k=m\) 时，贡献为（加号前为 \(l\) 的情况，加号后为 \(l-1\) 的情况）
\[\sum^{m-1}_{s=0}(l) \times (-1)^s \times {m \choose s} + (l-1) \times (-1)^m \times {m \choose m}
\]
\(f(a)=l\) 的情况不会把最大幂指数的质因数都取走。
所以我们枚举到 \(m-1\)，中间的 \((-1)^s\) 和 \((-1)^m\) 实质上就是莫比乌斯函数，最右边的二项式系数是枚举选取的方案。
将上面的式子加号右边的部分拆开，把带 \(l\) 的部分合并到左面，得到
\[\sum^{m}_{s=0}(l) \times (-1)^s \times {m \choose s} + 1 \times (-1)^m
\]
二项式定理，得
\[-1 \times (-1)^m=(-1)^{m+1}
\]
当 \(k \neq m\) 时，
当 \(f(a)=l\) 时，设在 \(k\) 个 \(c_i=l\) 的质数中选了 \(t\) 个，另外 \(m-k\) 个质数中选了 \(s\) 个，贡献为：
\[\sum^{m-k}_{s=0}\sum^{k-1}_{t=0} {k \choose t} \times l \times (-1)^{s+t} \times {m-k \choose s}
\]
中间的 \((-1)^{s+t}\) 仍是莫比乌斯函数。因为 \(f(a)=l\)，所以 \(k\) 个质数不能被选完，因此枚举到 \(k-1\)。
左右拆开，分别二项式定理，得
\[\begin{aligned}
\sum^{k-1}_{t=0} {k \choose t} \times l \times (-1)^t \times \sum^{m-k}_{s=0}(-1)^s \times 1 ^{m-k-s} \times {m - k \choose s}
=0
\end{aligned}
\]
当 \(f(a)=l-1\) 时，\(k\) 个 \(c_i=l\) 的质数一定全部被选择，设在另外 \(m-k\) 个质数中选择 \(s\) 个，贡献为：
\[\begin{aligned}
& \ \ \ \ \ \sum^{m-k}_{s=0}(l-1) \times (-1)^k \times (-1)^s \times {m-k \choose s} \\
&= \sum^{m-k}_{s=0}(l-1) \times (-1)^k \times (-1)^s \times 1^{m-k-s} \times {m-k \choose s} \\
&= (l-1) \times (-1)^k \times \sum^{m-k}_{s=0} (-1)^s \times 1^{m-k-s} \times {m-k \choose s} \\
&= 0
\end{aligned}
\]
所以最终 \(h(n)\) 的表达式为
\[h(n)=
\begin{cases}
(-1)^{m + 1} & k = m \\
0 & k \neq m
\end{cases}
\]
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1e7;
const int M = 10500;
int T,n,m;
bool p[N + M];
int pri[N >> 2],cnt,low[N + M];
int h[N + M],a[N + M];
void Seive() {
p[1] = 1;
for(int i = 2;i <= N; i++) {
if(!p[i]) {
cnt ++;
low[i] = pri[cnt] = i;
a[i] = h[i] = 1;
}
for(int j = 1;j <= cnt && i * pri[j] <= N; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0) {
a[i * pri[j]] = a[i] + 1;
low[i * pri[j]] = low[i] * pri[j];
if(i == low[i])
h[i * pri[j]] = 1;
else if(a[i / low[i]] == a[i * pri[j]])
h[i * pri[j]] = -h[i / low[i]];
else
h[i * pri[j]] = 0;
break;
}
a[i * pri[j]] = 1;
low[i * pri[j]] = pri[j];
if(a[i] == 1)
h[i * pri[j]] = -h[i];
else
h[i * pri[j]] = 0;
}
}
for(int i = 1;i <= N; i++)
h[i] += h[i - 1];
return ;
}
signed main() {
ios::sync_with_stdio(false);
Seive();
cin >> T;
while(T --> 0) {
cin >> n >> m;
if(n > m)
swap(n,m);
int l = 1,r = 0,ans = 0;
while(l <= n) {
r = min(n / (n / l),m / (m / l));
ans += (n / l) * (m / l) * (h[r] - h[l - 1]);
l = r + 1;
}
cout << ans << "\n";
}
return 0;
}
SDOI2015 约数个数和
题目大意
设 \(d(x)\) 为 \(x\) 的约数个数（即 \(\sigma_0(x)\)），给定 \(n,m\)，求
\[\sum_{i=1}^n\sum_{j=1}^md(ij)
\]
推导
先给出一个 \(d(x)\) 的性质
\[d(ij)=\sum_{x \mid i}\sum_{y \mid j} \left[ \gcd(x,y)=1 \right] \]
证明
设 $i$ 中有因子 $p^a$，$j$ 中有因子 $p^b$，$ij$ 的因子 $k$ 中有因子 $p^c$。显然 \(ij\) 的每个因数对答案的贡献都是 \(1\)，我们考虑构建一种策略使得 \(ij\) 的因数与数对 \((x,y)\) 构成双射。
设数对为 \((x,y)\)。
- 如果 \(c \leq a\)，我们只在 \(i\) 中选择，即 \(x=c,y=0\)；
- 否则，我们在 \(i\) 中选择 \(p^a\)，在 \(j\) 中选择 \(p^{c-a}\)（表示为 \(p^{a+y}\)）。
当 \(c > a\) 时，我们把取满的部分赋为 \(0\)，即 \(x=0,y=c-a\)。
显然，\(x,y\) 中有且仅有一个数为 \(0\)。
那么枚举 \(ij\) 的约数，就转化为了枚举我们选择策略的数对 \((x,y)\)。
\[\gcd(p^x,p^y)=p^{\min(x,y)}=p^0=1
\]
所以
\[\gcd(x,y)=1 \Longleftrightarrow x = 0 \lor y = 0
\]
即满足我们的选择策略。
这样就不重不漏地枚举了 \(ij\) 的每一个约数，故等式成立。
那么利用上面的式子转化式子。
\[\sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{y \mid j} \left[ \gcd(x,y)=1 \right]
\]
把 \(x,y\) 提前，
\[\sum^{n}_{x=1}
\sum^{m}_{y=1}
\left\lfloor \dfrac{n}{x} \right\rfloor
\left\lfloor \dfrac{m}{y} \right\rfloor
\left[ \gcd(x,y)=1 \right]
\]
把 \(x,y\) 换成 \(i,j\)，
\[\sum^{n}_{i=1}
\sum^{m}_{j=1}
\left\lfloor \frac{n}{i} \right\rfloor
\left\lfloor \frac{m}{j} \right\rfloor
\sum_{d \mid \gcd(i,j)} \mu(d)\]
把 \(d\) 提前，
\[\sum^{n}_{d=1}
\mu(d)
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}
\left\lfloor \frac{n}{di} \right\rfloor
\left\lfloor \frac{m}{dj} \right\rfloor
\]
设函数 \(g(n)=\sum^{n}\limits_{i=1}\left\lfloor \frac{n}{i} \right\rfloor\)，那么有
\[\sum^{n}_{d=1}
\mu(d)
g(\left\lfloor \frac{n}{d} \right\rfloor)
g(\left\lfloor \frac{m}{d} \right\rfloor)
\]
显然，\(g(n)=\sum^{n}\limits_{i=1} d(i)\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int M = 50050;
int T,n,m;
int d[100500],num[100500],mu[100500],g[100500];
int pri[100500],cnt = 0;
bool p[100500];
void Seive() {
d[1] = 1;
mu[1] = 1;
for(int i = 2;i <= M; i++) {
if(!p[i]) {
cnt ++;
pri[cnt] = i;
d[i] = 2;
num[i] = 1;
mu[i] = -1;
}
for(int j = 1;j <= cnt && i * pri[j] <= M; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0) {
num[i * pri[j]] = num[i] + 1;
d[i * pri[j]] = d[i] / num[i * pri[j]] * (num[i * pri[j]] + 1);
break;
}
num[i * pri[j]] = 1;
d[i * pri[j]] = d[i] * 2;
mu[i * pri[j]] = -mu[i];
}
}
for(int i = 1;i <= M; i++) {
g[i] = g[i - 1] + d[i];
mu[i] += mu[i - 1];
}
return ;
}
signed main() {
ios::sync_with_stdio(false);
Seive();
cin >> T;
while(T --> 0) {
cin >> n >> m;
if(n > m)
swap(n,m);
int ans = 0;
for(int l = 1,r = 0;l <= n; l = r + 1) {
r = min(n / (n / l),m / (m / l));
ans += (mu[r] - mu[l - 1]) * g[n / l] * g[m / l];
}
cout << ans << "\n";
}
return 0;
}
SDOI2017 数字表格
题目大意
记 \(f(i)\) 表示斐波那契数列的第 \(i\) 项，求
\[\prod^{n}_{i=1}\prod^{m}_{j=1}f(\gcd(i,j))
\]
答案对 \(10^9+7\) 取模。
推导
\[\prod^{n}_{i=1}\prod^{m}_{j=1}f(\gcd(i,j))
\]
首先是经典枚举技巧，
\[\begin{aligned}
& \ \ \ \ \
\prod^{n}_{d=1}f(d)
\prod^{n}_{i=1}
\prod^{m}_{j=1}
\left[ \gcd(i,j)=d \right] \\
&=
\prod^{n}_{d=1}f(d)^{\sum^{n}\limits_{i=1} \sum^{m}\limits_{j=1}\left[ \gcd(i,j)=d \right]}
\end{aligned}
\]
指数部分单独取出来处理，直接跳步，得到
\[\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}
\mu(t)
\left\lfloor \frac{n}{dt} \right\rfloor
\left\lfloor \frac{m}{dt} \right\rfloor
\]
代回原式，
\[\prod^{n}_{d=1}f(d)^
{
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}
\mu(t)
\left\lfloor \frac{n}{dt} \right\rfloor
\left\lfloor \frac{m}{dt} \right\rfloor
}
\]
老套路，设 \(T=dt\)，
\[\prod^{n}_{T=1}\prod_{d \mid T} f(T)^
{
\mu(\frac{T}{d})
\left\lfloor \frac{n}{T} \right\rfloor
\left\lfloor \frac{m}{T} \right\rfloor
}
\]
把中间一部分提出来单独处理，
设
\[g(n)=\prod_{d \mid n} f(d)^{\mu(\frac{n}{d})}
\]
代入原式，就得到了
\[\prod^{n}_{T=1}g(T)^
{
\left\lfloor \frac{n}{T} \right\rfloor
\left\lfloor \frac{m}{T} \right\rfloor
}
\]
\(g(T)\) 直接暴力求。
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
int T,n,m;
const int N = 1005000;
const int M = 1000000;
const int Mod = 1e9 + 7;
bool p[N];
int pri[N],cnt;
int mu[N],inv[N];
int f[N],g[N];
long long Pow(long long a ,long long b) {
long long ans = 1 ;
long long base = a % Mod;
while( b > 0 ) {
if( b&1 != 0 )
ans = ans * base % Mod;
base = base * base % Mod;
b = b >> 1;
}
return ans;
}
void Seive() {
mu[1] = p[1] = f[1] = inv[1] = 1;
for(int i = 2;i <= M; i++) {
f[i] = (f[i - 1] + f[i - 2]) % Mod;
inv[i] = Pow(f[i],Mod - 2) % Mod;
if(!p[i]) {
mu[i] = -1;
cnt ++;
pri[cnt] = i;
}
for(int j = 1;j <= cnt && i * pri[j] <= M; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
}
mu[i * pri[j]] = -mu[i];
}
}
for(int i = 0;i <= M; i++)
g[i] = 1;
for(int i = 1;i <= M; i++) {
if(mu[i] == 0)
continue;
for(int j = i;j <= M; j += i) {
if(mu[i] == 1)
g[j] = g[j] * f[j / i] % Mod;
else
g[j] = g[j] * inv[j / i] % Mod;
}
}
for(int i = 2;i <= M; i++)
g[i] = g[i] * g[i - 1] % Mod;
return ;
}
signed main() {
ios::sync_with_stdio(false);
Seive();
cin >> T;
while(T --> 0) {
cin >> n >> m;
if(n > m)
swap(n,m);
int ans = 1,Inv = 0;
for(int l = 1,r = 0;l <= n; l = r + 1) {
// cerr << "a";
r = min(n / (n / l),m / (m / l));
Inv = g[r] * Pow(g[l - 1],Mod - 2) % Mod;
ans = ans * Pow(Inv,(n / l) * (m / l) % (Mod - 1)) % Mod;
}
ans = (ans + Mod) % Mod;
cout << ans << "\n";
}
return 0;
}
BZOJ4407 于神之怒加强版
题目大意
给定 \(n,m,k\)，计算
\[\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k
\]
对 \(10^9 + 7\) 取模的结果。
推导
经典枚举
\[\begin{aligned}
& \ \ \ \ \
\sum^{n}_{d=1}d^k
\sum^{n}_{i=1}
\sum^{m}_{j=1}
\left[ \gcd(i,j)=d \right] \\
&=
\sum^{n}_{d=1}d^k
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}
\sum_{t \mid \gcd(i,j)}\mu(t) \\
&=
\sum^{n}_{d=1}d^k
\left\lfloor \dfrac{n}{dt} \right\rfloor
\left\lfloor \dfrac{m}{dt} \right\rfloor
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}\mu(t) \\
&=
\sum^{n}_{T=1}
\sum_{d \mid T}d^k \mu(\frac{T}{d})
\left\lfloor \dfrac{n}{T} \right\rfloor
\left\lfloor \dfrac{m}{T} \right\rfloor
\end{aligned}
\]
令 \(g(n)=\sum_{d \mid n}d^k\mu(\frac{n}{d})\)，即 \(g(n)= \operatorname{id}_k * \mu\)，显然这是一个积性函数。
当 \(p\) 为质数时，\(d \in \left\{ 1,p \right\}\)，容易得到 \(g(p)=p^k-1\)。
考虑 \(p\) 为质数时，\(g(p^c)\) 的值：
\[\begin{aligned}
& \ \ \ \ \
\sum_{d \mid p^c}d^k\mu(\frac{p^c}{d}) \\
&=
\sum^{c}_{i=0}(p^i)^k\mu(\frac{p^c}{p^i}) \\
&=
\sum^{c}_{i=0}(p^i)^k\mu(p^{c-i})
\end{aligned}
\]
当 \(c-i>1\) 时，有 \(\mu(p^{c-i})=0\)，不产生贡献。产生贡献的只有 \(i=c\) 和 \(i=c-1\) 两种情况。
容易得到 \(g(p^c)=p^{kc}-p^{kc-k}\)。
考虑 \(\gcd(p,y)=1\)，\(g(p^cy)\) 的值，根据积性函数的性质即可。
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
int T,k,n,m;
const int N = 5005000;
const int M = 5000000;
const int Mod = 1e9 + 7;
bool p[N];
int pri[N],cnt;
int g[N],sum[N];
long long Pow(long long a ,long long b) {
long long ans = 1 ;
long long base = a % Mod;
b = b % (Mod - 1);
while(b) {
if(b & 1)
ans = ans * base % Mod;
base = base * base % Mod;
b >>= 1;
}
return ans;
}
void Seive() {
memset(g,1,sizeof(g));
g[1] = 1;
for(int i = 2;i <= M; i++) {
if(!p[i]) {
cnt ++;
pri[cnt] = i;
g[i] = (Pow(i,k) - 1 + Mod) % Mod;
}
for(int j = 1;j <= cnt && i * pri[j] <= M; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0) {
g[i * pri[j]] = g[i] * (g[pri[j]] + 1) % Mod;
break;
}
g[i * pri[j]] = g[i] * g[pri[j]] % Mod;
}
}
for(int i = 1;i <= M; i++)
sum[i] = (sum[i - 1] + g[i]) % Mod;
return ;
}
signed main() {
ios::sync_with_stdio(false);
cin >> T >> k;
Seive();
while(T --> 0) {
cin >> n >> m;
if(n > m)
swap(n,m);
int ans = 0;
for(int l = 1,r = 0;l <= n; l = r + 1) {
r = min(n / (n / l),m / (m / l));
ans = (ans + (sum[r] - sum[l - 1]) % Mod * (n / l) % Mod * (m / l)) % Mod;
ans = (ans + Mod) % Mod;
}
ans = (ans + Mod) % Mod;
cout << ans << "\n";
}
return 0;
}
Crash 的数字表格 / JZPTAB
题目大意
给定 \(n,m\)，求
\[\sum^{n}_{i=1}\sum^{m}_{j=1}\operatorname{lcm}(i,j)
\]
答案对 \(20101009\) 取模。
推导
\[\sum^{n}_{i=1}\sum^{m}_{j=1}\frac{i \times j}{\gcd(i,j)}
\]
经典枚举
\[\sum^{n}_{d=1}\sum^{n}_{i=1}\sum^{m}_{j=1}\frac{i \times j}{d} \left[ \gcd(i,j)=d \right]
\]
提出 \(d\)，\(i,j\) 都除以了 \(d\)，多除的一个放到前面去，
\[\sum^{n}_{d=1}d
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}
i \times j \left[ \gcd(i,j)=1 \right]
\]
然后得
\[\sum^{n}_{d=1}d
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}
\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}
i \times j
\sum_{t \mid \gcd(i,j)}\mu(t)
\]
设 \(g(n)=\frac{n \times(n + 1)}{2}\)，
\[\sum^{n}_{d=1}d
\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}t^2 \cdot \mu(t)
\times
g(\left\lfloor \frac{n}{dt} \right\rfloor)
\times
g(\left\lfloor \frac{m}{dt} \right\rfloor)
\]
设 \(T=dt\)，有
\[\sum^{n}_{T=1}g(\left\lfloor \frac{n}{dt} \right\rfloor)
\times
g(\left\lfloor \frac{m}{dt} \right\rfloor)
\times
T
\times
\sum_{d \mid T}d \cdot \mu(d)
\]
设 \(f(n)=\sum_{d \mid n}d \cdot \mu(d)\)，显然是积性函数。
当 \(p\) 为质数时，有 \(f(p)=1-p,f(p^k)=f(p)\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 10005000;
const int M = 10000000;
const int Mod = 20101009;
int n,m;
bool p[N];
int cnt,pri[N >> 4];
int f[N];
void Seive() {
f[1] = 1;
for(int i = 2;i <= M; i++) {
if(!p[i]) {
cnt ++;
pri[cnt] = i;
f[i] = (1 - i) % Mod;
}
for(int j = 1;j <= cnt && i * pri[j] <= M; j++) {
p[i * pri[j]] = 1;
if(i % pri[j] == 0) {
f[i * pri[j]] = f[i];
break;
}
f[i * pri[j]] = f[i] * f[pri[j]] % Mod;
}
}
for(int i = 2;i <= M; i++)
f[i] = (f[i] * i % Mod + f[i - 1]) % Mod;
return ;
}
int g(int x) {
return (x * (x + 1) / 2) % Mod;
}
signed main() {
ios::sync_with_stdio(false);
Seive();
cin >> n >> m;
if(n > m)
swap(n,m);
int ans = 0;
for(int l = 1,r = 0;l <= n; l = r + 1) {
r = min(n / (n / l),m / (m / l));
ans = (ans + g(n / l) * g(m / l) % Mod * (f[r] - f[l - 1] + Mod) % Mod + Mod) % Mod;
}
cout << ans << "\n";
return 0;
}
\[\text{END}
\]
