原文链接www.cnblogs.com/zhouzhendong/p/UOJ221.html
题解
首先把题目转化为求
[sum_{x=1}^n sum_{y=1}^m [gcd(x,y) = 1] [ gcd(y,k) = 1]
]
推式子：
[sum_{x=1}^n sum_{y=1}^m [gcd(x,y) = 1] [ gcd(y,k) = 1]\ = sum_{i = 1}^{min(n,m)} mu(i) sum_{x=1}^{left lfloor frac n i
ight
floor } sum_{y=1}^{left lfloor frac m i
ight
floor } [gcd(iy,k) = 1]\ = sum_{i = 1}^{min(n,m)}left lfloor frac n i
ight
floor mu(i) [gcd(i,k) = 1] sum_{y=1}^{left lfloor frac m i
ight
floor } [gcd(y,k) = 1]
]
整除分块一下，变成求
[sum_{i = 1} ^ n mu(i) [gcd(i,k) = 1]
]
推式子：
[sum_{i = 1} ^ n mu(i) [gcd(i,k) = 1] \ = sum_{i = 1} ^ n mu(i) - sum_{i = 1} ^ n mu(i) [gcd(i,k) > 1] \ = sum_{i = 1} ^ nmu(i) - sum_{d>1,d|k} sum_{i = 1} ^ {lfloor frac n d
floor} [gcd(id,k) = d] mu(id)
]
其中
[[gcd(id,k) = d] mu(id) \ = [gcd(id,k) = d] mu(i) mu(d) [gcd(i,d) = 1] \ = mu(d) [gcd(i,k) = 1] mu(i)
]
所以，原式 =
[sum_{i = 1} ^ nmu(i) - sum_{d>1,d|k} mu(d) sum_{i = 1} ^ {lfloor frac n d
floor} mu(i) [gcd(i,k) = 1]
]
大力杜教筛即可。
代码
#include <bits/stdc++.h>
#define clr(x) memset(x,0,sizeof x)
#define For(i,a,b) for (int i=(a);i<=(b);i++)
#define Fod(i,b,a) for (int i=(b);i>=(a);i--)
#define fi first
#define se second
#define pb(x) push_back(x)
#define mp(x,y) make_pair(x,y)
#define outval(x) cerr<<#x" = "<<x<<endl
#define outtag(x) cerr<<"---------------"#x"---------------"<<endl
#define outarr(a,L,R) cerr<<#a"["<<L<<".."<<R<<"] = ";
For(_x,L,R)cerr<<a[_x]<<" ";cerr<<endl;
using namespace std;
typedef long long LL;
LL read(){
LL x=0,f=0;
char ch=getchar();
while (!isdigit(ch))
f|=ch=='-',ch=getchar();
while (isdigit(ch))
x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
return f?-x:x;
}
const int N=1000005;
int gcd(int a,int b){
return b?gcd(b,a%b):a;
}
int n,m,k;
vector <int> v;
int p[N],pcnt,u[N];
void getprime(int n){
static int vis[N];
clr(vis),pcnt=0;
u[1]=1;
For(i,2,n){
if (!vis[i])
p[++pcnt]=i,u[i]=-1;
for (int j=1;j<=pcnt&&i*p[j]<=n;j++){
vis[i*p[j]]=1;
if (i%p[j])
u[i*p[j]]=-u[i];
else {
u[i*p[j]]=0;
break;
}
}
}
}
unordered_map <int,int> fu,U,fu2,U2;
int GetU(int n){
if (fu[n])
return U[n];
int res=1;
for (int i=2,j;i<=n;i=j+1){
j=n/(n/i);
res-=GetU(n/i)*(j-i+1);
}
fu[n]=1;
return U[n]=res;
}
int GetU2(int n){
if (fu2[n])
return U2[n];
int res=GetU(n);
for (auto i : v)
if (i>1)
res-=(LL)u[i]*GetU2(n/i);
fu2[n]=1;
return U2[n]=res;
}
int calcK(int n){
int ans=0;
for (auto i : v)
ans+=u[i]*(n/i);
return ans;
}
int main(){
n=read(),m=read(),k=read();
getprime(1000000);
v.clear();
For(i,1,k)
if (k%i==0&&u[i])
v.pb(i);
U[0]=0,U2[0]=0,fu[0]=fu2[0]=1;
For(i,1,1000000){
fu[i]=1,U[i]=U[i-1]+u[i];
fu2[i]=1,U2[i]=U2[i-1]+(gcd(i,k)==1?u[i]:0);
}
LL ans=0;
for (int i=1,j;i<=n&&i<=m;i=j+1){
j=min(n/(n/i),m/(m/i));
ans+=(LL)(GetU2(j)-GetU2(i-1))*(n/i)*calcK(m/i);
}
cout<<ans<<endl;
return 0;
}