神犇YY虐完数论后给傻×kAc出了一题

给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对

kAc这种傻×必然不会了,于是向你来请教……

链接🔗

「Luogu-P2257」YY的GCD

题解

这题公式巨长,理解了好久才搞懂

先列出我们所需要求的东西是什么
pprimei=1nj=1m[gcd(i,j)=p]\sum_{p \in prime} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j) = p]

首先定义两个函数
f(d)=i=1nj=1m[gcd(i,j)=d]f(d) = \sum_{i=1}^{n} \sum_{j=1}^{m} [gcd(i,j) = d]
g(d)=i=1nj=1m[dgcd(i,j)]g(d) = \sum_{i=1}^{n} \sum_{j=1}^{m} [d | gcd(i,j)]

根据 g(d)g(d) 的值根据整除定义,可以得出 ndmd\left \lfloor \frac{n}{d} \right \rfloor \left \lfloor \frac{m}{d} \right \rfloor

根据 g(d)g(d) 的定义,可以得出 g(d)=ndf(d)g(d) = \sum_{n | d}f(d)

然后通过莫比乌斯反演可以得到 f(n)=ndμ(dn)g(d)f(n) = \sum_{n|d} \mu (\left \lfloor \frac{d}{n} \right \rfloor) g(d)


接下来 我们就开始对原公式进行化简

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

ans=pprimef(p)ans = \sum_{p \in prime} f(p)

ans=pprimepdμ(dp)g(d)ans = \sum_{p \in prime} \sum_{p|d} \mu (\left \lfloor \frac{d}{p} \right \rfloor) g(d)

这里设 n<mn < m ,然后将枚举项设为dd

ans=pprimednpμ(d)g(dp)ans = \sum_{p \in prime} \sum_{d}^{\left \lfloor \frac{n}{p} \right \rfloor} \mu (d) g(dp)

然后把g(dp)g(dp) 拆开来

ans=pprimednpμ(d)ndpmdpans = \sum_{p \in prime} \sum_{d}^{\left \lfloor \frac{n}{p} \right \rfloor} \mu (d) \left \lfloor \frac{n}{dp} \right \rfloor \left \lfloor \frac{m}{dp} \right \rfloor

dp=t令 dp = t

ans=pprimednpμ(d)ntmtans = \sum_{p \in prime} \sum_{d}^{\left \lfloor \frac{n}{p} \right \rfloor} \mu (d) \left \lfloor \frac{n}{t} \right \rfloor \left \lfloor \frac{m}{t} \right \rfloor

因为 dp=tdp = t ,所以后面部分提出来

t=1nntmtpt,pprimeμ(tp)\sum_{t = 1}^n \left \lfloor \frac{n}{t} \right \rfloor \left \lfloor \frac{m}{t} \right \rfloor \sum_{p|t,p \in prime} \mu(\frac{t}{p})

公式到这就化简完毕了

很明显,前面部分可以用数论分块做到 O((n))O( \sqrt(n) ) ,我们来讨论一下后面部分

μ\mu是可以线性筛出来的,对于任意一个质数pp,我们可以也可以通过筛出 pt,pprimeμ(tp)\sum_{p|t,p \in prime} \mu(\frac{t}{p}) 时间复杂度大概为 O(n loglog n)O(n \ log log \ n)

最后做个前缀和,就可以配合前面的式子进行数论分块了

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxd = 1e7+10;
int m,n,ans,k,t;
int mu[maxd], prime[maxd],summu[maxd],cnt,f[maxd];
bool vis[maxd];
void init() 
{
    mu[1] = 1;
    for(int i=2;i<maxd;i++)
    {
        if(!vis[i]) prime[++cnt] = i, mu[i] = -1;
        for(int j=1;j<=cnt && i*prime[j] <= maxd;j++)
        {
            vis[i*prime[j]] = 1;
            if(i%prime[j] == 0) break;
            mu[i*prime[j]] = -mu[i]; 
        }
    }
    for(int i=1;i<=cnt;i++)
        for(int j=1;prime[i]*j<maxd;j++)
            f[j*prime[i]] += mu[j];
		//前缀和
    for(int i=1;i<maxd;i++) summu[i] = summu[i-1] + f[i];
}
long long sum(int m,int n)  //数论分块
{
    long long t = min(n, m), ans = 0;
    for (int i = 1, j = 1; i <= t; i = j + 1) 
    {
        j = min(n / ( n / i), m / (m / i));
        ans += (summu[j] - summu[i - 1]) * (long long)(n / i) * (m / j);
    }
    return ans;
}
int main()
{
    // freopen("a.in","r",stdin);
    // freopen("k.out","w",stdout);
    init();scanf("%d",&t);
    for(int i=1;i<=t;i++)
    {
        scanf("%d %d",&n,&m);
        printf("%lld\n",sum(n,m));
    }
    return 0;
}