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

求 $sum_{i=1}^{n}sum_{j=1}^{m}lcm(i,j)mu(gcd(i,j))^2$

 

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

 

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

 

$Rightarrowsum_{d=1}^{n}mu(d)^2dsum_{i=1}^{left lfloor frac{n}{d} ight floor}sum_{j=1}^{left lfloor frac{m}{d} ight floor}ij[gcd(i,j)==1]$

 

省略一些反演的套路.....

 

$Rightarrow sum_{d=1}^{n}mu(d)^2dsum_{b=1}^{left lfloor frac{n}{d} ight floor}mu(b)b^2sum_{i=1}^{left lfloor frac{n}{db} ight floor}isum_{j=1}^{left lfloor frac{m}{d} ight floor}j$

 

令 $Sum(n,m)=frac{n(n+1)}{2} imesfrac{m(m+1)}{2}$

 

原式 $= sum_{d=1}^{n}mu(d)^2dsum_{b=1}^{left lfloor frac{n}{d} ight floor}mu(b)b^2Sum(left lfloor frac{n}{db} ight floor,left lfloor frac{m}{db} ight floor)$ 

 

令 $D=db$

 

$Rightarrowsum_{D=1}^{n}Sum(left lfloor frac{n}{D} ight floor,left lfloor frac{m}{D} ight floor)sum_{d|D}mu(d)^2d imesmu(frac{D}{d}) imes(frac{D}{d})^2$

 

令 $h(D)=sum_{d|D}mu(d)^2d imesmu(frac{D}{d}) imes(frac{D}{d})^2$

 

上式 $=sum_{D=1}^{n}Sum(left lfloor frac{n}{D} ight floor,left lfloor frac{m}{D} ight floor) imes h(D)$

 

如果能知道 $h(D)$ 的前缀和，就可以用 $sqrt n$ 的时间复杂度完成每一个询问

 

求 $h(D)$ 有两种方式，第一种是暴力枚举因数，直接更新，时间复杂度是正确的.

另一种方式是线筛 $h$ 函数，因为我们发现这个是积性的.

这里讲一下线性筛的做法（题解里写的是第一种）

 

(1) $D$ 为质数，$h(D)=D-D^2$

 

(2) $i,p$ 互质，那么 $h(i imes p)=h(i) imes h(p)$

 

(3) $i\%p==0$

 

先列出 $h(i)=sum_{d|i}mu(d)^2d imesmu(frac{i}{d}) imes(frac{i}{d})^2$

 

令 $t=frac{i}{p}$

 

(1)若 $t$ 中也含有 $p$ 这个因子，$i imes p$ 中 $p$ 的次幂至少为 $3$，无论如何 $h(i)$ 都是 $0$

 

(2)若 $t$ 中没有 $p$ 这个因子，考虑 $sum_{d|t}mu(d)^2d imesmu(frac{i}{d}) imes(frac{i}{d})^2$ 这个式子 在乘一个 $p$ 会导致这个式子变成原来的 $-p^2$ 倍 .
再考虑 $i$ 中至少有一个 $p$ 对 $h(i imes p)$ 的贡献：我们发现无论如何 $d$ 中至少会有 $2$ 个 $p$ 的因子，所有贡献为 $0$ . 这个情况下 $h(i imes p)=-h(t) imes 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
",re); 
	}
	return 0; 
}