莫比乌斯反演例题集 ^_^

最新推荐文章于 2023-02-08 19:39:32 发布

weixin_34268753

最新推荐文章于 2023-02-08 19:39:32 发布

阅读量208

点赞数

原文链接：http://www.cnblogs.com/PaperCloud/p/10010260.html

版权

【luogu 2257】YY的GCD

Problem Here

预备知识

除法分块、莫比乌斯反演

最终公式：

\[ans= \sum_{T=1}^{n} (\tfrac{n}{T}) (\tfrac{m}{T}) \sum_{p|t,isprime[p]=1} \mu ( \tfrac{T}{p} )\]

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 10000000
int T,n,m;
int prime[MN+5],cnt,mu[MN+5];
bool mark[MN+5];
ll sum[MN+5],ans;
inline void init(){
    mark[1]=1;mu[1]=1;register int i,j;
    for(i=2;i<=MN;i++){
        if(!mark[i]) {mu[i]=-1;prime[++cnt]=i;}
        for(j=1;j<=cnt&&prime[j]*i<=MN;j++){
            mark[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(i=1;i<=cnt;i++)
    for(j=1;j*prime[i]<=MN;j++) sum[j*prime[i]]+=mu[j];
    for(i=2;i<=MN;i++) sum[i]+=sum[i-1];
}
int main(){
    T=read();init();
    register int i,r;
    while(T--){
        n=read();m=read();int MIN=min(n,m);
        ans=0ll;
        for(i=1;i<=MIN;i=r+1){
            r=min(n/(n/i),m/(m/i));
            ans+=(1ll)*(n/i)*(m/i)*(sum[r]-sum[i-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

【HAOI2011】Problem b

Problem Here

预备知识

容斥

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 50005
int prime[MN+5],cnt,mu[MN+5];
bool mark[MN+5];
inline void init()
{
    mark[1]=1;mu[1]=1;register int i,j;
    for(i=2;i<=MN;i++)
    {
        if(!mark[i]) {mu[i]=-1;prime[++cnt]=i;}
        for(j=1;j<=cnt&&prime[j]*i<=MN;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(i=2;i<=MN;i++) mu[i]+=mu[i-1];
}
ll q(int n,int m,int k)
{
    ll ans=0ll;
    register int i,r;
    for(i=1;i*k<=n&&i*k<=m;i=r+1){
        r=min(n/(n/(i*k))/k,m/(m/(i*k))/k);
        ans+=1ll*(mu[r]-mu[i-1])*(n/(i*k))*(m/(i*k));
    }
    return ans;
}
int main()
{
    int n=read(),a,b,c,d,k;
    init();
    while(n--)
    {
        a=read(),b=read();
        c=read(),d=read();
        k=read();
        printf("%lld\n",q(b,d,k)-q(a-1,d,k)-q(b,c-1,k)+q(a-1,c-1,k));
    }
    return 0;
}

【NOI2010】能量采集

Problem Here

预备知识

求gcd的新姿势

\[\sum_{d|n}^{} \phi(d) =n\]

prove: 考虑小于等于n的数中，与n的gcd为\(\frac{n}{d}\)的数的个数是\(\phi(d)\),而所有的\(\phi(d)\)相加正好是n

所有满足 i|x,i|y的数，均满足i|gcd(x,y)

//一个正整数，等于它所有因数的欧拉函数之和。
#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 100000
int prime[MN+5],cnt;
ll phi[MN+5],ans;
bool mark[MN+5];
inline void init()
{
    mark[1]=1;phi[1]=1;register int i,j;
    for(i=2;i<=MN;i++)
    {
        if(!mark[i]) {phi[i]=i-1;prime[++cnt]=i;}
        for(j=1;j<=cnt&&prime[j]*i<=MN;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            else phi[i*prime[j]]=(prime[j]-1)*phi[i];
        }
    }
    for(i=2;i<=MN;i++) phi[i]+=phi[i-1];
}
int main(){
    register int n,m;
    n=read(),m=read();init();
    register int i,r;
    for(i=1;i<=n&&i<=m;i=r+1){
        r=min(n/(n/i),m/(m/i));
        ans+=1ll*(n/i)*(m/i)*(phi[r]-phi[i-1]);
    }
    (ans<<=1)-=1ll*n*m;
    printf("%lld\n",ans);
    return 0;
}

【luogu 1829】Crash的数字表格

Problem Here

预备知识

有一些莫比乌斯反演的的常见套路：

对于求和gcd的，把gcd提到最前面枚举
对于判断一个数是否为1，用\[\sum_{d|n}^{} \mu(d)=[n==1]\]来实现，之后在想办法把\(\mu(d)\)给提出来
把能够进行除法分块的部分提出来

最终公式:

\[ans=\sum_{T=1}^{n} Sum(\tfrac{n}{T}) Sum(\tfrac{m}{T}) T \sum_{d|T} d \mu(d)\]

而函数 \[g(x)= \sum_{d|x}^{} d \mu(d) \]是个积性函数，可以用欧拉筛来求。

//把能够除法分块的部分拖到最前面
//积性函数可以用欧拉筛来求 
#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 10000007
#define mod 20101009
int prime[MN],cnt;
ll f[MN],ans;
bool mark[MN];
inline void init()
{
    mark[1]=1;f[1]=1;register int i,j;
    for(i=2;i<=MN-7;i++)
    {
        if(!mark[i]==1){prime[++cnt]=i;f[i]=1-i+mod;}
        for(j=1;j<=cnt&&i*prime[j]<=MN-7;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j]==0){f[i*prime[j]]=f[i];break;}
            else f[i*prime[j]]=f[i]*f[prime[j]]%mod;
        }
    }
    for(i=1;i<=MN-7;i++) f[i]=f[i]*i%mod;
    for(i=2;i<=MN-7;i++) (f[i]+=f[i-1])%=mod; 
}
inline ll sum(int x){return (1ll*x*(x+1)/2)%mod;}
int main()
{
    register int i,r,n,m,M;
    n=read(),m=read();M=min(n,m);
    init();
    for(i=1;i<=M;i=r+1){
        r=min(n/(n/i),m/(m/i));
        (ans+=(sum(n/i)*sum(m/i)%mod*(mod+f[r]-f[i-1])%mod)%mod)%=mod;
    }
    printf("%lld\n",ans);
    return 0;
}

【luogu 3327】[SDOI2005] 约数个数和

Problem Here

预备知识

\[d(nm)=\sum_{i|n}^{} \sum_{j|m} [gcd(i,j)==1]\]

prove:

考虑分别理解左式和右式：

取\(nm\)的约数个数和相当于各个质因子的（次数\(+1\)）的累乘
试想该如何分配\(i\)，\(j\)中包含质因子\(p\)的次数？显然，只有（次数\(+1\)）种做法，累乘即为答案
\[\sum _{i=1}^{n} \frac{n}{i}= \sum_{j=1}^{n} d(j)\] ，从右往左理解要方便的多

具体可以参见约束研究
d(n)是一个积性函数。

最终公式:

令 \[g(n)=\sum _{i=1}^{n}\frac{n}{i}\]

\[ans=\sum_{d=1}^{min(n,m)} \mu(d) \ g(\tfrac{n}{d}) g (\tfrac{m}{d})\]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 50005
int mu[MN],cnt,prime[MN];
ll d[MN],p[MN];
bool mark[MN];
inline void init()
{
    mark[1]=true;
    p[1]=d[1]=mu[1]=1;
    register int i,j;
    for(i=2;i<=MN-5;++i)
    {
        if(!mark[i]){prime[++cnt]=i;mu[i]=-1;d[i]=2;p[i]=2;}
        for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
        {
            mark[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;d[i*prime[j]]=d[i]/p[i]*(p[i]+1);
                p[i*prime[j]]=p[i]+1;break; 
            }
            else
            {
                mu[i*prime[j]]=-mu[i];d[i*prime[j]]=d[i]*2;
                p[i*prime[j]]=2;
            }
        }
    }
    for(i=2;i<=MN-5;i++) mu[i]+=mu[i-1],d[i]+=d[i-1];
} 
int main()
{
    register int i,r,n,m,T;
    T=read(),init();
    while(T--)
    {
        register ll ans=0ll;
        n=read(),m=read();
        for(i=1;i<=n&&i<=m;i=r+1)
        {
            r=min(n/(n/i),m/(m/i));
            ans+=1ll*d[n/i]*d[m/i]*(mu[r]-mu[i-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

【SDOI2014】数表

Problem Here

预备知识

莫比乌斯反演+树状数组！

“同时整除i和j的所有自然数之和”其实可以理解为\(\sum_{d|gcd(i,j)}^{} d\),也就是\(\gcd(i,j)\)的约数和

而每个数的约数和是可以在O(n log n) 的时间内算出的
记F[i]为i的约数和
根据莫比乌斯反演的基本操作，我们要求的是\[\sum_{d=1}^{n} \left \lfloor \frac{n}{d} \right \rfloor\left \lfloor \frac{m}{d} \right \rfloor\sum_{x|d}^{} F(x)\mu (\frac{d}{x})\]
因为加上了a的限制，先对询问按a排序，采用树状数组维护前缀和

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 100005
#define MQ 20005
struct ques
{
    int n,m,a,id;
    bool operator < (const ques&o) const{return a<o.a;}
}q[MQ];
struct pac
{
    int a,id;
    bool operator < (const pac&o) const{return a<o.a;}
}F[MN];
//F[x]:x的约数和 
//mu[x]：莫比乌斯函数 
int mu[MN],prime[MN],cnt;
bool mark[MN];
inline void init()
{
    mark[1]=1;mu[1]=1;
    register int i,j;
    for(i=1;i<=MN-5;++i) for(j=1;j*i<=MN-5;++j) F[i*j].a+=i;
    for(i=1;i<=MN-5;++i) F[i].id=i; 
    std::sort(F+1,F+MN-4);
    for(i=2;i<=MN-5;++i)
    {
        if(!mark[i]){prime[++cnt]=i;mu[i]=-1;}
        for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
        {
            #define k i*prime[j] 
            mark[k]=1;
            if(i%prime[j]==0){mu[k]=0;break;}
            else mu[k]=-mu[i];
            #undef k
        }
    }
}
#define lowbit(x) (x&(-x))
int t[MN];
inline void C(int x,int val){for(;x<=MN-5;x+=lowbit(x)) t[x]+=val;}
inline int G(int x){int res=0;for(;x;x-=lowbit(x)) res+=t[x];return res;}
int ans[MQ];
int main()
{
    register int i,r,Q,now=1;
    Q=read();init();
    for(i=1;i<=Q;++i) q[i].n=read(),q[i].m=read(),q[i].a=read(),q[i].id=i;
    std::sort(q+1,q+Q+1);
    for(int t=1;t<=Q;++t)
    {
        #define N q[t].n
        #define M q[t].m
        #define A q[t].a
        for(;F[now].a<=A&&now<=MN-5;++now)
            for(i=1;i*F[now].id<=MN-5;++i)
                C(i*F[now].id,F[now].a*mu[i]);
        for(i=1;i<=N&&i<=M;i=r+1)
        {
            r=min(N/(N/i),M/(M/i));
            ans[q[t].id]+=(M/i)*(N/i)*(G(r)-G(i-1));
        }
        #undef N
        #undef M
        #undef A 
    }
    for(i=1;i<=Q;++i) printf("%d\n",ans[i]&0x7fffffff);
    return 0; 
}

【SDOI2017】数字表格

Problem Here

预备知识

其实没有预备知识der~，都是套路

最终公式:

\[ans=\prod_{T=1}^{n} \left ( \prod_{d|T}^{ } f[d]^{\mu (\frac{n}{T})}\right )^{\left [ \frac{n}{T} \right ]\left [ \frac{m}{T} \right ]}\]

还有就是要小心别T了......

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
#define MN 1000005
#define mod 1000000007
#define swap(x,y) (x^=y^=x^=y)
int mu[MN],prime[MN],cnt/*,fib[MN]*/,f[MN],invf[MN];
bool mark[MN];
inline int fpow(int x,int m)
{
    int res=1;
    while(m)
    {
        if(m&1) res=1ll*res*x%mod;
        x=1ll*x*x%mod;m>>=1;
    }
    return res;
}
inline void init()
{
    register int i,j,k;
    //fib[1]=fib[2]=1;
    //for(i=3;i<=MN-5;++i) fib[i]=fib[i-1]+fib[i-2],fib[i]>mod?fib[i]-=mod:0;
    mark[1]=1;mu[1]=1;
    for(i=2;i<=MN-5;++i)
    {
        if(!mark[i]){prime[++cnt]=i;mu[i]=-1;}
        for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
        {
            #define K i*prime[j] 
            mark[K]=1;
            if(i%prime[j]==0){mu[K]=0;break;}
            else mu[K]=-mu[i];
            #undef K
        }
    }
    for(i=1;i<=MN-5;++i) f[i]=invf[i]=1;
    int A=1,B=0;
    for(i=1;i<=MN-5;++i){
        B=(A+B)%mod;A=(B-A+mod)%mod;
        int G[3]={fpow(B,mod-2),1,B};
        for(j=i,k=1;j<=MN-5;j+=i,++k){
            f[j]=(ll)f[j]*G[mu[k]+1]%mod,
            invf[j]=(ll)invf[j]*G[1-mu[k]]%mod;
        }
    }
    f[0]=invf[0]=1;
    for(i=1;i<=MN-5;++i)
        f[i]=(ll)f[i-1]*f[i]%mod,
        invf[i]=(ll)invf[i-1]*invf[i]%mod; 
}
int main()
{
    register int t=read(),i,r,n,m;init();
    while(t--)
    {
        int ans=1;
        n=read(),m=read();if(n>m) swap(n,m);
        for(i=1;i<=n;i=r+1)
        {
            r=min(n/(n/i),m/(m/i));
            #define Fall 1ll*f[r]*invf[i-1]%mod 
            //注意!幂应该对mod-1取模 
            ans=1ll*ans*fpow(Fall,1ll*(n/i)*(m/i)%(mod-1))%mod;
            #undef Fall
        }
        printf("%d\n",ans);
    }
    return 0;
}

Blog来自PaperCloud，未经允许，请勿转载，TKS！

转载于:https://www.cnblogs.com/PaperCloud/p/10010260.html