Luogu 2257

Luogu 2257 - YY的GCD

题意

求下式的值

[sum_{pin prime}sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=p] ]

限制

(Case=10^4, n,mleq 10^7)

思路

令(nleq m)

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

[sum_{pin prime}sum_{i=1}^{lfloorfrac n p floor}sum_{j=1}^{lfloorfrac m p floor}[gcd(i,j)=1] ]

由莫比乌斯反演得

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

所以原式可得

[sum_{pin prime}sum_{i=1}^{lfloorfrac n p floor}sum_{j=1}^{lfloorfrac m p floor}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)的倍数有(lfloorfrac a b floor)个

所以原式便等同于

[sum_{pin prime}sum_{d=1}^{lfloorfrac n p floor}mu(d)lfloorfrac{lfloorfrac n p floor}d floorlfloorfrac{lfloorfrac m p floor}d floor\ =sum_{pin prime}sum_{d=1}^{lfloorfrac n p floor}mu(d)lfloorfrac n {pd} floorlfloorfrac m {pd} floor ]

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

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

那么原式将会变成

[sum_{pin prime}sum_{d=1}^{lfloorfrac n p floor}mu(lfloorfrac x p floor)lfloorfrac n x floorlfloorfrac m x floor\ =sum_{x=1}^nlfloorfrac n x floorlfloorfrac m x floorsum_{pin prime, p|x}mu(lfloorfrac x p floor) ]

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

再对其求次前缀和，对原式进行分块求解，即可在(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<<'
';
}
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    init(N);
    int T;
    cin>>T;
    while(T--)
        solve();
    return 0;
}