莫比乌斯算法浅析_莫比乌斯算法-CSDN博客

转载自：http://blog.csdn.net/herodeathes/article/details/77932208
【Problem Description】

Given 5 integers: a, b, c, d, k, you’re to find x in a…b, y in c…d that GCD(x, y) = k. GCD(x, y) means the greatest common divisor of x and y. Since the number of choices may be very large, you’re only required to output the total number of different number pairs.

Please notice that, (x=5, y=7) and (x=7, y=5) are considered to be the same.

Yoiu can assume that a = c = 1 in all test cases.

【Input】
The input consists of several test cases. The first line of the input is the number of the cases. There are no more than 3,000 cases.
Each case contains five integers: a, b, c, d, k, 0 < a <= b <= 100,000, 0 < c <= d <= 100,000, 0 <= k <= 100,000, as described above.

【Output】
For each test case, print the number of choices. Use the format in the example.

【Sample Input】
2
1 3 1 5 1
1 11014 1 14409 9

【Sample Output】
Case 1: 9
Case 2: 736427

【Hint】
For the first sample input, all the 9 pairs of numbers are (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 5), (3, 4), (3, 5).

        颓废了很久，浪子回头了......收拾一下状态， 还得加油上路。

好吧回到正题，最近看到莫比乌斯反演（好吧是老师布置的），感觉很有意思（好吧还是老师说的），于是开始认（zao）认（meng）真（xi）真（you）的学习。感觉展示我自己英语水平的时刻又到了。真相只有一个——

对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。对于100%的数据，a=c=1。

学完这个专题我觉得我在学业上又（shi）进（fen）一（tong）步（ku）。
我现在都不想多说话了，来介绍下莫比乌斯反演吧。

懵逼钨丝反演是个……啊不，莫比乌斯反演简单来说是一个公式，用这个公式可以通过一个已知函数，推导出另一个与这个函数有关的函数的值。好像有点拗口？？没事，程序员是不需要在意这些细节的。在这里只简略讲一下。它起到的简化数论运算的效果是其他算法很难比拟的。

公式如下：G(n)=sigma(F(d)) (其中d|n) （集合包含）
其实反演还有另一个集合公式，但是那啥…是废的，ACM没出现过
所以先学这个。
通过这个公式，我们可以得出F与G的关系：
G(1) = F(1)
G(2) = F(1)+F(2)
G(3) = F(1)+F(3)
G(6) = F(1)+F(2)+F(3)+F(6)
然后就是容斥原理，因为用G来表达F的值的时候，会用到他（可以在草稿纸上列一下，例如F(2)=G(2)-G(1)）。所以每个式子前面的G都有一个系数（-G和G，哦对了，还有0*G，也就是不出现）,然后我们通过欧拉函数得到求解F(n)时，G(d)的系数。
那么剩下的问题就是每个G的系数！！！
我们以求解 F(6)为例子来说明，并定义一个系数函数 H(d,n).
其中 H(d,n)表示求解F(n)时，G(d)的系数（其中d|n）
所以可以得到这个式子
F(6) = H(6,6)* G(6)+H(2,6) * G(2)+H(3,6)G(3)+H(1,6) G(1)
我们用 a,b,c,d分别替代四个H(6,6),H(2,6),H(3,6),H(1,6)，并且把对应的G用F表示出来，得到
F(6)=a * (F(6)+F(3)+F(2)+F(1))+b* (F(2)+F(1))+c* (F(3)+F(1))+d*F(1),再变形一下，又有
F(6) * (a-1)+F(3) * (a+c)+F(2)* (a+b)+F(1)*(a+b+c+d)=0，把F(6),F(3),F(2),F(1)当作不同的元，则得到了下面的方程组！！！
a-1==0
a+c==0
a+b==0
a+b+c+d==0

这样总能找到解，别问我为什么，那啥方程式懂吗？？

那么再给出一个结论，大家可以自行证明，计算过程不算复杂。
H(a,b) 都可以写成 H(1,b/a) ，于是我们把H的第一个元素略去，简写为 H(x)
说到这里，就可以把H和U联系起来了，其实 U(x) = H(x) = H(1,x)
再来，我们就可以给U(x)赋予一个更具体的意义， U(x)表示在计算 F(x)时，G(1)的系数！！(因为U(x)==H(1,x))

接下来，我们来尝试一下，如何用上面那个U(x)的新意义，来计算U(x)的值！！
首先需要明确2点！
一是G(x)中,一定包含一个F(1)，因为 1|x
二是，F(1)==G(1)

(0).如果 x==1
因为 F(1)==G(1) 所以 U[1]=1;

