P2257 YY的GCD

题目描述

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

给定 N, M ，求 1 \leq x \leq N，1 \leq y \leq M 且 \gcd(x, y) 为质数的 (x, y) 有多少对。

多组询问

T=10^4，N, M \leq 10^7

Solution

题目求：(N\leq M)
\sum_{p\in prime}\sum_{x=1}^{N}\sum_{y=1}^M[(x,y)=p]

提 p 出来

\sum_{p\in prime}\sum_{x=1}^{\frac{N}{p}}\sum_{y=1}^{\frac{M}{p}}[(x,y)=1]

令 N’=\frac{N}{p} ，M’=\frac{M}{p}

令 f(d)=\sum_{p\in prime}\sum_{x=1}^{N’}\sum_{y=1}^{M’}[(x,y)=d]

F(n)=\sum\limits_{n|d}f(d)=f(d)=\sum_{p\in prime}\sum_{x=1}^{N’}\sum_{y=1}^{M’}[n|(x,y)]

即有

F(n)=\sum_{p\in prime}\lfloor{\frac{N}{np}\rfloor}\lfloor\frac{M}{np}\rfloor

利用莫比乌斯反演：

f(d)=\sum_{d|n}\mu(n)F(n)

f(d)=\sum_{d|n}\mu(n)\sum_{p\in prime}\lfloor{\frac{N}{np}\rfloor}\lfloor\frac{M}{np}\rfloor

令 T=np

f(1)=\sum_{n=1}^N\sum_{p\in prime}\mu(\frac{T}{p})\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor

改变枚举方式，枚举 T

f(1)=\sum_{T=1}^N\sum_{p\in prime , p|T}\mu(\frac{T}{p})\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor=\sum_{p\in prime ,p|T}\mu(\frac{T}{p})\sum_{T=1}^N\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor

对于每一个 T ，我们预处理出 \sum_{p\in prime ,p|T}\mu(\frac{T}{p})

然后利用整除分块求解即可

#include<bits/stdc++.h>
#define ll long long
const int N = 1e7+10;
bool vis[N];
int prime[N];
int Mob[N];
int g[N],pre[N];
ll h[N];
void Mobius_sieve(){
    int cnt = 0;
    vis[1] = 1;
    Mob[1] = 1;
    for(int i = 2; i < N; i++){
        if(!vis[i])
            prime[cnt++] = i, Mob[i] = - 1 , g[i] = 1;
        for(int j = 0; j < cnt && 1LL * prime[j] * i < N; j++){
            vis[prime[j] * i] = 1;
            g[prime[j]*i] = g[i] + 1;
            Mob[i * prime[j]] = (i % prime[j] ? -Mob[i]: 0);
            if(i % prime[j] == 0)
                break;
        }
    }
    for(int i=2;i<N;++i){
        if(!vis[i]) pre[i] = 1;
        if(Mob[i]!=0){
            for(int j=0;prime[j]*i<N;++j){
                pre[prime[j]*i] += Mob[i]; 
            }
        }
    }
    for(int i=1;i<N;++i)
        pre[i]+=pre[i-1];
}
int main(){
    Mobius_sieve();
    int t = 0 ;
    scanf("%d",&t);
    while(t--){
        int n,m;
        long long ans = 0 ;
        scanf("%d%d",&n,&m);
        if(n>m) std::swap(n,m);
        for(int l=1,r=0;l<=n;l=r+1){
            r = std::min(n,std::min(n/(n/l),m/(m/l)));
            ans += 1ll*(pre[r]-pre[l-1])*(n/l)*(m/l);
        }
        printf("%lld\n",ans);
    }
}