P3455 [POI2007]ZAP-Queries

题目

求：

\[\large \sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k] \]

其中 \(n,m,k\le 5\times 10^4\) ，\(t\le 5\times 10^4\) 组数据。

P3455 [POI2007]ZAP-Queries

分析

一般有这样的套路，设函数 \(g(d)\) 为：

\[\large g(d)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d] \]

这时我们可以观察我们的函数 \(f\) ：

\[\large f(d)=\sum\limits_{d\mid n}g(n)=\sum\limits_{d\mid n}\sum\limits_{i=1}^N\sum\limits_{j=1}^M[\gcd(i,j)=d] \]

发现这时的 \(f\) 函数性质非常好，就相当于取两个数使得他们都是 \(d\) 的倍数，可以直接写出：

\[\large f(d)=\lfloor\dfrac Nd\rfloor\lfloor \dfrac Md\rfloor \]

然后我们尝试通过莫比乌斯反演来求出 \(g\) 就是利用 \(f\) ：

\[\large g(n)=\sum\limits_{n|d} \mu(\frac{d}{n}) f(d) \]

那么就是：

\[\large g(n)=\sum\limits_{n|d} \mu(\frac{d}{n}) \lfloor\dfrac Nd\rfloor\lfloor \dfrac Md\rfloor \]

然后我们换成直接枚举这个倍数 \(\large t=\frac dn\) ：

\[\large g(n)=\sum\limits_{t=1}^{\lfloor\frac{\min(N,M)}{n}\rfloor} \mu(t) \lfloor\dfrac N{tn}\rfloor\lfloor \dfrac M{tn}\rfloor \]

那么对于这个柿子可以应用整除分块来处理后半部分，然后预处理一下 \(\mu\) 的前缀和即可在 \(\sqrt{N}\) 的复杂度做到 \(O(n)\) 预处理的单次询问

注意这里的整除分块：

其实我们不是按照 \(tn\) 来分块（这样我就不会处理 \(\mu\) 里面的东西了），而是 \(t\) 来分块。

要把柿子转化成这样：

\[\large g(n)=\sum\limits_{t=1}^{\min(N,M)} \mu(t) \lfloor\dfrac{\lfloor\frac N{t}\rfloor}{n}\rfloor\lfloor\dfrac{\lfloor\frac M{t}\rfloor}{n}\rfloor \]

然后这个时候以 \(t\) 来分块看上去就很没有问题了，因为变量只有 \(t\) ，和 \(n\) 没有关系，我们只关注分子的变化即可。

接下来可以注意这样的分块，其实就等价于 \(N,M\) 都来分块，分了很多个隔板，然后块的个数也就是 \(2\sqrt{N}\) 级别了。

也就是每次取 \(r\) 的时候取两者算出来较小的右端点。

代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define ll long long
#define ull unsigned long long
#define inc(x,y,mod) (((x)+(y))>=(mod)?(x)+(y)-(mod):(x)+(y))
#define dec(x,y,mod) ((x)-(y)<0?(x)-(y)+(mod):(x)-(y))
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define dep(i,y,x) for(int i=(y);i>=(x);i--)
const int N=1e6+5,M=2e5+5,MOD=1e9+7;
int n,m,k,Ans;
int prime[N],mu[N],cnt,pre[N];
bool isprime[N];
void GetPrimes(int v){
	mu[1]=pre[1]=1;
	for(int i=2;i<=v;i++){
		if(!isprime[i]) prime[++cnt]=i,mu[i]=-1;
		pre[i]=pre[i-1]+mu[i];
		for(int j=1;j<=cnt&&prime[j]*i<=v;j++){
			isprime[i*prime[j]]=true;
			if(i%prime[j]==0) break;
			mu[i*prime[j]]=-mu[i];
		}
	}
	return ;
}
signed main(){
	ios::sync_with_stdio(false);
	GetPrimes(5e4);
	int t;cin>>t;
	while(t--){
		cin>>n>>m>>k;Ans=0;
		if(n<m) swap(n,m);
		for(int l=1,r;l<=m;l=r+1){
			r=min(n/(n/l),m/(m/l));
			Ans+=(pre[r]-pre[l-1])*(n/(l*k))*(m/(l*k));
		}
		cout<<Ans<<endl;
	}

	system("Pause");
	return 0;
}