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)即为莫比乌斯函数,从而可以迅速得解。

C++语言: 高亮代码由发芽网提供

#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 = 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之后会加深对莫比乌斯函数的理解。题解可以在网上搜,这里不再赘述

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 = 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;
}

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

时间: 2024-10-11 17:08:16

M?bius反演(转)的相关文章

bzoj 2820 / SPOJ PGCD 莫比乌斯反演

那啥bzoj2818也是一样的,突然想起来好像拿来当周赛的练习题过,用欧拉函数写掉的. 求$(i,j)=prime$对数 \begin{eqnarray*}\sum_{i=1}^{n}\sum_{j=1}^{m}[(i,j)=p]&=&\sum_{p=2}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{p}\rfloor}[i⊥j]\newline&=&\sum_{p=

【计算几何】【圆反演】hdu6097 Mindis

给你一个中心在原点的圆,再给你俩在圆内且到原点距离相等的点P,Q,让你在圆上求一点D,最小化DP+DQ. http://blog.csdn.net/qq_34845082/article/details/77099332 附:过反演中心的圆反演后变成一条和该圆正交的直线. 不过反演中心的圆反演后是一个与原圆关于反演中心位似的圆. 不过反演中心的直线反演后变成一个过反演中心且与其正交的圆. #include<cstdio> #include<cmath> using namespac

hdu1695(莫比乌斯反演)

题目链接: http://acm.hdu.edu.cn/showproblem.php?pid=1695 题意: 对于 a, b, c, d, k . 有 x 属于 [a, b],  y 属于 [c, d], 求 gcd(x, y) = k 的 x, y 的对数 . 其中 a = b = 1 . 注意: (x, y), (y, x) 算一种情况 . 思路: 莫比乌斯反演 可以参考一下: http://blog.csdn.net/lixuepeng_001/article/details/5057

算法学习——莫比乌斯反演(1)

.. 省选GG了,我果然还是太菜了.. 突然想讲莫比乌斯反演了 那就讲吧! 首先我们看一个等式-- (d|n表示d是n的约束) 然后呢,转换一下 于是,我们就发现! 没错!F的系数是有规律的! 规律is here! 公式: 这个有什么卵用呢? 假如说有一道题 F(n)可以很simple的求出来而求f(n)就比较difficult了,该怎么办呢? 然后就可以用上面的式子了 是莫比乌斯函数,十分有趣 定义如下: 若d=1,则=1 若d=p1*p2*p3...*pk,且pi为互异素数,则=(-1)^k

bzoj2301 [HAOI2011]Problem b【莫比乌斯反演 分块】

传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2301 很好的一道题.首先把每个询问转化为4个子询问,最后的结果就是这四个子询问的记过加加减减,类似二维前缀和.那么问题转化为在1 <= x <= lmtx, 1 <= y <= lmty时gcd(x, y) == k的对数,这个问题在转化一下,转化成1 <= x <= lmtx / k,1 <= y <= lmty / k时x与y互质的对数.莫比乌斯反

BZOJ2301: [HAOI2011]Problem b 莫比乌斯反演

分析:对于给出的n个询问,每次求有多少个数对(x,y),满足a≤x≤b,c≤y≤d,且gcd(x,y) = k,gcd(x,y)函数为x和y的最大公约数. 然后对于求这样单个的gcd(x,y)=k的,我们通常采用莫比乌斯反演 但是,时间复杂度是O(n*(n/k))的,当复杂度很坏的时候,当k=1时,退化到O(n^2),超时 然后进行分块优化,时间复杂度是O(n*sqrt(n)) #include<cstdio> #include<cstring> #include<queue

BZOJ2005: [Noi2010]能量采集 莫比乌斯反演的另一种方法——nlogn筛

分析:http://www.cnblogs.com/huhuuu/archive/2011/11/25/2263803.html 注:从这个题收获了两点 1,第一象限(x,y)到(0,0)的线段上整点的个数是gcd(x,y) 2,新学了一发求gcd(x,y)=k有多少对的姿势,已知0<x<=n,0<y<=m 令x=min(n,m),令f[i]代表gcd(x,y)=i的对数, 那么通过O(xlogx)的复杂度就可以得到f[1]到f[n](反着循环) 普通的容斥(即莫比乌斯反演)其实也

容斥原理与莫比乌斯反演的关系

//容斥原理,c[i]表示i当前要算的次数,复杂度和第二层循环相关 O(nlogn~n^2) LL in_exclusion(int n,int *c) { for(int i=0;i<=n;i++) c[i]=1; //不一定是这样初始化,要算到的才初始化为1 LL ans=0; for(int i=0;i<=n;i++) if(i要算) { ans+=(统计数)*c[i]; for(int j=i+1;j<=n;j++) if(i会算到j) c[j]-=c[i];//j要算的次数减去

二项式反演

问:给你k种颜色,你必须用上所有颜色去涂满n个相邻的格子,并且要求相邻格子的颜色不同,求方案数. 我们设必须用 i 种颜色两两不相邻的涂格子的方案数为 b(i) ; 很明显: ,我们令 a(k)=k·(k-1)n-1 , 然后有. 如果你知道二项式反演的话,那么这个问题就已经解决了,因为. 是不是觉得二项式反演很厉害,下面我将给出它的证明. 二项式反演公式: 证明: 然后让我们对进行分析: 我们预热一下: 有A,B,C,D,E,F,G 7个人,我们要先从中选出4个候选人,再从中选出3个作为mas