注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

OI之路,漫漫人生

只为梦想,没有理由

 
 
 

日志

 
 

莫比乌斯反演小结  

2014-05-09 20:54:51|  分类: 算法杂类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

莫比乌斯反演小结:

有两个函数,F[X]和f[X],定义莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生,那么我们现在要根据F[N]推出小f[N]。首先对于F[P]=f[1]+f[P],将f[1]移动到右边变成F[P]-f[1]=f[P];这就是为什么叫它反演(*^__^*) 。。比如再对于f[6]=F[6]-f[1]-f[2]-f[3]=F[6]-F[2]-F[3]+F[1],因为F[2]=f[1]+f[2],F[3]=f[1]+f[3],f[1]=F[1].然后我们隐约的发现其实也是和F[N]的算法一样的,只是我们不太确定符号。那么我们可以先得到:莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生。这就是莫比乌斯反演公式!!;至于符号交给莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生函数了。

然后我们开始证明了:首先根据F(n)=f(n)得到u(1)=1

f(p)=u(1)*F(p)+u(p)*F(1)=F(p)-F(1),u(1)=1容易得到u(p)=-1;

f(p1*p2)=F(P1*P2)-f(1)-f(p1)-f(p2)=u(1)*F(p1*p2)+u(p1)*F(p2)+u(p2)*F(p1)+u(p1*p2)*F(1);

明显都可以抵消后u(p1*p2)=1;(p为质数)

以此类推下去。

对于双平方因子:容易推导u(p1^2)=0;又因为莫比乌斯函数是积性函数,所以只要都平方因子的数都是0;

总结一下就是这样:

(1)若莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生,那么莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生

(2)若莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生均为互异素数,那么莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生

    

(3)其它情况下莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生

而莫比乌斯函数还有几个有趣的性质:

(1)对任意正整数有莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生

其实可以用杨辉三角组合数枚举质因数种类证明当N>1时都为0,因为答案=-NC0+NC1-NC2+NC3-NC4+。。NCN=0;

(2)对任意正整数有莫比乌斯反演小结 - 山哥 - OI之路,漫漫人生

此证明略过。

然后给出以下一个莫比乌斯题的通用公式:记住下面这个基础公式,许多墨乌题转换为这个东西可以秒杀了。

如果题目让你求类似这样的:

I:[1..n] j:[1..m]

∑        ∑   e(gcd(I,j)=1)         e(x){  1    (x=1)}{0    x!=1}

可以发现就是 i=1 to n    j=1 to m   +u(d)  (满足d|gcd(i,j))

这样我们其实可以枚举d,也就是把d提出;

就只变成了 i=1 to n   +u(i)*(n/i)*(m/i);设n为较小的那个数。

证明: 因为在所有1-N和1-M的范围中,GCD(A,B)%d==0 这样的数一定有N/I*M/I个;所以得证。

那么效率就转换为oN),而枚举d的效率还可以用分块降为sqrt(N)

一开始要线性处理好u(d),关于线性筛莫比乌斯函数我有一份自己的模板:

inline void mobius()

{int i,j;

for (i=2;i<=N;i++){

if (!p[i])p[i]=prime[++tot]=i;

for (j=1;j<=tot&&prime[j]*i<=N;j++)

{p[i*prime[j]]=prime[j];

if (!(i%prime[j]))break;

}

}

mo[1]=1;

for (i=2;i<=N;i++)

if (p[i]==i)mo[i]=-1;

else if (p[i]==p[i/p[i]])mo[i]=0;

else mo[i]=-mo[i/p[i]];

}

现在基础的都讲完了,开始切题了:

1.http://www.lydsy.com/JudgeOnline/problem.php?id=2818

此题用欧拉函数可以搞定,用莫比乌斯也可以:枚举一个质数p,然后就转换为在[1..N/P][1..M/P]中求gcd(I,J)=1,即solve(n/p,m/p),转换为上面讲的公式了,所以套下模板即可了。代码略过;

2. 题意:给定两个数,其中,求为质数的有多少对?其中的范围是

此题题解转载:

本题与上题不同的是n和m不一定相同。于是新定义了一个g(x)函数,推导请看http://blog.csdn.net/acdreamers/article/details/8542292代码略过;

3. http://www.lydsy.com/JudgeOnline/problem.php?id=2440

此题就是考察了莫比乌斯的性质,u(d)就已经包含了容斥的思想。本来我们对于一个n,算平方数一定就I=1sqrt(N)然后 N/(I*I)就是有平方数因子的个数,然而我们会多减,比如3^2*2^2多算了。所以在减的时候乘上u(d)就行了,加上二分枚举即可。

代码:

#include<cstdio>

#include<cmath>

#define nn 50010

using namespace std;

int p[nn],prime[nn],mo[nn];

int l,r,tot,t,n,mid;

inline void read(int &x){

    char ch;int bo=0;x=0;

    for (ch=getchar();ch<'0'||ch>'9';ch=getchar())if (ch=='-')bo=1;

    for (;ch>='0'&&ch<='9';x=x*10+ch-'0',ch=getchar());

    if (bo)x=-x;

}

