BZOJ 2694: Lcm 莫比乌斯反演 + 积性函数 + 线性筛 + 卡常

时间:2021-08-22 16:47:40
求 $\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)\mu(gcd(i,j))^2$
 
$\Rightarrow \sum_{d=1}^{n}\mu(d)^2\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{d}[gcd(i,j)==d]$
 
$\Rightarrow \sum_{d=1}^{n}\mu(d)^2\frac{1}{d}\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)==d]$
 
$\Rightarrow\sum_{d=1}^{n}\mu(d)^2d\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}ij[gcd(i,j)==1]$
 
省略一些反演的套路.....
 
$\Rightarrow \sum_{d=1}^{n}\mu(d)^2d\sum_{b=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(b)b^2\sum_{i=1}^{\left \lfloor \frac{n}{db} \right \rfloor}i\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}j$
 
令 $Sum(n,m)=\frac{n(n+1)}{2}\times\frac{m(m+1)}{2}$
 
原式 $= \sum_{d=1}^{n}\mu(d)^2d\sum_{b=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(b)b^2Sum(\left \lfloor \frac{n}{db} \right \rfloor,\left \lfloor \frac{m}{db} \right \rfloor)$ 
 
令 $D=db$
 
$\Rightarrow\sum_{D=1}^{n}Sum(\left \lfloor \frac{n}{D} \right \rfloor,\left \lfloor \frac{m}{D} \right \rfloor)\sum_{d|D}\mu(d)^2d\times\mu(\frac{D}{d})\times(\frac{D}{d})^2$
 
令 $h(D)=\sum_{d|D}\mu(d)^2d\times\mu(\frac{D}{d})\times(\frac{D}{d})^2$
 
上式 $=\sum_{D=1}^{n}Sum(\left \lfloor \frac{n}{D} \right \rfloor,\left \lfloor \frac{m}{D} \right \rfloor)\times h(D)$
 
如果能知道 $h(D)$ 的前缀和,就可以用 $\sqrt n$ 的时间复杂度完成每一个询问
 
求 $h(D)$ 有两种方式,第一种是暴力枚举因数,直接更新,时间复杂度是正确的.
另一种方式是线筛 $h$ 函数,因为我们发现这个是积性的.
这里讲一下线性筛的做法(题解里写的是第一种)
 
(1) $D$ 为质数,$h(D)=D-D^2$
 
(2) $i,p$ 互质,那么 $h(i\times p)=h(i)\times h(p)$
 
(3) $i\%p==0$
 
先列出 $h(i)=\sum_{d|i}\mu(d)^2d\times\mu(\frac{i}{d})\times(\frac{i}{d})^2$
 
令 $t=\frac{i}{p}$
 
(1)若 $t$ 中也含有 $p$ 这个因子,$i\times p$ 中 $p$ 的次幂至少为 $3$,无论如何 $h(i)$ 都是 $0$
 
(2)若 $t$ 中没有 $p$ 这个因子,考虑 $\sum_{d|t}\mu(d)^2d\times\mu(\frac{i}{d})\times(\frac{i}{d})^2$ 这个式子 在乘一个 $p$ 会导致这个式子变成原来的 $-p^2$ 倍 .
再考虑 $i$ 中至少有一个 $p$ 对 $h(i\times p)$ 的贡献:我们发现无论如何 $d$ 中至少会有 $2$ 个 $p$ 的因子,所有贡献为 $0$ . 这个情况下 $h(i\times p)=-h(t)\times p^3$ 
#include<bits/stdc++.h>
#define maxn 4020002
#define ll long long
#define M 4000001
using namespace std;
#define O2 __attribute__((optimize("-O2")))
char *p1,*p2,buf[100000];
#define nc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
int rd() {int x=0; char c=nc(); while(c<48) c=nc(); while(c>47) x=(((x<<2)+x)<<1)+(c^48),c=nc(); return x;}
void setIO(string s)
{
string in=s+".in";
string out=s+".out";
freopen(in.c_str(),"r",stdin);
freopen(out.c_str(),"w",stdout);
}
int cnt;
bool vis[maxn];
int prime[maxn], mu[maxn];
ll g[maxn],mod=1,sum[maxn];
O2 inline void Init()
{
int i,j;
mu[1]=1;
for(i=1;i<=30;++i) mod=mod*2;
for(i=2;i<=M;++i)
{
if(!vis[i]) prime[++cnt]=i, mu[i]=-1;
for(j=1;j<=cnt&&1ll*prime[j]*i<=M;++j)
{
vis[prime[j]*i]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(i=1;i<=M;++i)
for(j=1;1ll*i*j<=M;++j)
if(mu[j]) g[i*j]+=mu[i]*i, g[i*j]%mod;
for(i=1;i<=M;++i) g[i]=(1ll*g[i]*i+g[i-1]+mod)%mod;
for(i=1;i<=M;++i) sum[i]=(1ll*(i+1)*i/2)%mod;
}
O2 int main()
{
// setIO("input");
Init();
int T,n,m;
T=rd();
while(T--)
{
n=rd(), m=rd();
if(n>m) swap(n,m);
int i,j;
ll re=0;
for(i=1;i<=n;i=j+1)
{
j=min(n/(n/i), m/(m/i));
re=(re + (1ll*sum[n/i]*sum[m/i]*(g[j]-g[i-1])%mod)%mod + mod)%mod;
}
printf("%lld\n",re);
}
return 0;
}