(1).假设 x 是一个质数
F(x) = U(1)G(x)+U(x) G(1)
带入U(1) == 1, 因为G(x)中含有一个F(1),而左边不含F(1)，所以我们需要利用G(1)来消去F(1)
所以得到 U(x)=-1

(2).假设 x 可以写成2个不同质数的乘积 x=p*q
那么 F(x)=U(1)*G(pq)+U(q)*G(p)+U(p)*G(q)+U(x)*G(1)
这里 U(1)，U(p)，U(q) 就是前面2种情况
带入系数，因为左边没有 F(1)，所以为了抵消右边的F(1)，我们需要令 U(x)=1;

(3).假设 x 可以写成3个不同质数的乘积 x=p*p1*p2 我们令 z = p1*p2
F(x) = U(1) * G(pz)+ U(z) * G(p)+U(p) G(z)+U(x) G(1);
其中 U(1),U(p),U(z) 分别为前面几种情况，带入过后，为抵消F(1) 得到 U(x)=-1
由此可以相同的方式向下递推，得到第一条结论
如果 x = p1* p2…* pr , 其中pi是互异的质数，那么
U[x] = (-1)^r ———————– 1！！

(4).假设 x 可以写成一个质数的平方 x=p^2
F(x) = U(1)* G(x)+U(p)* G(p)+U(x)*G(1)
带入系数得到 U(x)=0;

(5).假设 x 可以写成一个质数的三次方 x=p^3
F(x) = U(1)*G(x)+U(p)*G(p^2)+U(p^2)*G(p)+U(x)*G(1)
带入系数后 U(x)=0;
由此可用相同方式向下递推，得到第二条结论
如果 x = p^e (e>1) U[x] = 0; ————————– 2！！

(6).假设 x 可写成 x = p^e*q 其中p,q为不同质数，e>1
F(x) = U(1)*G(x)+U(q)*G(p^e)+U(p^e)*G(q)+U(x)*G(1)
带入系数后 U(x) = 0;

由此可继续向下递推，得到第二条结论的加强版！！
如果 x = p^e*z 其中p为质数, z为任意数,e>1 那么 U[x] = 0

差不多了。

以为结束了？要介绍代码了？嘿嘿，来，继续折磨你
这道题怎么做
先说坑点吧。这道题目有比较恶心的一点，数据有k==0的，这时显然答案是0，没有2个数的gcd为0。首先，gcd是没啥用的。因为约掉gcd后两个数互质。于是我们可以让x/=k y/=k并且假设 x<=y
然后题目变成了 2个数分别在区间[1..x]和[1..y]中的互质数有多少对。
这道题目还要求1-3 和 3-1 这种情况算成一种，因此只需要限制x

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm> using namespace std; const int maxn=1000000; bool check[maxn+10]; //表示这个数是不是素数 int prime[maxn+10]; //储存素数的数组 int U[maxn+10];//表示系数  void get_Mobius() {
    memset(check,0,sizeof(check));
    U[1]=1;
    int tot=0;
    for(int i=2;i<=maxn;i++)
    {
        if(check[i]==0)
        {
            prime[tot++]=i;
            U[i]=-1;//本身是个质数，1个质数乘积 
        }
        for(int j=0;j<tot;j++)
        {
            if(i*prime[j]>maxn)break;
            check[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                U[i*prime[j]]=0;//除开质数的两种情况 
                break;
            }
            else
            {
                U[i*prime[j]]=-U[i];//多了一个数的乘积，与前一个相反 
            }
        }
    } } /* 只要让F（t）=满足gcd(x,y)%t==0的数对个数  f(t)=满足gcd(x,y)=t的数对个数，则F（t）和f(t)就存在莫比乌斯反演的关系了。显然F（t）=(b/t)*(d/t)  因为如果gcd(x,y)=1,则gcd(x?k,y?k)=k，所以我们把b和d同时除以k， 得到的f(1)再去重就是答案。令lim=min(b/k,d/k),就根据我给出的莫比乌斯反演第二个公式， 
*/ int main() {;
    int T;
    int a,b,c,d,k;
    get_Mobius();
    scanf("%d",&T);
    for(int ii=1;ii<=T;ii++)
    {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k==0)
        {
            printf("0\n");
            continue;
        }
        b/=k;a/=k;
        d/=k;c/=k;
        if(b>d)swap(b,d);
        long long ans1=0;
        for(int i=1;i<=b;i++)
            ans1+=(long long)U[i]*(b/i)*(d/i);//带入 G(n) 的公式 有  F(1) = sigma(U[i]*(a/i)*(b/i)) (1<=i<=N)
        long long ans2=0;
        for(int i=1;i<=b;i++)
            ans2+=(long long)U[i]*(b/i)*(b/i);//相同的数啊如此重复的状况 
        ans1-=ans2/2;
        printf("%lld\n",ans1);
    }
    return 0; }