莫比乌斯反演

来源：互联网发布：辽宁联通网络测速编辑：程序博客网时间：2024/04/29 23:42

莫比乌斯反演

分类：数论2013-01-25 15:26 111人阅读评论(0) 收藏举报

莫比乌斯函数值：

[cpp] view plaincopy
int mobi(int n)  
{  
    int m=1;  
    for(int i=2;i*i<=n;i++)  
    {  
        if(n%i==0)  
        {  
            m*=-1;  
            int k=0;  
            do  
            {  
                k++;  
                if(k>1)  
                {  
                    m=0;  
                    break;  
                }  
                n/=i;  
            }while(n%i==0);  
        }  
    }  
    if(n>1)  m*=-1;  
    return m;  
}  

题目：GCD

对于给定的N(<=10^7)，求使得gcd(x,y)为质数的有序数对(x,y)(1<=x,y<=N)的个数。

显然ans=sigma(f[N/p])，1<p<=N且p为质数，f[k]表示1<=x,y<=k中互质的有序数对的个数。

因为f[k]=1+ 2*sigma(phi[i]) (1<i<=k)，phi是欧拉函数，线性筛搞一搞即可。

[cpp] view plaincopy
#include <cstdio>  
#include <algorithm>  
#include <cstdlib>  
#include <cstring>  
#include <cmath>  
const int fin=1,maxn=10000001;  
int N,tot;  
int pr[maxn];  
bool flag[maxn];  
long long f[maxn];  
int main()  
{  
    scanf("%d",&N);  
    int i,j;  
    long long ans=0;  
    for(i=2;i<=N;++i)  
    {  
        if(!flag[i])  
        {  
            pr[tot++]=i;  
            f[i]=i-1;  
        }  
        for(j=0;j<tot&&i*pr[j]<=N;++j)  
        {  
            flag[i*pr[j]]=true;  
            if(i%pr[j]==0)  
            {  
                f[i*pr[j]]=f[i]*pr[j];  
                break;  
            }  
            else  
                f[i*pr[j]]=f[i]*(pr[j]-1);  
        }  
    }  
    for(i=2;i<=N;++i)  
        f[i]+=f[i-1];  
    for(i=1;i<=N;++i)  
        f[i]=1+(f[i]<<1);  
    for(i=0;i<tot;++i)  
        ans+=f[N/pr[i]];  
    printf("%I64d",ans);  
    return 0;  
}  

题目：YY的GCD

题意：给定N,M，求满足1<=i<=N,1<=j<=M,gcd(i,j)为质数的有序数对(i,j)的对数，N,M<=10^7 测试数据组数<=10^4

莫比乌斯反演好神奇...

一直只知道百度百科上说的那一种形式：若F(n)=Σ(d|n)f(d)则f(n)=Σ(d|n)miu(d)F(n/d)

查了题解才知道有另一种形式...也就是这道题会用到的形式...

若F(x)=Σ(x|d)f(d)则f(x)=Σ(x|d)miu(d/x)F(d) 其中d在某个限定范围内，接下来的问题就是：

对于给定的N,M，求1<=i<=N,1<=j<=M,gcd(i,j)=1的有序数对(i,j)的对数。再枚举质数p后，把N'=N/p,M'=M/p时这个问题的解累加起来就可以了。

对于给定的N,M，设f(x)为1<=i<=N,1<=j<=M,gcd(i,j)=x的(i,j)的对数，F(x)为1<=i<=N,1<=j<=M,x|gcd(i,j)的(i,j)的对数。

显然有F(x)=(N/x)*(M/x) (N/x和M/x都是下取整) F(x)=Σ(x|d)f(d)

反演得到
f(x)=Σ(x|d)miu(d/x)F(d)=Σ(x|d)miu(d/x)*(N/d)*(M/d)

这个问题的解是f(1)，把x=1代入得到：f(1)=Σ(d)miu(d)*(N/d)*(M/d)

做完了?当然没有。不TLE见鬼了。考察答案的表达式：Ans=Σ(质数p) [ Σ(d)miu(d)*(N/pd)*(M/pd) ]

换个角度，先枚举T=pd，再去枚举p，则d=T/p。得到Ans=∑(T) [ ∑(质数p|T)miu(T/p)*(N/T)*(M/T) ]

这么做的好处立刻就显现出来了，N/T和M/T与p无关!

Ans=∑(T) [ (N/T)*(M/T) *∑(质数p|T)miu(T/p) ]

再令g(x)=sigma(质数p|x)miu(x/p) 改写上式为

Ans=∑(T)(N/T)*(M/T)*g(T)

恩，非常和谐。

接下来全副精力对付g(x)。仍然考虑线性筛搞定他。

假设枚举到i，质数枚举到p(程序里的prime[j])，要更新A=i*p的信息。

1. p|i
这时A的素数分解式中，p这一项的次数>=2。

考虑g(A)的求和式：

