转载至http://blog.csdn.net/clove_unique
题目链接
https://vjudge.net/problem/HYSBZ-2005
题解
首先证明对于某个点(x,y),k=gcd(x,y)-1:
设gcd(x,y)=t,令x=at,y=bt,那么在这条直线上的整数点可以表示为(a,b)(2a,2b)(3a,3b)……(x,y),由于不算x,y,则答案为gcd(x,y)-1
那么总损耗2k+1=2×gcd(x,y)-1。
我们最终要求的式子为:
∑i=1n∑j=1m(gcd(i,j)∗2−1)
=2∗∑i=1n∑j=1mgcd(i,j)−n∗m
那么我们只需要算出
∑i=1n∑j=1mgcd(i,j)
这个式子就可以了
推导如下:
∑i=1n∑j=1mgcd(i,j)
=∑i=1n∑j=1m∑d|gcd(i,j)ϕ(d)
=∑i=1n∑j=1m∑d=1n[d|i][d|j]ϕ(d)
=∑d=1n∑i=1n[d|i]∑j=1m[d|j]ϕ(d)
=∑d=1n⌊nd⌋⌊md⌋ϕ(d)
实际上
⌊nd⌋⌊md⌋
只有
(n√+m−−√)
个取值。
可以用分块来求。
需要预处理phi的前缀和。
#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
#define LL long long
const int max_n=1e5+5;
LL n,m,ans;
LL p[max_n],phi[max_n],prime[max_n];
inline void get_phi(){
phi[1]=1;
for (int i=2;i<=n;++i){
if (!p[i]){
prime[++prime[0]]=i;
phi[i]=i-1;
}
for (int j=1;j<=prime[0]&&i*prime[j]<=n;++j){
p[i*prime[j]]=1;
if (i%prime[j]==0){
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
phi[i]+=phi[i-1];
}
}
int main(){
scanf("%lld%lld",&n,&m);
if (n>m) swap(n,m);
get_phi();
for (LL i=1,j;i<=n;i=j+1){
j=min(n/(n/i),m/(m/i));
ans+=(LL)(phi[j]-phi[i-1])*(n/i)*(m/i);
}
printf("%lld\n",ans*2-n*m);
}
另外,本题还有一种程序更简便的解法。令f(k)为gcd(i, j) == k的数对的个数,令F(k)为gcd(i, j)为k的倍数的数对的个数。则F(k) = ∑f(i*k),其中 i*k <= min(n, m),另外也有F(k) = (n/k) * (m/k),所以f(k) = (n/k)*(m/k) - f(2*k) - f(3*k) - f(4*k) ... ... ,逆推出来即可。
#include<cstdio>
#include<cstring>
#include<string>
#include<cctype>
#include<iostream>
#include<set>
#include<map>
#include<cmath>
#include<sstream>
#include<vector>
#include<stack>
#include<queue>
#include<algorithm>
#define fin(a) freopen("a.txt","r",stdin)
#define fout(a) freopen("a.txt","w",stdout)
typedef long long LL;
using namespace std;
const int INF = 1e8 + 10;
const int maxn = 1e5 + 10;
LL f[maxn]; //f[i] : gcd(x, y) == i的数的个数
int main() {
int n, m;
scanf("%d%d", &n, &m);
int t = min(n, m);
LL ans = 0;
for(int i = t; i >= 1; i--) {
f[i] = (LL)(m/i)*(n/i);
for(int j = i+i; j <= t; j += i)
f[i] -= f[j];
ans += f[i] * (2*i-1);
}
printf("%lld\n", ans);
return 0;
}