LOJ#2085. 「NOI2016」循环之美 莫比乌斯反演+杜教筛

求：$sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1][(j,k)=1]$

这个时候可以拆前面的，也可以拆后面的.

由于后面的 $k$ 是一个定值，考虑拆解后面的部分.

得：$sum_{d|k} mu(d) sum_{i=1}^{n} sum_{j=1}^{frac{m}{d}} [(i,jd)=1]$

$Rightarrow sum_{d|k} mu(d) sum_{i=1}^{n} sum_{j=1}^{frac{m}{d}} [(i,j)=1][(i,d)=1]$

令 $S(n,m,k)=sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1][(j,k)=1]$

则有 $S(n,m,k)=sum_{d|k} mu(d)S(frac{m}{d},n,d)$

其中边界条件是 $n,m leqslant 0$ 和 $k=1$.

$k=1$ 时求的是 $sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1]$

也就是 $sum_{d=1}^{min(n,m)} mu(d)frac{n}{d}frac{m}{d}$

这个需要用到杜教筛来计算 $mu$ 的前缀和.

$Rightarrow S(n)=frac{sum_{i=1}^{n}(f*g)(i)-sum_{i=2}^{n}S(frac{n}{i})g(i)}{g(1)}$

我们选 $g=I(x)$，则 $(f*g)(i)=[i=1]$.

则 $S(n)=1-sum_{i=2}^{n} S(frac{n}{i})$

时间复杂度不会计算，但是跑的挺快.

代码：

#include <map>
#include <cstdio>
#include <cmath>
#include <vector>
#include <cstring> 
#define N 20000008 
#define ll long long 
#define pb push_back  
#define setIO(s) freopen(s".in","r",stdin) 
using namespace std;
map<int,ll>ansmu;  
map<int,int>vised;    
vector<int>G[2003];  
int vis[N],mu[N],prime[N],sum[N],tot;  
void init() { 
	mu[1]=1;   
	for(int i=2;i<N;++i) {  
		if(!vis[i]) prime[++tot]=i,mu[i]=-1;  
		for(int j=1;j<=tot&&i*prime[j]<N;++j) { 
			vis[i*prime[j]]=1;    
			if(i%prime[j]) { 
				mu[i*prime[j]]=-mu[i];      
			} 
			else { 
				mu[i*prime[j]]=0; 
				break; 
			}
		}
	}   
	for(int i=1;i<N;++i) { 	
		sum[i]=sum[i-1]+mu[i];  
	}
	for(int i=1;i<2003;++i) { 	
		for(int j=1;j<=i;++j)    
			if((i%j)==0&&mu[j]) G[i].pb(j); 
	}
}  
ll gcd(ll a,ll b) { 
	return b?gcd(b,a%b):a; 
}   
ll getmu(int n) { 
	if(n<N) return sum[n];  
	if(vised[n]) return ansmu[n];      
	ll ans=0; 
	for(int i=2,j;i<=n;i=j+1) {   
		j=n/(n/i),ans+=(j-i+1)*getmu(n/i);    
	}	 
	vised[n]=1,ansmu[n]=1ll-ans;  
	return 1ll-ans;  
}
ll calc(int n,int m) {  
	if(n>m) swap(n,m);  
	ll ans=0; 
	for(int i=1,j;i<=n;i=j+1)  { 
		j=min(n/(n/i),m/(m/i));  
		ans+=1ll*(n/i)*(m/i)*(getmu(j)-getmu(i-1));  
	}   
	return ans;  
}
ll solve(int n,int m,int k) {  
	if(n<=0||m<=0) return 0;   
	if(k==1) {   
		return calc(n,m);  
	}  
	ll ans=0;
	for(int i=0;i<G[k].size();++i) {  
		int cur=G[k][i];    
		ans+=1ll*mu[cur]*solve(m/cur,n,cur);  
	}  
	return ans;  
}
int main() {
	// setIO("input");    	
	int n,m,k;  
	scanf("%d%d%d",&n,&m,&k);  	      
	ll ans=0; 
	init();          
	printf("%lld
",solve(n,m,k));    
	return 0; 
}