Luogu P2257 YY的GCD

题意

(T) 组数据，每组数据给定 (n,m)，求 (sumlimits_{pin P}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=p])。（这里 (P) 是质数集）

( exttt{Data Range:}T=10^4,1leq n,mleq 10^7)

题解

第一个求和符号里面的和式很平凡，直接写出来大概是这样：

[sumlimits_{pin P}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=p]=sumlimits_{pin P}sumlimits_{d=1}^{leftlfloorfrac{n}{p} ight floor}mu(d)leftlfloorfrac{n}{dp} ight floorleftlfloorfrac{m}{dp} ight floor ]

然后我们考虑将枚举 (d) 改成枚举 (T=dp)，那么有

[sumlimits_{pin P}sumlimits_{d=1}^{leftlfloorfrac{n}{p} ight floor}mu(d)leftlfloorfrac{n}{dp} ight floorleftlfloorfrac{m}{dp} ight floor=sumlimits_{pin P}sumlimits_{T=1}^{n}muleft(frac{T}{p} ight)leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floor ]

注意到我们可以交换一下求和变量

[sumlimits_{pin P}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=p]=sumlimits_{T=1}^{n}leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floorsumlimits_{pmid T,pin P}muleft(frac{T}{p} ight) ]

注意到在筛出 (mu) 之后可以考虑通过枚举来预处理出 (sumlimits_{pmid T,pin P}muleft(frac{T}{p} ight))，然后整除分块就做完了。

代码

#include<bits/stdc++.h>
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=1e7+51;
ll test,ptot,n,m;
ll np[MAXN],prime[MAXN],mu[MAXN];
li f[MAXN],prefix[MAXN];
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
}
inline void sieve(ll limit)
{
    np[1]=mu[1]=1;
    for(register int i=2;i<=limit;i++)
    {
        if(!np[i])
        {
            prime[++ptot]=i,mu[i]=-1;
        }
        for(register int j=1;i*prime[j]<=limit;j++)
        {
            np[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}
inline li calc(ll n,ll m)
{
    li res=0;
    for(register int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        res+=(prefix[r]-prefix[l-1])*(n/l)*(m/l);
    }
    return res;
}
int main()
{
    test=read(),sieve(1e7+10);
    for(register int i=1;i<=ptot;i++)
    {
        for(register int j=1;prime[i]*j<=1e7;j++)
        {
            f[prime[i]*j]+=mu[j];
        }
    }
    for(register int i=1;i<=1e7;i++)
    {
        prefix[i]=prefix[i-1]+f[i];
    }
    for(register int i=0;i<test;i++)
    {
        n=read(),m=read(),printf("%lld
",calc(n,m));
    }
}