[BZOJ 3309] DZY Loves Math

题意

求下式的值

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

其中 \(f(x)\) 为 \(x\) 的质因子的最大幂次, \(n,m\le 1\times 10^7, q\le10000\)

题解

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

\[\text{Ans}=\sum_x\left\lfloor \frac n x \right \rfloor\left\lfloor \frac m x \right \rfloor \sum_{d|x}\mu(d)f\left(\frac x d \right) \]

剩下的问题就是怎么搞出 \(g(x)=\sum\limits_{d|x}\mu(d)f\left(\frac x d\right)\) 也就是 \(g=\mu*f\) 的前缀和了

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

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

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

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

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

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

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

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

因为每个满足 \(g(x)\ne0\) 的 \(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\n",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];
	}
}