YY的GCD
题目描述
神犇YY虐完数论后给傻×kAc出了一题
给定(N), (M) ,求(1 le x le N,1 le y le M)且(gcd(x, y))为质数的((x, y))有多少对
kAc这种傻×必然不会了，于是向你来请教……
多组输入
输入输出格式
输入格式：
第一行一个整数(T)表述数据组数
接下来T行，每行两个正整数，表示(N,M)
输出格式：
(T)行，每行一个整数表示第(i)组数据的结果
说明
(T = 10000)
(N, M le 10000000)
Solution
思路：按反演的套路把式子列出来以后，通过想办法操作求和，把好算的放在一起，(mu)一个人呆。
以下(f(x)=sum_{i=1}^asum_{j=1}^b[gcd(a,b)=x])，(g(x)=lfloorfrac{a}{x} floor lfloor frac{b}{x} floor)，(p)代表质数集合
[sum_{i=1}^asum_{j=1}^b[gcd(a,b) in p]
]
[=sum_{n in p}f(n)
]
[=sum_{n in p}sum_{n|d}^{min(a,b)}mu(frac{d}{n})g(d)
]
[=sum_{n in p}sum_{k=1}^{min(lfloor frac{a}{n}
floor, lfloor frac{b}{n}
floor)}mu(k)lfloorfrac{a}{kn}
floor lfloorfrac{b}{kn}
floor
]
现在我们不想让(mu)和别人呆在一起，因为两个求和之间也没夹个什么东西，不妨把它们换一换，直接换肯定不行，得找一波操作。
令(T=kn)，则有(k=frac{T}{n})，上式就等于
[=sum_{T=1}^{min(a,b)}sum_{n|T,nin p}mu(frac{T}{n})lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floor
]
[=sum_{T=1}^{min(a,b)}lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floorsum_{n|T,nin p}mu(frac{T}{n})
]
后面那个感觉预处理还挺慢的，但根据调和级数不会慢于(O(nlnn))，而且跑一下发现才(3e7)多运算。
预处理完后面那个就可以除法分块回答询问了。
Code:
#include <cstdio>
#define ll long long
const int N=1e7;
int mu[N+10],f[N+10],pri[N+10],ispri[N+10],cnt,a,b,T;
void init()
{
mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!ispri[i])
{
mu[i]=-1;
pri[++cnt]=i;
}
for(int j=1;j<=cnt&&pri[j]*i<=N;j++)
{
ispri[pri[j]*i]=1;
if(i%pri[j]==0) break;
else mu[pri[j]*i]=-mu[i];
}
}
for(int i=1;i<=cnt;i++)
for(int j=pri[i];j<=N;j+=pri[i])
f[j]+=mu[j/pri[i]];
for(int i=1;i<=N;i++)
f[i]+=f[i-1];
}
int min(int x,int y){return x<y?x:y;}
int main()
{
init();
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&a,&b);
ll ans=0;
for(int l=1,r;l<=min(a,b);l=r+1)
{
r=min(a/(a/l),b/(b/l));
ans+=1ll*(a/l)*(b/l)*(f[r]-f[l-1]);
}
printf("%lld
",ans);
}
return 0;
}
2018.10.19