UOJ#221. 【NOI2016】循环之美 数论,杜教筛

原文链接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;
}