题面
https://www.luogu.com.cn/problem/P3327
前置知识
- 线性筛积性函数:https://www.cnblogs.com/zhoushuyu/p/8275530.html
- 莫比乌斯反演:《具体数学:计算机科学基础》4.9节
- 莫比乌斯反演的基础应用: https://www.cnblogs.com/xh092113/p/12268356.html
- 数论分块(接下来会介绍)
数论分块
用来在(O(sqrt{n}))的时间内解决对({sum_{i=1}^{n}}f({lfloor}{frac{n}{i}}{ floor}))求和或类似问题的一种技巧。
性质:({lfloor}{frac{n}{i}}{ floor}(1{leq}i{leq}n)) 的取值最多有(O(sqrt{n}))个。
证:(1{leq}i{leq}{sqrt{n}})时,取值最多为i的个数({sqrt{n}})。
({sqrt{n}}{leq}i{leq}n)时,有(1{leq}{lfloor}{frac{n}{i}}{ floor}{leq}{sqrt{n}}),又({lfloor}{frac{n}{i}}{ floor}{in}Z),所以最多也是({sqrt{n}})。
性质得证。
故如果将1~n这n个数,({lfloor}{frac{n}{i}}{ floor})取值相同的合并为一段,最终我们会得到(O({sqrt{n}}))段。
而如果我们知道了某一段的左端点L,就可以通过
的方式,得到这一段的右端点R。之所以这个式子成立,是因为R是使得(x*{{lfloor}{frac{n}{L}}{ floor}}{leq}n)的最大的x。
之后,再使L=R+1,就可以得到下一段的左端点。直到R=n终止。
这样,就可以用O(1)的时间对每一段求和,所以总的时间复杂度也就是段数(O({sqrt{n}}))。
题解
引理:(d(ij) = {sum_{p|i}}{sum_{q|j}}[(p,q)=1])。
证:设(i={prod_{i=1}^{k}}p_k^{alpha_k},j={prod_{i=1}^{k}}p_k^{eta_k})。其中(p_1{leq}p_2{leq}…{leq}p_k)为质数,(max({alpha_i},{eta_i}){geq}1)。
由公式可得(d(ij)={prod_{i=1}^{k}}({alpha_i}+{eta_i}+1))。(1)
另一方面,
对于(0{leq}{gamma}{leq}{alpha},0{leq}{ heta}{leq}{eta}),要满足(min({gamma},{ heta})=0)只能有((0,1)、(0,2)…(0,{eta})),((1,0)、(2,0)…(alpha,0))以及((0,0))这({alpha}+{eta}+1)组解。
故上式(={prodlimits_{i=1}^{k}}(alpha_i+eta_i+1))=(1)式。引理得证。
回原题
由引理,
为化解gcd,我们使用技巧,此处(f(x)=[x=1]=e(x),g(x)=f{ imes}{mu}={mu})。
则原式=
进而,我们观察后两个括号内的数,其实就是(h({lfloor}{frac{n}{d}}{ floor}))和(h({lfloor}{frac{m}{d}}{ floor})),其中h(n)是为1~n的约数个数和,也就是({sum_{i=1}^{n}}d(i))。(这里的d不是前面推式子的时候的变量,而是题目中就已定义的那个函数)
过程中,之所以({sum_{i=1}^{n}}d(i)={sum_{i=1}^{n}}{lfloor}{frac{n}{i}}{ floor}),是因为
我们可以用埃筛预处理出所有的d(i),然后求前缀和得到h(i)。
目前的时间复杂度是O(Tn),不能满足要求。
注意到h的形式满足数论分块的要求,故此处可以用数论分块进行优化。
但是需要注意的是,这里要对数论分块进行一个小拓展,即({lfloor}{frac{n}{d}}{ floor})与({lfloor}{frac{m}{d}}{ floor})变化的时候都需要分段,故这里应该分成最多(O({sqrt{n}}+{sqrt{m}}))段,求右端点的时候也应该从本来的
R = n / (n / L);
拓展为
R = min(n / (n/L),m / (m/L));
分块结束后,同一段的一起处理,只需要先行预处理({mu(i)})的前缀和即可O(1)统计一段。
至此总时间复杂度为(O(T{sqrt{n}}))。
代码
#include<bits/stdc++.h>
using namespace std;
#define N 50000
#define rg register
#define ll long long
inline ll read(){
ll s = 0,ww = 1;
char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
return s * ww;
}
inline void write(ll x){
if(x < 0)putchar('-'),x = -x;
if(x > 9)write(x / 10);
putchar('0' + x % 10);
}
ll pn;
ll pri[10000+5],mu[N+5],sm[N+5]; //sm是mu的前缀和
bool isp[N+5];
inline void Eular(){
pn = 0;
mu[1] = 1;
for(rg ll i = 2;i <= N;i++)isp[i] = 1;
for(rg ll i = 2;i <= N;i++){
if(isp[i])pri[++pn] = i,mu[i] = -1;
for(rg ll j = 1;i * pri[j] <= N;j++){
isp[i*pri[j]] = 0;
if(i % pri[j])mu[i*pri[j]] = -mu[i];
else{
mu[i*pri[j]] = 0;
break;
}
}
}
}
ll d[N+5],h[N+5];
inline void prepro(){ //预处理d[]和h[](h即d的前缀和)
for(rg ll i = 1;i <= N;i++){
for(rg ll j = i;j <= N;j += i)d[j]++;
}
for(rg ll i = 1;i <= N;i++)h[i] = h[i-1] + d[i];
}
int main(){
Eular();
for(rg ll i = 1;i <= N;i++)sm[i] = sm[i-1] + mu[i];
prepro();
ll T = read();
while(T--){
ll n = read(),m = read();
ll limit = min(n,m);
ll ans = 0;
ll L = 1,R = min(n / (n/L),m / (m/L));
while(1){
ans += (sm[R] - sm[L-1]) * h[n/L] * h[m/L];
if(R == limit)break;
L = R + 1,R = min(n / (n/L),m / (m/L));
}
write(ans);putchar('
');
}
return 0;
}