如果枚举的质数p'不等于p，A/p'就也会有p这一项，次数>=2，这时候miu(A/p')=0

如果枚举的质数p'=p，A/p=i，这一项就是miu(i)

因此g(A)=miu(i)
2. p!|i (即i%p!=0）

这时候A比i多一个质因子p，对miu(i)分情况讨论。

2.1 miu(i)==0 (即i有大于1次的项)

这时A除去任何一个p'都会留下i的那个大于1次的项，除非是下面这一种非常特殊的情况：

2.1.1 i的素数分解式中，大于1次的项只有一个，且这一项为2次。记这一项为p0。

这时除去任何一个p'!=p0都会留下这一项，但是除去p0则会得到A/p0——这个数所有的项都是1次的。因此g(A)=miu(A/p0)

2.1.2 i的素数分解式大于1次的项不止一个或者大于1次的项唯一，但次数高于2次。易见g(A)=0

2.2 miu(i)!=0 (即i全是1次) 这个时候A的项也全是1次。设r(x)为x的质因子个数。

则可以得到g(A)=r(A)*(-1)^(r(A)-1)。因为除以任何一个p'，miu(A/p')都是一样的。

同理g(i)=r(i)*(-1)^(r(i)-1)，且有r(A)=r(i)+1。利用r(A)=r(i)+1可以方便地得到：g(A)和g(i)异号，且绝对值比g(i)多1。

亦即g(A)=(g(i)>0)?-1:1 -g(i)

g(A)的维护讨论完了。

完了?没完，看情况2.1.1，我们有这么个遗留问题：

如果x的大于1次的项唯一，且这一项为2次，则令f(x)为这个项，否则f(x)=1。

事实上f(x)=1包含3种情况：

1. 大于1的项不唯一

2. 大于1次的项唯一但大于2次。

3. 全为1次

1和2利用现有的结果无法区分，但事实上不需要区分。3则可以用miu(x)判出来。

好，我们来对付f(x)，仍然是线性筛，变量意义同g(x)的讨论。

1. p|i

A由i把最小因子p的次数加1得到，显然这一项的次数>=2。

1.1 f(i)!=1

1.1.1 如果f(i)=p，那么A中p的次数就是3次了，f(A)=1。

1.1.2 如果f(i)!=p，那么A中大于1次的项就不唯一了，仍有f(A)=1

因此f(i)!=1必然有f(A)=1

1.2 i全为1次即f(i)=1且miu(i)!=0 这时显然f(A)=p

1.3 i不全为1次即f(i)=1且miu(i)=0 这时显然f(A)=1

2. p!|i

A比i多一个1次的质因数p，那么应有f(A)=f(i)

f(A)的递推也讨论完了。

完了?虽然剩下的工作很简单但是也是必不可少的..

回去看求和的式子：Ans=∑(T)(N/T)*(M/T)*g(T)

直接做是O(min(M,N))的，别忘了有1W组数据啊。

但是我们有个结论：对于给定的N，(N/T)的取值只有sqrt(N)个

那么给定的N,M，(N/T)*(M/T)就只有sqrt(N)+sqrt(M)个了，而且相同的取值当然是连成一段的。

由此，分段来计算。令gs(x)为g(x)的前缀和。

对于枚举到的i，我们希望找到最大的使得(N/i)*(M/i)=(N/j)*(M/j)的j值。

这可能有点困难，我们退而求其次，求出最大的使得N/i=N/j且M/i=M/j的j值，效果不会变差太多。

这个就很简单了，可以得到j=min(N/(N/i),M/(M/i))。于是把i~j这一段整体计算即可。

[cpp] view plaincopy
#include <cstdio>  
#include<algorithm>  
#include<cstdlib>  
#include<cstring>  
const int maxp=10000001;  
int pr[maxp],miu[maxp],g[maxp],f[maxp];  
long long gs[maxp];  
bool flag[maxp];  
int gcd(int x,int y)  
{  
    int t;  
    if(x<y)  
    {  
        t=x;  
        x=y;  
        y=t;  
    }  
    while(y)  
    {  
        t=x%y;  
        x=y;  
        y=t;  
    }  
    return x;  
}  
  
void init()  
{  
    int i,j,A;  
    miu[1]=1;  
    for(i=2;i<maxp;++i)  
    {  
        if(!flag[i])  
        {  
            pr[++pr[0]]=i;  
            miu[i]=-1;  
            f[i]=g[i]=1;  
        }  
        for(j=1;j<=pr[0]&&i*pr[j]<maxp;++j)  
        {  
            flag[A=i*pr[j]]=true;  
            if(i%pr[j]!=0)  
            {  
                miu[A]=-miu[i];  
                f[A]=f[i];  
                if(miu[i]==0)  
                     g[A]=(f[i]!=1)?miu[A/f[i]]:0;  
                else  
                     g[A]=(g[i]>0)?(-g[i]-1):(-g[i]+1);  
            }  
            else  
            {  
                miu[A]=0;  
                if(f[i]==1)  
                     f[A]=(miu[i]==0)?1:pr[j];  
                else  
                     f[A]=1;  
                g[A]=miu[i];  
                break;  
            }  
        }  
    }  
    for(i=2;i<maxp;++i)  
         gs[i]=gs[i-1]+g[i];  
}  
  
int main()  
{  
    init();  
    int T,N,M,i,j,t,t1,t2;  
    long long ans;  
    scanf("%d",&T);  
    while(T--)  
    {  
        scanf("%d%d",&N,&M);  
        if(N<M){t=N;N=M;M=t;}  
        ans=0;  
        for(i=2;i<=M;i=t+1)  
        {  
            t1=N/(N/i);  
            t2=M/(M/i);  
            t=(t1<t2)?t1:t2;  
            ans+=(gs[t]-gs[i-1])*(N/i)*(M/i);  
        }  
        printf("%lld\n",ans);  
    }  
    return 0;  
}  