LOJ572 Misaka Network 与求和

Misaka Network 与求和

一方通行成功接入了 Misaka Network。

现在他要使用超能力，自然计算式被送到了御坂网络进行处理。这次的计算式是这样子的：

\[\sum_{i=1}^{N}\sum_{j=1}^{N}f(\gcd(i,j))^k \mod 2^{32} \]

其中 \(f(x)\) 表示 \(x\) 次大的质因数，重复的质因数计算多次，例如 \(f(6)=2,f(4)=2\)。规定 \(f(1)=0,f(p)=1\)，其中 \(p\) 为质数。

但是妹妹们都不会算这个式子……所以御坂 20001 号找到了你，希望你帮她算一下。

对于所有数据 \(1 \leq N,k \leq 2 \times 10^9\)。

题解

http://jklover.hs-blog.cf/2020/04/19/Loj-572-Misaka-Network-与求和/#more

为了方便，记\(f(i)\)就表示次大质因子的\(k\)次幂。

\[\begin{aligned} \text{ans}&=\sum_{i=1}^n\sum_{j=1}^n f(\gcd(i,j)) \\ &=\sum_{d=1}^n f(d) \sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{^{\lfloor\frac n d\rfloor}} [\gcd(i,j)=1] \\ &=\sum_{d=1}^nf(d)\sum_{t=1}^{\lfloor\frac n d\rfloor}\mu(t)\lfloor\frac n{dt}\rfloor^2\\ &=\sum_{d=1}^n \lfloor\frac{n}{d}\rfloor^2\sum_{t|d}\mu(t)f(\frac{d}{t})\\ &=\sum_{d=1}^n \lfloor\frac{n}{d}\rfloor^2\cdot (f*\mu)(d) \end{aligned} \]

对\(d\)整除分块，就需要多次询问\(f* \mu\)的前缀和。记\(F(n)=\sum_{i=1}^n (f* \mu)(i)\)。

用杜教筛的套路，可以得出

\[\sum_{i=1}^n f(n) = \sum_{i=1}^n (f * \mu * I) (i)\\ = \sum_{i=1}^n (F * I) (i) = \sum_{i=1}^n F(\lfloor \frac n i \rfloor)\\ \]

\[F(n)=\sum_{i=1}^n f(i)-\sum_{i=2}^n F(\lfloor \frac n i \rfloor) \]

计算\(F\)时整除分块，递归下去计算。

每次还会询问\(f\)的前缀和，可以把Min_25筛魔改一下，变成

• \(G(n,m)\)表示\(\leq n\)的合数，次大质因子\(\geq p_m\)的\(f\)的和。

• \(H(n)\)表示\(\leq n\)的质数的\(f\)的和。

这样就比较好处理了。

因为我偷懒没有预处理杜教筛，所以时间复杂度\(O(n^{\frac 34})\)，相当的慢。

IN uint fpow(uint a,int b){
	uint ans=1;
	for(;b;b>>=1,a=a*a)
		if(b&1) ans=ans*a;
	return ans;
}

CO int N=2e5+10;
int n,m,prime[N],num; uint kth[N];
int pos[N],idx,ref1[N],ref2[N]; uint H[N];

uint G(int x,int j){
	if(x<=1 or prime[j]>x) return 0;
	uint ans=0;
	for(int i=j;i<=num and prime[i]*prime[i]<=x;++i)
		for(int64 pw=prime[i];pw*prime[i]<=x;pw*=prime[i]){
			int k=x/pw;
			k=k<=m?ref1[k]:ref2[n/k];
			ans+=G(x/pw,i+1)+kth[prime[i]]*(H[k]-i+1);
		}
	return ans;
}

uint _F[N];

uint F(int x){
	if(x==0) return 0;
	int k=x<=m?ref1[x]:ref2[n/x];
	if(_F[k]) return _F[k];
	uint ans=G(x,1)+H[k];
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=F(x/l)*(r-l+1);
	}
	return _F[k]=ans;
}

int main(){
	read(n),m=ceil(sqrt(n));
	int K=read<int>();
	kth[1]=1;
	for(int i=2;i<=m;++i){
		if(!prime[i]) prime[++num]=i,kth[i]=fpow(i,K);
		for(int j=1;j<=num and i*prime[j]<=m;++j){
			prime[i*prime[j]]=1;
			if(i%prime[j]==0) break;
		}
	}
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		pos[++idx]=n/l;
		n/l<=m?ref1[n/l]=idx:ref2[n/(n/l)]=idx;
	}
	for(int i=1;i<=idx;++i) H[i]=pos[i]-1;
	for(int j=1;j<=num;++j)
		for(int i=1;i<=idx and prime[j]*prime[j]<=pos[i];++i){
			int k=pos[i]/prime[j];
			k=k<=m?ref1[k]:ref2[n/k];
			H[i]-=H[k]-(j-1);
		}
	uint ans=0;
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans+=(F(r)-F(l-1))*(n/l)*(n/l);
	}
	printf("%u\n",ans);
	return 0;
}