[洛谷P2257]YY的GCD 莫比乌斯反演+分块

题目描述

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

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

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

多组输入

输入输出格式

输入格式：

第一行一个整数T 表述数据组数

接下来T行，每行两个正整数，表示N, M

输出格式：

T行，每行一个整数表示第i组数据的结果

输入输出样例

输入样例#1： 

2
10 10
100 100

输出样例#1：

30
2791

说明

T = 10000

N, M <= 10000000

这里是一道莫比乌斯反演的板子题。虽然说是板子题，但你如果第一次做这种题目并且不知道解题套路的话是很难做出来的。

我们先将问题简化为 这个式子应该不需要过多赘述。如果是老玩家，那么直接设$f(d)$ 表示满足 $\mathrm{g}\mathrm{c}\mathrm{d}(x,y)=d$ 且 $1\le x\le n,1\le y\le m$这样的 $(x,y)$ 有多少对；设 $g(d)$ 表示满足 $d\mid \mathrm{g}\mathrm{c}\mathrm{d}(x,y)$ 且 $1\le x\le n,1\le y\le m$ 这样的 (x, y)

$(x,y)$ 有多少对。这个g(x)我们是能直接O(1)求出来的

同时它与f(d)还有莫比乌斯变化的关系

因此莫比乌斯反演得到的结果就是

回过头来看题目的要求是让我们求出所有质数p的f(p)之和，即

因为枚举pi的倍数不太方便变形，于是我们改为枚举pi的系数。

由于dpi的值是从1到min(n,m)的因此我们把dpi合成为一个T。

观察式子我们发现，后面的可以在求出莫比乌斯函数u以后顺便求出来。因此单次查询的时间复杂度可以做到O(n)

可是良（xian）心（e）的出题人一定不会这样轻易让你AC，所以要求你算出10000组数据的答案。这个时候分块就闪亮登场了。

我们可以发现这个式子的结果在很多情况下都是一样的，再推一推我们发现一共只有种不同的取值。于是可以将一样的值合并处理。这时候我们维护的前缀和，并把一样的一连串T利用乘法分配律求出即可。

AC代码（PS:标程为了更加清楚明白，常数比较大，所以开了O2优化依然比较卡，望大家见谅）：

#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
const int Maxn=10000005;
ll pri[Maxn/10],u[Maxn],mu[Maxn],sum[Maxn],k,t;
bool vis[Maxn];
inline ll mn(ll x,ll y){
    return x<y?x:y;
}
inline ll read(){//快速读入，卡常必备 
    ll x=0;char c=getchar();
    while(c<'0'||c>'9'){c=getchar();}
    while(c<='9'&&c>='0'){x=x*10;x=x+c-'0';c=getchar();}
    return x;
}
void prep(){//预处理 
    u[1]=1;
    for(int i=2;i<=Maxn;++i){//筛素数+求莫比乌斯函数 
        if(!vis[i]){pri[++k]=i;u[i]=-1;}
        for(int j=1;j<=k&&pri[j]*i<=Maxn;++j){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0)break;
            u[pri[j]*i]=-u[i];
        }
    }
    for(int i=1;i<=k;++i){
        for(int j=1;j*pri[i]<=Maxn;++j){
            mu[j*pri[i]]+=u[j];//求sigma(k|T)(u(T/k)) 
        } 
    }
    for(int i=1;i<=Maxn;++i){
        sum[i]=sum[i-1]+mu[i];//前缀和 
    }
}
void work(){
    ll n,m,x,ans=0,lim;
    n=read();m=read();lim=mn(n,m);
    for(int i=1;i<=lim;i=x+1){
        x=mn(n/(n/i),m/(m/i));//分块 
        ans+=(n/i)*(m/i)*(sum[x]-sum[i-1]);
    }
    printf("%lld\n",ans);
}
int main(){
    prep();
    scanf("%lld",&t);
    while(t--){
        work();
    }
    return 0;
}