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)
Solution
根据shenlao告诉我的一个推题目贼爽的公式
[sum_{d|gcd(x,y)}mu(d)=[gcd(x,y)=1]
]
原题意,求
[Ans=sum_{pin prime}sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=p]
]
然后就开始上瘾痛苦的推结论
[Ans=sum_{pin prime}sum_{i=1}^{lfloorfrac{n}{p}
floor}sum_{j=1}^{lfloorfrac{m}{p}
floor}[gcd(i,j)=1]
]
[Ans=sum_{pin prime}sum_{i=1}^{lfloorfrac{n}{p}
floor}sum_{j=1}^{lfloorfrac{m}{p}
floor}sum_{d|gcd(i,j)}mu(d)
]
更换枚举项,由枚举gcd(i,j)的约数改为直接枚举d,这中间对于i有(frac{n}{p})个数,其中d的倍数的个数为(frac {n}{d imes p}),对于j同理,相乘之后就保证了(d|gcd(i,j))这个条件
[Ans=sum_{pin prime}sum _{d=1}^{min({lfloorfrac{n}{p}
floor},{lfloorfrac{m}{p}
floor})}mu(d){lfloorfrac{n}{d imes p}
floor}{lfloorfrac{m}{d imes p}
floor}
]
令(d imes p=C)
[Ans=sum_{C=1}^{min(n,m)}sum_{p|C,pin prime}mu(frac {C}{p}){lfloorfrac{n}{C}
floor}{lfloorfrac{m}{C}
floor}
]
[Ans=sum_{C=1}^{min(n,m)}{lfloorfrac{n}{C}
floor}{lfloorfrac{m}{C}
floor}sum_{p|C,pin prime}mu(frac {C}{p})
]
至此公式推导结束,(O(n))就可以过了,可是由于多组询问,所以我们预处理一下.然后整除分块(O(sqrt n))做,不会的可以看本蒟蒻的博客整除分块
Code
#include<bits/stdc++.h>
#define rg register
#define il inline
#define lol long long
#define NN 10000000
#define Min(a,b) (a)<(b)?(a):(b)
using namespace std;
const int N=1e7+10;
void in(int &ans) {
ans=0; char i=getchar();
while(i<'0' || i>'9') i=getchar();
while(i>='0' && i<='9') ans=(ans<<1)+(ans<<3)+i-'0',i=getchar();
}
int n,m,T,tot;
int mu[N],prime[N];
lol f[N];
bool vis[N];
void print(lol x)
{
if(x<0) putchar('-'),x=-x;
if(x>9) print(x/10);
putchar(x%10+'0');
}
il void init() {
mu[1]=1;
for(rg int i=2;i<=NN;i++) {
if(!vis[i]) prime[++tot]=i,mu[i]=-1;
for(rg int j=1;j<=tot && i*prime[j]<=NN;j++) {
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
else mu[i*prime[j]]=-mu[i];
}
}
for(rg int i=1;i<=tot;i++)//注意枚举质数在前
for(rg int j=1;j*prime[i]<=NN;j++)//常数在后
f[prime[i]*j]+=1ll*mu[j];//对后面的数有贡献的是常数而不是质数
for(rg int i=1;i<=NN;i++) f[i]+=f[i-1];
}
int main()
{
in(T); init();
while(T--) {
in(n),in(m);
lol ans=0; if(n>m) swap(n,m);
for(rg int l=1,r;l<=n;l=r+1) {
r=Min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*(f[r]-f[l-1]);
}
print(ans);putchar('
');
}
return 0;
}