[BZOJ 3309] DZY Loves Math

题意

求下式的值

[sum_{i=1}^nsum_{j=1}^mf(gcd(i,j)) ]

其中 (f(x)) 为 (x) 的质因子的最大幂次, (n,mle 1 imes 10^7, qle10000)

题解

首先按照以前反演的套路容易推出这个鬼式子:

[ ext{Ans}=sum_xleftlfloor frac n x ight floorleftlfloor frac m x ight floor sum_{d|x}mu(d)fleft(frac x d ight) ]

剩下的问题就是怎么搞出 (g(x)=sumlimits_{d|x}mu(d)fleft(frac x d ight)) 也就是 (g=mu*f) 的前缀和了

然而这个 (f) 不满足互质积性...不过它满足互质取 (max) !

(看完题目一开始以为不用卷上 (mu) 于是以为直接线筛就星了)

然而卷上一个积性函数之后就没互质类性质了...

我们重新来看这个函数: 因为系数上有个 (mu) , 所以只要 (d) 里面有 (p^2) 因子那整个项就没贡献了, 于是我们把这个和式的枚举条件转化为 (x) 的质因子的集合的子集枚举. 这个时候 (frac xd) 里面对应的质因子的幂次就会 (-1). 于是显然这个 (fleft(frac x d ight)) 只有两种取值: (f(x)) 和 (f(x)-1). 而且显然最高次的质因子和 (fleft(frac x d ight)) 的取值直接相关, 我们将 (x) 的质因子分为两个集合: (M) 代表最高次质因子, (R) 代表非最高次质因子. 然后我们考虑什么时候会取到这两种取值:

• 如果 (fleft(frac x d ight)=f(x)), 那么肯定最高次的质因子至少存在一个, 而不是最高次的质因子任意取. 这样的话共有 ((2^{|M|}-1) imes2^{|R|}) 个项. 此时作为系数的 (mu) 的取值根据选中的质因子个数的奇偶性也有(pm 1)两种各一半. 如果项数是偶数的话显然它就真的对半消掉了, 否则会留下一个 (pm f(x)) 没消. 此时 (|R|) 必定为 (0).

• 如果 (fleft(frac x d ight)=f(x)-1), 那么肯定是最高次的质因子都被选中了, 那么此时只有非最高次质因子还可以任意取, 于是项数是 (2^{|R|}), (mu) 的正负也一样是对半分的. 而只有当 (|R|=0) 的时候才可能是奇数, 此时剩下一个 (f(x)-1).

而对于整个集合 (P=Mcup R) 来说, (mu) 的取值必定是一半 (1) 和一半 (-1), 那么上面两种情况中没有消掉的 (mu) 一定会是一个 (1) 和一个 (-1), 于是 (g(x)) 只可能有三种取值: 如果所有质因子幂指数相同(都为最大), 则若有奇数个质因子则 (g(x)=1), 否则 (g(x)=-1), 否则 (g(x)=0).

然后就是欢脱的筛 (g) 的过程辣~

因为每个满足 (g(x) e0) 的 (x) 都满足是若干不同质数的乘积或这种乘积的某个整次幂, 于是我们只需要筛出不同质数的乘积然后一直算它的整次幂就好了...这个不知道什么筛的复杂度好像很正常甚至比线筛跑得要快...

代码实现

#include <bits/stdc++.h>

const int MAXN=1e7+10;

int cnt;
int g[MAXN];
int pr[MAXN];
bool npr[MAXN];
bool valid[MAXN];

void PowerSieve(int);

int main(){
	int T;
	scanf("%d",&T);
	PowerSieve(1e7);
	while(T--){
		int n,m;
		scanf("%d%d",&n,&m);
		if(n>m)
			std::swap(n,m);
		long long ans=0;
		for(int i=1,j;i<=n;i=j+1){
			j=std::min(n/(n/i),m/(m/i));
			ans+=1ll*(n/i)*(m/i)*(g[j]-g[i-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}

void PowerSieve(int n){
	npr[0]=npr[1]=true;
	for(int i=2;i<=n;i++){
		if(!npr[i]){
			pr[cnt++]=i;
			g[i]=1;
			valid[i]=true;
		}
		for(int j=0,t;j<cnt&&(t=i*pr[j])<=n;j++){
			npr[t]=true;
			if(i%pr[j]==0)
				break;
			else if(valid[i]){
				g[t]=-g[i];
				valid[t]=true;
			}
		}
		if(valid[i]){
			for(long long x=1ll*i*i;x<=n;x*=i)
				g[x]=g[i];
		}
	}
	for(int i=1;i<=n;i++){
		g[i]+=g[i-1];
	}
}