容斥原理求gcd为k的数对个数
定义:
\(f(k)\)=\(\sum_{i}^{n}\)\(\sum_{j}^{m}\)\([gcd(i,j)=k]\)
即\(f(k)\)表示\(gcd(i,j)=k\)的数对个数
根据容斥原理,\(f(k)=\) 以k为公约数的数对个数\(-\)以k的倍数为最大公约数的数对个数
即\(f(k)=\lfloor{n/x}\rfloor*\lfloor{m/x}\rfloor-\sum f[k*i]\),其中\(k*i<=n\)
从\(min(n,m)\)开始倒着计算,即可用\(nlogn\)复杂度求出所有的\(f(k)\)
三道题:
1、洛谷P2398 GCD SUM
https://www.luogu.com.cn/problem/P2398
求出\(f(k)\)之后,答案就是\(\sum{k*f(k)}\)
#include<cstdio>
#define N 100001
long long f[N];
int main()
{
int n;
scanf("%d",&n);
for(int i=n;i;--i)
{
f[i]=1ll*(n/i)*(n/i);
for(int j=i+i;j<=n;j+=i) f[i]-=f[j];
}
long long ans=0;
for(int i=1;i<=n;++i) ans+=i*f[i];
printf("%lld",ans);
}
2、SDOI2008 仪仗队
https://www.luogu.com.cn/problem/P2158
他能看到它右面的一个和上面的一个,然后除去他所在的那一行一列,若\(gcd(i,j)=1\),则能看到,所以答案就是\(f(1)+2\)
\(n=1\)特判
#include<cstdio>
#define N 100001
long long f[N];
int main()
{
int n;
scanf("%d",&n);
--n;
if(!n)
{
printf("0");
return 0;
}
for(int i=n;i;--i)
{
f[i]=1ll*(n/i)*(n/i);
for(int j=i+i;j<=n;j+=i) f[i]-=f[j];
}
printf("%lld",f[1]+2);
}
3、NOI2010 能量采集
https://www.luogu.com.cn/problem/P1447
若\(gcd(i,j)=k\),则到点\((i,j)\)的过程中会经过k-1个点
因为
假设\(gcd(i,j)=k\),\(i^{'}=i/k\),\(j^{'}=j/k\),那么从\((0,0)\)到\((i,j)\)的过程中会经过\((i^{'},j^{'})\),\((2i^{'},2j^{'})\),\((3i^{'},3j^{'})\)……\(((k-1)i^{'},(k-1)j^{'})\)
所以答案就是\(\sum{f(k)*(2k-1)}\)
#include<cstdio>
#define N 100001
long long f[N];
int main()
{
int n,m;
scanf("%d%d",&n,&m);
int nm=n<m ? n : m;
for(int i=nm;i;--i)
{
f[i]=1ll*(n/i)*(m/i);
for(int j=i+i;j<=nm;j+=i) f[i]-=f[j];
}
long long ans=0;
for(int i=1;i<=n;++i) ans+=f[i]*(2*i-1);
printf("%lld",ans);
}