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\)

又已知在\(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;
}