缘起
数论千万条, 反演第一条, 反演不会做, 队友两行泪... 所以我们来学一下莫比乌斯反演吧~
分析
《莫比乌斯函数入门》中我们已经学习了mobius函数. 现在要来进一步学习mobius反演.
首先思考一个有趣的问题
代码语言：javascript复制给定 N 和 M, 求满足 1<=x<=N, 1<=y<=M 且 gcd(x, y) 为质数的点对(x,y) 的数量.
数据范围 1<=N,M<=1e6显然，暴力的 O(NM) 算法是绝对不能接受的.
此题的前置技能就是mobius反演.
何谓反演? 圆的反演
首先，我们先理解一下什么叫做反演? 反演的"演"的意思是演绎（deduction） 的意思. 初中学习平面几何的时候就学过圆的反演公式. 如下图所示, 如果 , 其中 r 是圆的半径. 那么 称为 关于圆的反演点.
所以我们就知道什么叫做反演了——知道了 的位置, 我们可以很轻松的定位 的位置, 反之亦然.
反演是一种逻辑演绎的技巧, 如果一个问题（例如上面的 ）比较难处理，但是它的反演问题（例如上面的 ）却比较好处理, 则可以转换为其反演问题进行处理.
讲清楚了什么是反演，我们就来讲什么是mobius反演，为了讲清楚mobius反演就不得不提及dirichlet卷积.
Dirichlet 卷积
dirichlet卷积定义如下
显然，dirichlet卷积满足交换律和结合律
关于dirichlet卷积小例子如下
分别是不变映射和恒等映射, 那么
分别表示因子个数函数和因子和函数.
既然dirichlet卷积是一个运算，那么我们自然关心它的单位元 ，正如加法运算 的 0 和乘法运算的1一样重要.
dirichlet卷积的单位元是
或者记做 , 其中 [] 我们称之为示性函数（数学中一般记做 ）.
这非常好证明. 只需要证明 有 即可
既然讲到了dirichlet卷积的单位元, 所以不得不提及dirichlet卷积的逆元. 事实上，我们一般关心mobius函数 和欧拉函数 的逆元. 而 nonvariant 函数就是 的逆元. 也即
Proof: 按照定义, 只需要证明 , 我们有
而 假设 n 有 k（k > 0） 个不同的素因子. 则
证毕
至于欧拉函数，我们有
Proof: 和证明欧拉函数的计算公式类似. 令
我们证明 f 是一个积性函数. 事实上我们有
有了此积性公式，剩下的事情就很简单了. 令 ， 则
其中用到了
所以我们就得到了
证毕.
一个值得注意的事实是, 上述欧拉函数的式子两边同时卷积上 , 我们有
展开写我们就得到了 和 这两个重要的积性函数之间的关系如下
好了, 介绍完毕dirichlet卷积，我们开始介绍mobius反演.
所谓mobius反演, 就是如下两个反演推导式子（回想一下本文开头对反演一词的解释）
站在反演的角度， g 是比较容易解决的问题, f 是比较难解决的问题.
上面两个式子的
- 相同点都是将比较难解决的 f 用比较容易解决的 g 表示了出来.
- 不同点是第一个式子是约数关系，第二个式子是倍数关系.
Proof:
第一个式子的证明
站在dirichlet卷积的角度, , 所以
其中， 我们要注意到dirichlet卷积的交换律和结合律即可.
第二个式子的证明
令 则考虑到 的逆元是 nonvariant, 我们有
证毕.
至此，mobius反演讲解完毕. 我们看看如何使用mobius反演解决本文开头提出的问题.
首先，我们考虑如下定义
那么，显然的, 我们有
而且，幸运的是 , g 是一个非常好解决的问题
所以，我们利用mobius反演就可以计算 f 了
于是答案为
令
则答案进一步化简为
如果我们能 O(N) 预处理出 sum(d) 的话, 则本题的复杂度就是 O(N), 这是完全可以接受的.
事实上, 我们有如下预处理 sum 的伪代码
代码语言：javascript复制for (re i = 1; i <= tot; i ) // prime[1,...,tot] 是所有质数, m是考察上限
{
for (re j = 1, t; (t = prime[i] * j) <= m; j )
{
sum[t] = mu[j];
}
}于是就开心的解决了此问题了, 复杂度是 O(N) 的.
奉上完整的代码
代码语言：javascript复制// 此处省略一些头部文件...
const int maxn = 1e6 5;
int n, m, prime[maxn], isp[maxn], tot, mu[maxn], sum[maxn];
ilv mobius(int m)
{
mu[1] = 1;
for (re i = 2; i <= m; i )
{
if (!isp[i])
{
prime[ tot] = i;
mu[i] = -1;
}
for (re j = 1, t; j <= tot && (t = prime[j] * i) <= m; j )
{
isp[t] = 1;
if (i % prime[j])
{
mu[t] = -mu[i];
}
else
{
mu[t] = 0;
break;
}
}
}
}
signed main()
{
#ifdef LOCAL
freopen("d:\data.in", "r", stdin);
// freopen("d:\my.out", "w", stdout);
#endif
read(n), read(m);
if (n > m)
{
swap(n, m);
}
mobius(m);
for (re i = 1; i <= tot; i )
{
for (re j = 1, t; (t = prime[i] * j) <= m; j )
{
sum[t] = mu[j];
}
}
int ans = 0;
for (re i = 1; i <= n; i )
{
ans = (n / i) * (m / i) * sum[i];
}
write(ans);
flush();
return 0;
}应用上述程序可以快速计算出当 n = m = 1e6 时的答案为 274934018418， 也就两千多亿——要是我有这么多钱就好了~
为了证明上面给出的mobius反演的正确性. 我们来看一道mobius反演的入门题 hdu 1695 gcd
代码语言：javascript复制给定 a, b, c, d, k, 求满足 a<=x<=b, c<=y<=d, (x,y)=k 的不同点对 (x, y) 的数量.
注意, (5, 7) 和 (7, 5) 被认为是相同点对
【输入】
多样例,第一行是样例数, 每个样例 有五个数.
0<a<=b<=1e5, 0<c<=d<=1e5, 0<=k<=1e5
输入保证 a = c = 1
【输出】
输出答案
【样例输入】
2
1 3 1 5 1
1 11014 1 14409 9
【样例输出】
Case 1: 9
Case 2: 736427
【hint】
对于第一个样例, 9对(x,y) 如下
(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 5), (3, 4), (3, 5).
【限制】
3s 32MB首先，使用上面一样的推导，我们有
其中 f(1) 是集合 {} 的互质点对个数. 但是题目要求的是不同的互质点对的数量. 所以还要减去多算的点对. 而多算的点对恰好就是集合 {} 的互质点对个数的一半.
于是我们可以写出如下代码
代码语言：javascript复制//#include "stdafx.h"
//#define LOCAL
#pragma GCC optimize(2)
#pragma G optimize(2)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <string>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
using namespace std;
#define int long long
#define re register int
#define ci const int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to, len; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l hjt[cur].r >> 1)
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
namespace fastio
{
const int BUF = 1 << 21;
char fr[BUF], fw[BUF], *pr1 = fr, *pr2 = fr;int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1 ; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw[pw ] = c; }
ili read(int &x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while(!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while(isdigit(c)) x = (x << 3) (x << 1) (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double &x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) (xx << 1) (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x = (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 48); }
ilv writeln(int x) { write(x);pc(10); }
ili read(char *x)
{
char c = gc(); if (!~c) return EOF;
while(!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x = c, c = gc();
*x = 0; return 1;
}
ili readln(char *x)
{
char c = gc(); if (!~c) return EOF;
while(c == 10) c = gc();
while(c >= 32 && c <= 126) *x = c, c = gc();
*x = 0; return 1;
}
ilv write(char *x) { while(*x) pc(*x ); }
ilv writeln(char *x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
const int maxn = 1e5 5;
int a, b, c, d, k, mu[maxn], isp[maxn], prime[maxn], tot;
ilv mobius(int m)
{
mu[1] = 1;
for (re i = 2; i <= m; i )
{
if (!isp[i])
{
prime[ tot] = i;
mu[i] = -1;
}
for (re j = 1, t; j <= tot && (t = i * prime[j]) <= m; j )
{
isp[t] = 1;
if (i % prime[j])
{
mu[t] = -mu[i];
}
else
{
mu[t] = 0;
break;
}
}
}
}
ili kk(int n, int m)
{
int ans = 0;
for (re i = 1; i <= n; i )
{
ans = (n / i) * (m / i) * mu[i];
}
return ans;
}
signed main()
{
#ifdef LOCAL
freopen("d:\data.in", "r", stdin);
// freopen("d:\my.out", "w", stdout);
#endif
mobius(maxn - 1);
int kase; scanf("%lld", &kase);
for (re kse = 1; kse <= kase; kse )
{
write("Case "), write(kse), write(": ");
read(a), read(b), read(c), read(d), read(k);
if (!k)
{
writeln("0");
continue;
}
if (b > d) swap(b, d);
b /= k, d /= k;
writeln(kk(b, d) - kk(b, b) / 2);
}
flush();
return 0;
}ac情况如下
代码语言：javascript复制Accepted 93ms 5460kB