inline void mobius()

 {int i,j;

  for (i=2;i<=50000;i++){

              if (!p[i])p[i]=prime[++tot]=i;

  for (j=1;j<=tot&&prime[j]*i<=50000;j++)

                          {p[i*prime[j]]=prime[j];

                        if (!(i%prime[j]))break;

                          }

                        }

  mo[1]=1;

  for (i=2;i<=50000;i++)

             if (p[i]==i)mo[i]=-1;

       else if (p[i]==p[i/p[i]])mo[i]=0;

        else mo[i]=-mo[i/p[i]];

 }

inline bool calc(int w)

 {int j,s=w;

 for  (j=2;j<=sqrt(s);j++)

          w+=mo[j]*s/(j*j);

 if (w>=n)return true;

 else return false;

 }

int main()

{

 read(t);tot=0;

 mobius();

 while (t--)

      {

      read(n);

      l=n-1;r=n<<1;

      while (l+1<r)

           {

           mid=(1ll+l+r-1)>>1;

           if (calc(mid)) r=mid;

                    else l=mid;

           }

      printf("%d\n",r);

      }

return 0;

}

4. http://www.lydsy.com/JudgeOnline/problem.php?id=2005

此题用容斥也可以做,不过用莫乌也是简单题。

我们要求的是I=1 TO N    J=1 TO M (gcd(I,J)=k的对数*k*2-1)

然后我们提出并枚举k,就变成了

1-    K,然后用莫乌求[1..N/K],[1..M/K]gcd(I,j)的对数,最后减去那个1n*m,发现又转换成了重要公式,于是又秒掉了。。。

分块的时候由于是一块一块处理,于是转换成前缀和的形式,last-i即可。

 #include<cstdio>

#define maxn 100100

using namespace std;

typedef long long ll;

ll c,n,m,i,j,ans,tot,sum[maxn],p[maxn],prime[maxn],mo[maxn];

ll min(ll a,ll b){return (a>b)?b:a;}

void mobius()

 {ll i;

  for (i=2;i<=n;i++)

             {

       if (!p[i]){p[i]=prime[++tot]=i;}

       for (j=1;(j<=tot)&&(prime[j]*i<=n);j++)

                      {

                      p[prime[j]*i]=prime[j];

                      if (i%prime[j]==0)break;

                      }

             }

 mo[1]=sum[1]=1;

 for (i=2;i<=n;i++)

                {

       if (p[i]==i)mo[i]=-1;

       else if (p[i]==p[i/p[i]])mo[i]=0;

       else mo[i]=-mo[i/p[i]];

                }

 for (i=2;i<=n;i++)sum[i]=sum[i-1]+mo[i];

 }

ll num(ll w)

 {ll nn=n/w;ll mm=m/w;

  ll tem=0;

  for (ll i=1,last=0;i<=nn;i=last+1)

      {last=min(nn/(nn/i),mm/(mm/i));

       tem+=(sum[last]-sum[i-1])*(nn/i)*(mm/i);

      }

  return tem;

 }

int main()

{

 scanf("%lld%lld",&n,&m);

 if (n>m){c=n;n=m;m=c;}

 mobius();

 for (i=1;i<=n;i++)ans=ans+num(i)*2*i;

 printf("%lld\n",ans-n*m);

 return 0;

}

5. http://www.lydsy.com/JudgeOnline/problem.php?id=2301

这道题k已经给你了,于是不需要枚举,直接求[1..N/K][1..M/K]gcd(I,j)的个数,(*^__^*) 嘻嘻……又转换为那个公式了。由于此题范围有限制,直接用容斥搞一下就行了。

#include<cstdio>

#include<algorithm>

#define nn 60000

using namespace std;

typedef long long ll;

ll ans,a,b,c,d,k,i,tot,t,j,mu[nn+2],prime[nn+2],sum[nn+2],p[nn+2];

ll work(ll n,ll m)

 {ll tem=0,i,last=0;

  if (n>m)swap(n,m);

  for (i=1,last;i<=n;i=last+1)

         {

         last=min(n/(n/i),m/(m/i));

         tem+=(sum[last]-sum[i-1])*(n/i)*(m/i);

         }

  return tem;

 }

int main()

{

for (i=2;i<=nn;i++)

             {

             if (!p[i])p[i]=prime[++tot]=i;

     for (j=1;j<=tot&&(prime[j]*i<=nn);j++)

                        {

                        p[prime[j]*i]=prime[j];

                        if (i%prime[j]==0)break;

                        }

             }

mu[1]=sum[1]=1;

for (i=2;i<=nn;i++){

   if (p[i]==i)mu[i]=-1;else

   if (p[i]==p[i/p[i]])mu[i]=0;

   else mu[i]=-mu[i/p[i]];

   sum[i]=sum[i-1]+mu[i];

                  }   

scanf("%lld",&t);

while (t--)

     {scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k);

ans=work(b/k,d/k)-work((a-1)/k,d/k)-work(b/k,(c-1)/k)+work((a-1)/k,(c-1)/k);

printf("%lld\n",ans);

     }  

return 0;

}

  评论这张
 
阅读(13030)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018