传送门:
https://www.luogu.com.cn/problem/P3327
https://www.acwing.com/problem/content/1360/
莫比乌斯反演 + 整除分块
分析
首先,我们给出一个结论:
证明:
设 (i) 的分解式为 (i = prod p_k^{alpha_k}) ,类似地,(j = prod p_k^{eta_k})
对于质数 (p_k) ,(ij) 对应的因数个数为 (alpha_k + eta_k + 1)
而对于右式,(p_k) 如果想要产生贡献,那么只可能是:(x, y) 不同时存在因数 (p_k),注意到对应的贡献数即为 (alpha_k + eta_k + 1)
因此,由乘法原理,所求证式成立。
为了防止混淆,我们将题面中的 (n, m) 记为 (N, M)
由上述结论,题目的式子转化为:
莫比乌斯反演:若 (F(n) = sum_{n|d}f(d)),则有 (f(n) = sum_{n|d} mu(frac{d}{n})F(d))
于是我们设 (f(n) = sum_{i=1}^N sum_{j=1}^M sum_{x|i} sum_{y|j} [(x, y) = n]) ,对应的 (F(n) = sum_{i=1}^N sum_{j=1}^M sum_{x|i} sum_{y|j} [n | (x, y)])
下面对 (F(n)) 进行变换:
(F(n) \= sum_{i=1}^N sum_{j=1}^M sum_{x|i} sum_{y|j} [n | (x, y)] \ = sum_{x=1}^N sum_{y=1}^M lfloor frac{N}{x} floor lfloor frac{M}{y} floor [n | (x, y)] \)
设 (x' = lfloor frac{x}{n} floor ~, y'=lfloor frac{y}{n} floor ~ , N'=lfloor frac{N}{n} floor ~ , M'=lfloor frac{M}{n} floor) ,进而有
(F(n) \= sum_{x=1}^{N'} sum_{y=1}^{M'} lfloor frac{N'}{x'} floor lfloor frac{M'}{y'} floor \ = sum_{x=1}^{N'} lfloor frac{N'}{x'} floor sum_{y=1}^{M'} lfloor frac{M'}{y'} floor)
记 (h(x) = sum_{i=1}^x lfloor frac{x}{i} floor) ,(h(x)) 可用整除分块处理。
由莫比乌斯反演,(f(n) = sum_{n|d} mu(frac{d}{n})F(d)) ,我们只需求 (f(1))
下面对 (f(1)) 进行变换:
由整除分块知,(h(lfloor frac{N}{d} floor)h(lfloor frac{M}{d} floor)) 取值最多为 (sqrt{N} + sqrt{M}) 块,而对于每一块,可以用前缀和 (O(1)) 处理出对应的 (mu) 值,因此整体的复杂度为 (O(T(sqrt{N} + sqrt{M})))
细节见代码:
#include<bits/stdc++.h>
using namespace std;
const int S=50050;
int primes[S], cnt;
bool vis[S];
int mu[S], sum[S];
int h[S];
int g(int b, int l){
return b/(b/l);
}
void init(){
// init sum[]
mu[1]=1;
for(int i=2; i<S; i++){
if(!vis[i]) primes[cnt++]=i, mu[i]=-1;
for(int j=0; i*primes[j]<S; j++){
vis[i*primes[j]]=true;
if(i%primes[j]==0) break;
mu[i*primes[j]]=-mu[i];
}
}
for(int i=1; i<S; i++) sum[i]=sum[i-1]+mu[i];
// init h[]
for(int x=1; x<S; x++){
int &v=h[x];
for(int l=1, r; l<=x; l=r+1){
r=min(x, g(x, l));
v+=x/l*(r-l+1);
}
}
}
int solve(int N, int M){
int res=0;
int n=min(N, M);
for(int l=1, r; l<=n; l=r+1){
r=min(n, min(g(N, l), g(M, l)));
res+=(sum[r]-sum[l-1])*h[N/l]*h[M/l];
}
return res;
}
int main(){
int T; cin>>T;
init();
while(T--){
int N, M; cin>>N>>M;
cout<<solve(N, M)<<endl;
}
return 0;
}