• 【洛谷2257】YY的GCD(莫比乌斯反演)


    点此看题面

    大致题意:\(\sum_{x=1}^N\sum_{y=1}^MIsPrime(gcd(x,y))\)

    莫比乌斯反演

    听说此题是莫比乌斯反演入门题?

    一些定义

    首先,我们可以定义\(f(d)\)\(F(d)\)如下:

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

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

    通过定义,不难发现:

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

    由于莫比乌斯反演的某些性质,我们又可以得到:

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

    公式化简

    首先,我们应该不难想到:

    \[answer=\sum_{IsPrime(p)}\sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)==p] \]

    然后就是一波化简。

    应该挺容易看出,由于\(f(p)\)的定义,上面的式子其实就相当于下面这个式子:

    \[answer=\sum_{IsPrime(p)}f(p) \]

    然后是莫比乌斯反演

    \[answer=\sum_{IsPrime(p)}\sum_{p|d}\mu(\lfloor\frac dp\rfloor)F(d) \]

    但是,这样有点难以处理。

    于是,我们改成枚举\(\lfloor\frac dp\rfloor\),于是原式就变成了这样:

    \[answer=\sum_{IsPrime(p)}\sum_{d=1}^{min(\lfloor\frac Np\rfloor,\lfloor\frac Mp\rfloor)}\mu(d)F(d·p) \]

    \(F(n)=\lfloor\frac Nn\rfloor\lfloor\frac Mn\rfloor\)代入进一步化简,可以得到:

    \[answer=\sum_{IsPrime(p)}\sum_{d=1}^{min(\lfloor\frac Np\rfloor,\lfloor\frac Mp\rfloor)}\mu(d)\lfloor\frac N{d·p}\rfloor\lfloor\frac M{d·p}\rfloor \]

    如果我们用\(G\)来表示\(d·p\),则\(d=\frac Gp\),原式就变成了这个样子:

    \[answer=\sum_{G=1}^{min(N,M)}\sum_{IsPrime(p),p|G}\mu(\frac Gp)\lfloor\frac NG\rfloor\lfloor\frac MG\rfloor \]

    通过乘法交换律乘法结合律,我们可以再一步转化,得:

    \[answer=\sum_{G=1}^{min(N,M)}\lfloor\frac NG\rfloor\lfloor\frac MG\rfloor(\sum_{IsPrime(p),p|G}\mu(\frac Gp)) \]

    然后就可以\(O(n)\)求解了。

    \(But\),多组数据,\(T\le 10000\)... ...

    求解答案

    首先,我们用定义一个\(g(n)\),它的定义如下:

    \[g(n)=\sum_{IsPrime(p),p|n}\mu(\frac np) \]

    于是,我们便能将上面的式子进一步转化:

    \[answer=\sum_{G=1}^{min(N,M)}\lfloor\frac NG\rfloor\lfloor\frac MG\rfloor g(n) \]

    然后,我们可以用除法分块

    不难发现,在一定范围内\(\lfloor\frac Ni\rfloor\)的值是保持不变的(\(\lfloor\frac Mi\rfloor\)同理),则\(\lfloor\frac NG\rfloor\lfloor\frac MG\rfloor\)其实最多只有\(\sqrt N+\sqrt M\),而对于\(\lfloor\frac NG\rfloor\lfloor\frac MG\rfloor\)相同的值,我们可以一起求解,于是就能想到用\(sum_i\)来表示\(\sum_{i=1}^i g(i)\),这样就能快速求解了。

    代码

    #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 10000000
    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:
            LL sum[N+5];
            Class_Mobius()//预处理
            {
                register int i,j;
                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(j=1;j<=Prime_cnt;++j) for(i=Prime[j];i<=N;i+=Prime[j]) sum[i]+=mu[i/Prime[j]];//计算g(i)
                for(i=1;i<=N;++i) sum[i]+=sum[i-1];//求前缀和
            }
    }Mobius;
    int main()
    {
        register int i,nxt,T;register LL ans;F.read(T);
        while(T--)
        {
        	if(F.read(n),F.read(m),n>m) swap(n,m);
        	for(ans=0,i=1;i<=n;i=nxt+1) nxt=min(n/(n/i),m/(m/i)),ans+=1LL*(n/i)*(m/i)*(Mobius.sum[nxt]-Mobius.sum[i-1]);//除法分块
        	F.write(ans),F.write_char('\n');//输出答案
        }
        return F.end(),0;
    }
    
    败得义无反顾,弱得一无是处
  • 相关阅读:
    openwrt编译出错处理记录
    转:小白编译openwrt固件教程
    转:openwrt中luci学习笔记
    (转)内存和外设的统一编址及独立编址
    (转)无法打开C盘,提示"本次操作由于这台计算机的限制被取消,请与你的管理员联系”
    (转)c# 解析JSON的几种办法
    计算机中断(转贴,很经典的解释)
    转:用WCAT进行IIS压力测试
    转:windows命令行下如何查看磁盘空间大小
    转:如何在32位程序中突破地址空间4G的限制
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/Luogu2257.html
Copyright © 2020-2023  润新知