前置知识
莫比乌斯函数
定义
设 \(n=\prod^{k}_{i=1}p_i^{c_i}\),其中 \(p_i\) 为 \(n\) 的质因子,\(c_i \geq 1\)。
性质
莫比乌斯函数是积性函数,有以下性质
常用于莫比乌斯反演中。
证明
- 当 \(n=1\) 时显然成立;
- 当 \(n \neq 1\) 时,设
那么有
考虑我们选 \(i\) 个本质不同的质因子,有 \(\binom{k}{i}\) 种方法,其中每个数的莫比乌斯函数值为 \((-1)^i\),所以这样计算出原式的贡献。
那么我们需要证明
考虑二项式定理
证毕。
莫比乌斯反演
反演公式
设 \(f(n),g(n)\) 是定义在正整数集上的两个函数。
形式 \(1\)
若有
则有
证明
考虑对下式进行变换
代入定义,有
可以交换求和号
详细解释
对于上面这一步变换,有以下考虑。由于 \(d \mid n\),存在 \(t\) 满足 \(t \mid \!\frac{n}{d}\)。那么对于每一个符合条件的 \(t\) 都有与之对应的 \(d\),那么就可以转化为
考虑原式,相当于先枚举符合条件的 \(d\),再枚举在此时 \(d\) 基础上符合条件的 \(t\),我们发现,换成现在的式子后,对答案的统计仍然不重不漏。
还可以换一种方式考虑,对于原式,我们可以转化为
显然,对于两个限制条件,我们可以直接消去第一个条件,再将艾弗森括号转换进和式,得到
得到这个式子,我们就可以直接交换 \(\mu(d)\) 和 \(g(t)\) 的位置了,并可以将和式中的一个条件拆开成为两个和式
又由 \(\sum_{d\mid n}\mu(d)=\left[ n=1 \right]\) 可以得到
把艾弗森括号里的部分转化一下
只有当 \(n=t\) 时,该式子才有值,即
证毕。
形式 \(2\)
若有
则有
证明与形式 \(1\) 的证明类似,不再展开说了。
例题
POI2007 ZAP-Queries
题目大意
给定 \(n\) 和 \(m\),求
证明
为了方便处理,我们钦定 \(n \leq m\)。
首先考虑 \(d\) 同时是 \(i\) 和 \(j\) 的最大公约数,所以满足 \(d \ | \ i\) 且 \(d \ | \ j\)。
于是我们可以把 \(d\) 提出来。
然后我们可以利用莫比乌斯函数的结论,变为
此时,实际上 \(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\) 的和式后面。我们也可以顺便把其他和式的限制条件放在她后面。
然后我们可以发现,观察后面两个和式。例如第二个和式,满足 \(\left[t \ | \ i \right]\) 的 \(i\) 最多有 \(\left\lfloor\\\frac{n}{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\),求
即求最大公约数为质数的数对个数。
思路
为了方便处理,我们还是钦定 \(n \leq m\)。
我们可以枚举一个 \(d\),把原式转化为
利用莫比乌斯函数性质,有
交换求和号,令 \(i=\frac{i}{t},j=\frac{j}{t}\),有
式子推到这里,如果我们去枚举 \(d\),肯定会 TLE,我们令 \(T = dt\),有
我们发现,后面那部分可以预处理。
考虑每个质数 \(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\),求
为了方便处理,我们还是钦定 \(n \leq m\)。
先不考虑 \(a\) 的限制,那么我们要求的就是
可以转化为
即
利用常规的变换技巧
换成整除分块形式
和上一道题类似,我们令 \(T=dx\),然后更改枚举顺序,有
这道题没有 \(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;
}
未完待续 \(\dots\)