码迷,mamicode.com
首页 > 其他好文 > 详细

M?bius反演(转)

时间:2014-05-13 10:04:15      阅读:318      评论:0      收藏:0      [点我收藏+]

标签:des   style   blog   code   c   tar   

一,反演理论的初识

 

以一个方程组的例子来引出主题:

a * x + b * y = p, c * x + d * y = q.

在已知 a, b, c, d, p, q 的前提下我们如何来求出呢?
在计算机解决这类问题的一个最基本方法就是高斯消元法。像这样通过已知的结果来求未知数的方法就叫做反演。

 

用数学语言来说: 为得到某个组合计数问题的解我们首先设法求出相应序列f(n)所满足的(累计)关系式 bubuko.com,布布扣 其中g(n)是已知序列然后从中解出bubuko.com,布布扣  的方法,就是反演。

 

 

二,反演算法具体过程

对于关系式 

 bubuko.com,布布扣  

bubuko.com,布布扣

其中c 和 d 都为下三角阵,且互逆下面分别简写为C和D。

由此可见,当 g(n) 和 C矩阵已知的前提下,我们可以求出C的逆矩阵,然后根据上式即可反演出f(x) 的解。说起来肯能比较抽象,下面我们通过一个简单的例子来说明这一过程。
对于等式  f(1) + f(2) +...+ f(n) =  (n + 1) * n / 2 = g(n

现在假设我们只知道g(n), 1n的和,我们该如何求出f(i) 呢?

显然f(i) =i 。 下面我们将用反演算法算出这个结果。

假设n = 3 根据上式我们知道C =  bubuko.com,布布扣, 从而可知其逆矩阵D = bubuko.com,布布扣

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) =   bubuko.com,布布扣

 

 满足 

u(i) = 1 iff  i == n , 

u(i) = -1 iff  i + 1 == n , 

其他情况u(i)=0

 

f(n) = bubuko.com,布布扣就是我们所求出来的反演公式。 总结起来说反演过程就是求逆矩阵而已,SO EASY

上面只是热身阶段,下面开始逐步进入这篇文章的主题!

 

三,二项式反演公式

问题: n个数字1,2,,n形成的n!个排列a1a2an满足ai != 的排列有多少个?

 

这个问题可以有很多种解法,递归法和容斥关系法等等。现在我们用二项式反演法去解决这个问题!

定义g(n) n个数字所组成的全排列,很显然g(n) = n! 

再定义f(i)i个数字都满足ai != i的排列个数。

则我们有 bubuko.com,布布扣 = g(n) , (其中com为组合数)

到这里我们就可以通过反演规则得出f(n) = bubuko.com,布布扣,问题就在于如何求u(i),即D矩阵第n行第i列元素的值。答案即为f(n)的值。

其实u(i) = (-1) ^(n - i) * com(n, i). 

即证明 

 bubuko.com,布布扣

 

 

bubuko.com,布布扣

通过上述公式我们可以反演出

 

 

 

 

 

 

bubuko.com,布布扣

 

 

四,M?bius反演公式

 

f(n)g(n)是定义在正整数集合上的两个函数.

bubuko.com,布布扣

 成立当且仅当下式成立

bubuko.com,布布扣

 

函数m(m)的定义如下:

 

bubuko.com,布布扣

 

例如: 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
#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 = 2i < Maxn++i)
   {
       if (!first[i]) miu[i] = -1, prim[top++] = i;
       for (int j = 0j < top && prim[j] * i < Maxn++j)
       {
           int x = prim[j] * i;
           first[x] = 1;
           if (i % prim[j] == 0break;
           miu[x] = -miu[i];
       }
       miu[i] += miu[i - 1];
   }
}
int n;
ll solve(ll n)
{
   ll ret = 0;
   for (ll i = 1i <= 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之后会加深对莫比乌斯函数的理解。题解可以在网上搜,这里不再赘述

C++语言
#include
#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 = 2i < Maxn++i)
   {
       if (!first[i]) miu[i] = 1, prim[top++] = i;
       for (int j = 0j < top++j)
       {
           int x = prim[j] * i;
           if (x >= Maxnbreak;
           first[x] = 1;
           two[x] = two[i];
           if (i % prim[j] == 0two[x]++;
           if (miu[i] == 0miu[x] = 0;
           else miu[x] = miu[i] > 0 ? -(miu[i] + 1: -(miu[i] - 1);
           if (two[x] == 1miu[x] = (miu[x] > 0 ? 1 : -1);
           else if (two[x] > 1miu[x] = 0;
           if (i % prim[j] == 0break;
       }
       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 = 1i <= a;)
   {
       vec[vec.size++] = i;
       i = min(a / (a / i), b / (b / i)) + 1;
   }
   vec[vec.size++] = a + 1;
   for (int i = 1i < vec.size++i)
   {
       ll t = vec[i] - 1;
       ll s = vec[i - 1];
       ret += (miu[t] - miu[- 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;
}

M?bius反演(转),布布扣,bubuko.com

M?bius反演(转)

标签:des   style   blog   code   c   tar   

原文地址:http://www.cnblogs.com/qiucz/p/3724636.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!