Luogu 2257 - YY的GCD （莫比乌斯反演、分块）

Luogu 2257 - YY的GCD

题意

求下式的值

\[\sum_{p\in prime}\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p] \]

限制

\(Case=10^4,\ n,m\leq 10^7\)

思路

令\(n\leq m\)

根据\(gcd(i,j)=k\iff gcd(\frac i k,\frac j k)=1\)，化简上式可得

\[\sum_{p\in prime}\sum_{i=1}^{\lfloor\frac n p\rfloor}\sum_{j=1}^{\lfloor\frac m p\rfloor}[\gcd(i,j)=1] \]

由莫比乌斯反演得

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

所以原式可得

\[\sum_{p\in prime}\sum_{i=1}^{\lfloor\frac n p\rfloor}\sum_{j=1}^{\lfloor\frac m p\rfloor}\sum_{d|\gcd(i,j)}\mu(d) \]

转化为枚举\(\gcd\)的因子\(d\)，由于\(i,j\)一定是\(d\)的倍数，所以不妨将\(i,j\)缩小\(d\)倍

设\(f(i,j)=\sum_{d|\gcd(i,j)}\mu(d)\)，可以知道对于某个\(\mu(d)\)，需要累加的次数等同于\(i\)范围内\(d\)倍数的个数与\(j\)范围内\(d\)倍数的个数之积

又已知在\(1\)~\(a\)范围内\(b\)的倍数有\(\lfloor\frac a b\rfloor\)个

所以原式便等同于

\[\sum_{p\in prime}\sum_{d=1}^{\lfloor\frac n p\rfloor}\mu(d)\lfloor\frac{\lfloor\frac n p\rfloor}d\rfloor\lfloor\frac{\lfloor\frac m p\rfloor}d\rfloor\\ =\sum_{p\in prime}\sum_{d=1}^{\lfloor\frac n p\rfloor}\mu(d)\lfloor\frac n {pd}\rfloor\lfloor\frac m {pd}\rfloor \]

此时直接处理复杂度是\(O(n\sqrt n)\)的

但我们可以先设\(x=pd\)，将其看作一个整体进行枚举

那么原式将会变成

\[\sum_{p\in prime}\sum_{d=1}^{\lfloor\frac n p\rfloor}\mu(\lfloor\frac x p\rfloor)\lfloor\frac n x\rfloor\lfloor\frac m x\rfloor\\ =\sum_{x=1}^n\lfloor\frac n x\rfloor\lfloor\frac m x\rfloor\sum_{p\in prime,\ p|x}\mu(\lfloor\frac x p \rfloor) \]

于是对于\(\sum_{p\in prime,\ p|x}\mu(\lfloor\frac x p \rfloor)\)，可以在求解莫比乌斯函数的线性筛里顺带求出

再对其求次前缀和，对原式进行分块求解，即可在\(O(\sqrt n)\)内求出答案

代码

Case Max (2840ms/4000ms)

O2 Case Max(2040ms/4000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=10000000;

int primcnt;
ll mu[N+50],prim[N+50];
ll f[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    primcnt=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[primcnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<primcnt;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
    for(int j=0;j<primcnt;j++)
        for(int i=1;i*prim[j]<=n;i++)
            f[i*prim[j]]+=mu[i];
    for(int i=2;i<=n;i++) //转为前缀和
        f[i]+=f[i-1];
}

void solve()
{
    int n,m;
    cin>>n>>m;
    if(n>m)
        swap(n,m);
    ll ans=0;
    for(int i=1;i<=n;)
    {
        int nxt=min(n/(n/i),m/(m/i));
        ans+=(f[nxt]-f[i-1])*(n/i)*(m/i);
        i=nxt+1; //分块跳转
    }
    cout<<ans<<'\n';
}
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    init(N);
    int T;
    cin>>T;
    while(T--)
        solve();
    return 0;
}