• 【BZOJ3994】[SDOI2015] 约数个数和(莫比乌斯反演)


    点此看题面

    大致题意:\(d(x)\)\(x\)的约数个数,求\(\sum_{i=1}^N\sum_{j=1}^Md(i·j)\)

    莫比乌斯反演

    这是一道莫比乌斯反演题。

    一个重要的性质

    首先我们要先了解\(d(i·j)\)这个函数的性质:

    \[d(i·j)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1] \]

    证明: 我也不知道,应该就是枚举\(i\)\(j\)的约数,求出其中不互质的约数对个数,避免重复计算。

    一些定义

    按照莫比乌斯反演的常见套路,我们可以定义\(f(d)\)\(F(n)\)如下:

    \[f(d)=\sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)==d] \]

    \[F(n)=\sum_{n|d}f(d) \]

    然后由莫比乌斯反演的某些性质,我们可以得到下面这个式子:

    \[f(n)=\sum_{n|d}\mu(\lfloor\frac dn\rfloor)F(d) \]

    公式化简

    首先,题目中已经给出:

    \[answer=\sum_{i=1}^N\sum_{j=1}^Md(i·j) \]

    由于上面提到的性质,我们可以得到:

    \[answer=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}[gcd(x,y)==1] \]

    根据莫比乌斯函数\(\mu\)的性质\(\sum_{d|n}\mu(d)=[n==1]\),所以,我们可以将\(\mu\)代入,将原式变成这个样子:

    \[answer=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}\sum_{d|gcd(x,y)}\mu(d) \]

    由于这个式子难以操作,因此,我们可以将原式略作修改,改成对\(d\)进行枚举,变成这个样子:

    \[answer=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}\sum_{d=1}^{min(n,m)}\mu(d)*[d|gcd(x,y)] \]

    不难发现,\(\mu(d)\)的值是与\(i,j,x,y\)无关的,因此可以将其单独提出,就变成了这样:

    \[answer=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}[d|gcd(x,y)] \]

    然后,我们可以从枚举\(i,j\)及其约数\(x,y\),转变为直接枚举约数\(x,y\),然后将其贡献乘上约数倍数的个数(这应该还是比较好理解的),于是就有了下面这个式子:

    \[answer=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^N\sum_{y=1}^M[d|gcd(x,y)]\lfloor\frac Nx\rfloor\lfloor\frac My\rfloor \]

    考虑由\(d|gcd(x,y)\)可以得到\(d|x\)\(d|y\),即\(x\)\(y\)\(d\)的倍数所以我们就可以通过直接枚举\(d\)的倍数\(d·x\)\(d·y\)来取代枚举\(x,y\),从而消去\(d|gcd(x,y)\)这个限制:

    \[answer=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac Nd\rfloor}\sum_{y=1}^{\lfloor\frac Md\rfloor}\lfloor\frac N{d·x}\rfloor\lfloor\frac M{d·y}\rfloor \]

    最后,我们可以将原式稍作变动,得到下面这个式子:

    \[answer=(\sum_{d=1}^{min(n,m)}\mu(d))(\sum_{x=1}^{\lfloor\frac Nd\rfloor}\lfloor\frac N{d·x}\rfloor)(\sum_{y=1}^{\lfloor\frac Md\rfloor}\lfloor\frac M{d·y}\rfloor) \]

    求解答案

    不难想到对这个式子中的后两部分进行除法分块(貌似需要二次除法分块),然后对\(\mu(d)\)用前缀和预处理一下,就能做到单次询问\(O(\sqrt N)\)的时间复杂度,这样就能过了。

    代码

    #include<bits/stdc++.h>
    #define max(x,y) ((x)>(y)?(x):(y))
    #define min(x,y) ((x)<(y)?(x):(y))
    #define uint unsigned int
    #define LL long long
    #define ull unsigned long long
    #define swap(x,y) (x^=y,y^=x,x^=y)
    #define abs(x) ((x)<0?-(x):(x))
    #define INF 1e9
    #define Inc(x,y) ((x+=(y))>=MOD&&(x-=MOD))
    #define ten(x) (((x)<<3)+((x)<<1))
    #define N 50000
    using namespace std;
    int n,m;
    class FIO
    {
        private:
            #define Fsize 100000
            #define tc() (FinNow==FinEnd&&(FinEnd=(FinNow=Fin)+fread(Fin,1,Fsize,stdin),FinNow==FinEnd)?EOF:*FinNow++)
            #define pc(ch) (FoutSize<Fsize?Fout[FoutSize++]=ch:(fwrite(Fout,1,FoutSize,stdout),Fout[(FoutSize=0)++]=ch))
            int f,FoutSize,OutputTop;char ch,Fin[Fsize],*FinNow,*FinEnd,Fout[Fsize],OutputStack[Fsize];
        public:
            FIO() {FinNow=FinEnd=Fin;}
            inline void read(int &x) {x=0,f=1;while(!isdigit(ch=tc())) f=ch^'-'?1:-1;while(x=ten(x)+(ch&15),isdigit(ch=tc()));x*=f;}
            inline void read_char(char &x) {while(isspace(x=tc()));}
            inline void read_string(string &x) {x="";while(isspace(ch=tc()));while(x+=ch,!isspace(ch=tc())) if(!~ch) return;}
            inline void write(LL x) {if(!x) return (void)pc('0');if(x<0) pc('-'),x=-x;while(x) OutputStack[++OutputTop]=x%10+48,x/=10;while(OutputTop) pc(OutputStack[OutputTop]),--OutputTop;}
            inline void write_char(char x) {pc(x);}
            inline void write_string(string x) {register int i,len=x.length();for(i=0;i<len;++i) pc(x[i]);}
            inline void end() {fwrite(Fout,1,FoutSize,stdout);}
    }F;
    class Class_Mobius//莫比乌斯反演
    {
        private:
            int Prime_cnt,mu[N+5],Prime[N+5];bool IsNotPrime[N+5];
        public:
            int sum[N+5];LL g[N+5];
            Class_Mobius()//预处理
            {
                register int i,j,l,r;
                for(mu[1]=1,i=2;i<=N;++i)//求出莫比乌斯函数
                {
                    if(!IsNotPrime[i]) Prime[++Prime_cnt]=i,mu[i]=-1;
                    for(j=1;j<=Prime_cnt&&i*Prime[j]<=N;++j) 
                        if(IsNotPrime[i*Prime[j]]=true,i%Prime[j]) mu[i*Prime[j]]=-mu[i];else break;
                }
                for(i=1;i<=N;++i) sum[i]=sum[i-1]+mu[i];//求出前缀和
                for(i=1;i<=N;++i) for(l=1;l<=i;l=r+1) r=i/(i/l),g[i]+=1LL*(r-l+1)*(i/l);//第一次除法分块
            }
    }Mobius;
    int main()
    {
        register int l,r,lim,T;register LL ans;F.read(T);
        while(T--)
        {
            for(F.read(n),F.read(m),ans=0,l=1,lim=min(n,m);l<=lim;l=r+1)//第二次除法分块
                r=min(n/(n/l),m/(m/l)),ans+=1LL*(Mobius.sum[r]-Mobius.sum[l-1])*Mobius.g[n/l]*Mobius.g[m/l];
            F.write(ans),F.write_char('\n');
        }
        return F.end(),0;
    }
    
    败得义无反顾,弱得一无是处
  • 相关阅读:
    文件重名问题
    文件上传
    回顾IO流
    Freemarker
    中文乱码问题
    Filter(过滤器)
    Ajax
    jQuery
    普华操作系统,开机无法进入桌面程序; 解决多次source /etc/profile的烦恼
    C++ 文件类型判别错误,将目录文件识别为普通文件
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3994.html
Copyright © 2020-2023  润新知