首页 > 代码库 > M?bius反演（转）
M?bius反演（转）
一，反演理论的初识
以一个方程组的例子来引出主题：
a * x + b * y = p, c * x + d * y = q.
在已知 a, b, c, d, p, q 的前提下我们如何来求出x 和y 呢？
在计算机解决这类问题的一个最基本方法就是高斯消元法。像这样通过已知的结果来求未知数的方法就叫做反演。
用数学语言来说： 为得到某个组合计数问题的解, 我们首先设法求出相应序列f(n)所满足的(累计)关系式 ，其中g(n)是已知序列, 然后从中解出
的方法，就是反演。
二，反演算法具体过程
对于关系式
其中c 和 d 都为下三角阵，且互逆下面分别简写为C和D。
由此可见，当 g(n) 和 C矩阵已知的前提下，我们可以求出C的逆矩阵，然后根据上式即可反演出f(x) 的解。说起来肯能比较抽象，下面我们通过一个简单的例子来说明这一过程。
对于等式 f(1) + f(2) +...+ f(n) = (n + 1) * n / 2 = g(n
现在假设我们只知道g(n), 即1到n的和，我们该如何求出f(i) 呢？
显然f(i) =i 。 下面我们将用反演算法算出这个结果。
假设n = 3 根据上式我们知道C = ， 从而可知其逆矩阵D =
。
用D和函数g相乘，我们可知
f(1) = g(1) = 1, f(2) = g(2) - g(1) = 2, f(3) = g(3) - g(2) + 0 *g(1) = 3.
其实对于f(n) 有 f(n) =
满足
u(i) = 1 iff i == n ,
u(i) = -1 iff i + 1 == n ,
其他情况u(i)=0。
f(n) = 就是我们所求出来的反演公式。 总结起来说反演过程就是求逆矩阵而已，SO EASY！
上面只是热身阶段，下面开始逐步进入这篇文章的主题！
三，二项式反演公式
问题： 在n个数字1,2,…,n形成的n!个排列a1a2…an中, 满足ai != i 的排列有多少个？
这个问题可以有很多种解法，递归法和容斥关系法等等。现在我们用二项式反演法去解决这个问题！
定义g(n) 为n个数字所组成的全排列，很显然g(n) = n! 。
再定义f(i)为i个数字都满足ai != i的排列个数。
则我们有 = g(n) , (其中com为组合数)。
到这里我们就可以通过反演规则得出f(n) = ，问题就在于如何求u(i)，即D矩阵第n行第i列元素的值。答案即为f(n)的值。
其实u(i) = (-1) ^(n - i) * com(n, i).
即证明
通过上述公式我们可以反演出
四，M?bius反演公式
设f(n)和g(n)是定义在正整数集合上的两个函数.
成立当且仅当下式成立
函数m(m)的定义如下：
例如: 30=2*3*5, 121=11 * 11.
u(30)=(-1) , u(3)=-1,mu(121)=0.
其证明比较繁琐，此处略去！
莫比乌斯反演公式应用非常之广，尤其在求和gcd相关的题目时，如果想弄明白就从AC题目入手吧！
题目： spoj7001
| 7001. Visible Lattice Points Problem code: VLATTICE |
Consider a N*N*N lattice. One corner is at (0,0,0) and the opposite one is at (N,N,N). How many lattice points are visible from corner at (0,0,0) ? A point X is visible from point Y iff no other lattice point lies on the segment joining X and Y.
Input :
The first line contains the number of test cases T. The next T lines contain an interger N
Output :
Output T lines, one corresponding to each test case.
Sample Input :
3
1
2
5
Sample Output :
7
19
175
Constraints :
T <= 50
1 <= N <= 1000000
题解：即求三维坐标最大公约数为1的点。设f(x) 为公约数为x的点数，g（x）为最大公约数为x的点数，则有f(x) = g(d) ,满足d是x的倍数且d<=N.其中f（x）= (n / x + 1) ^ 3 - 1. 通过莫比乌斯反演得到
g(x) = f(d) * u(d / x), 满足d是x的倍数且d<=n,而u(d)即为莫比乌斯函数，从而可以迅速得解。
#include
#include
#include
using namespace std;
const int Maxn = 1000010;
typedef long long ll;
int prim[Maxn], top;
int miu[Maxn];
bool first[Maxn];
void getPrim()
{
miu[1] = 1;
for (int i = 2; i < Maxn; ++i)
{
if (!first[i]) miu[i] = -1, prim[top++] = i;
for (int j = 0; j < top && prim[j] * i < Maxn; ++j)
{
int x = prim[j] * i;
first[x] = 1;
if (i % prim[j] == 0) break;
miu[x] = -miu[i];
}
miu[i] += miu[i - 1];
}
}
int n;
ll solve(ll n)
{
ll ret = 0;
for (ll i = 1; i <= n;)
{
ll j = n / (n / i);
ret += (miu[j] - miu[i - 1]) * ((n / i + 1) * (n / i + 1) * (n / i + 1) - 1);
i = j + 1;
}
return ret;
}
int main()
{
int T;
scanf("%d", &T);
getPrim();
while (T--)
{
scanf("%d", &n);
printf("%lld\n", solve(n));
}
return 0;
}
spoj4491
| 4491. Primes in GCD Table Problem code: PGCD |
Johnny has created a table which encodes the results of some operation -- a function of two arguments. But instead of a boring multiplication table of the sort you learn by heart at prep-school, he has created a GCD (greatest common divisor) table! So he now has a table (of height a and width b), indexed from (1,1) to (a,b), and with the value of field (i,j) equal to gcd(i,j). He wants to know how many times he has used prime numbers when writing the table.
Input
First, t ≤ 10, the number of test cases. Each test case consists of two integers, 1 ≤ a,b < 107.
Output
For each test case write one number - the number of prime numbers Johnny wrote in that test case.
Example
Input:
2
10 10
100 100
Output:
30
2791
这道题目较难，但是AC之后会加深对莫比乌斯函数的理解。题解可以在网上搜，这里不再赘述
#include
#include
#include
#include
using namespace std;
const int Maxn = 1e7 + 100;
int prim[Maxn], top;
typedef long long ll;
ll miu[Maxn];
bool first[Maxn];
int two[Maxn];
void getPrim()
{
miu[1] = 0;
for (int i = 2; i < Maxn; ++i)
{
if (!first[i]) miu[i] = 1, prim[top++] = i;
for (int j = 0; j < top; ++j)
{
int x = prim[j] * i;
if (x >= Maxn) break;
first[x] = 1;
two[x] = two[i];
if (i % prim[j] == 0) two[x]++;
if (miu[i] == 0) miu[x] = 0;
else miu[x] = miu[i] > 0 ? -(miu[i] + 1) : -(miu[i] - 1);
if (two[x] == 1) miu[x] = (miu[x] > 0 ? 1 : -1);
else if (two[x] > 1) miu[x] = 0;
if (i % prim[j] == 0) break;
}
miu[i] += miu[i - 1];
}
}
int a, b;
struct Vector
{
ll v[Maxn];
int size;
ll& operator[](int k)
{
return v[k];
}
}vec;
ll solve(ll a, ll b)
{
ll ret = 0;
vec.size = 0;
for (ll i = 1; i <= a;)
{
vec[vec.size++] = i;
i = min(a / (a / i), b / (b / i)) + 1;
}
vec[vec.size++] = a + 1;
for (int i = 1; i < vec.size; ++i)
{
ll t = vec[i] - 1;
ll s = vec[i - 1];
ret += (miu[t] - miu[s - 1]) * (a / s) * (b / s);
}
return ret;
}
int main()
{
getPrim();
int T;
scanf("%d", &T);
while (T--)
{
scanf("%d%d", &a, &b);
if (a > b) swap(a, b);
ll ans = 0;
ans = solve(a, b);
printf("%lld\n", ans);
}
return 0;
}