Loj#572. 「LibreOJ Round #11」Misaka Network 与求和

题目

有生之年我竟然能\(A\)

这个题求的是这个

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

\(f(i)\)定义为\(i\)的次大质因子，其中\(f(p)=1,f(1)=0\)

看到这道题的第一反应肯定是这东西TM还能求

习惯性反演

\[\sum_{d=1}^nF(d)f(d)^k \]

\[=\sum_{d=1}^nf(d)^k\sum_{d|i}\mu(\frac{i}{d})\left \lfloor \frac{n}{i} \right \rfloor^2 \]

\[=\sum_{i=1}^n\left \lfloor \frac{n}{i} \right \rfloor^2\sum_{d|i}\mu(\frac{i}{d})f(d)^k \]

现在处理后面那个\(\mu\times f\)的前缀和就可以整除分块了

看到数据范围不小，考虑杜教筛，在卷上一个\(I\)，把\(\mu\)消掉

于是就有

\[S(n)=\sum_{i=1}^nf(i)^k-\sum_{i=2}^nI(i)S(\left \lfloor \frac{n}{i} \right \rfloor) \]

发现只需要求出\(\sum_{i=1}^nf(i)^k\)就可以杜教筛了

这个东西显然可以被Min_25处理，就是算次大质因子的时候贡献算成\(k\)次方即可

代码

#include <tr1/unordered_map>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 3000005
#define re register
#define LL unsigned long long
using namespace std::tr1;
unordered_map<LL,LL> ma;
LL n,p[maxn],w[maxn],id1[maxn],id2[maxn],K,mu[maxn],Sqr;
LL m,f[maxn],tot,F[maxn][2],pre[maxn],D[maxn],g[maxn];
inline LL quick(LL a,LL b) {LL s=1;while(b) {if(b&1ll) s*=a;b>>=1ll;a*=a;}return s;}
inline LL S(LL x,int y) {
	if(x<=1&&p[y]>x) return 0;
	LL ans=0;
	for(re int k=y;k<=tot&&1ll*p[k]*p[k]<=x;k++) {
		LL p1=p[k];
		for(;p1<=x;p1*=p[k]) {
			ans+=S(x/p1,k+1);
			LL t=g[(x/p1<=Sqr)?id1[x/p1]:id2[n/(x/p1)]];
			if(t+1>k) ans+=(t+1-k)*D[p[k]];
		}
	}
	return ans;
}
LL calc(LL x) {
	if(x<=Sqr) return pre[x];
	if(ma[x]) return ma[x];
	LL cnt=S(x,1)+g[id2[n/x]];
	for(re LL l=2,r;l<=x;l=r+1) {
		r=x/(x/l);cnt-=(r-l+1)*calc(x/l);
	}
	return ma[x]=cnt;
}
int main()
{
	scanf("%u%u",&n,&K);Sqr=std::pow(n,0.666);
	for(re LL l=1,r;l<=n;l=r+1) {
		r=n/(n/l);w[++m]=n/l;
		if(w[m]<=Sqr) id1[w[m]]=m;else id2[n/w[m]]=m;
		g[m]=w[m]-1;
	}
	f[1]=1;D[1]=1;mu[1]=1;
	for(re int i=2;i<=Sqr;i++) {
		if(!f[i]) p[++tot]=i,D[i]=quick(i,K),mu[i]=-1;
		for(re int j=1;j<=tot&&p[j]*i<=Sqr;j++) {
			f[p[j]*i]=1;if(i%p[j]==0) break;mu[p[j]*i]=-1*mu[i];
		}
	}
	for(re int j=1;j<=tot&&1ll*p[j]*p[j]<=n;j++) 
        for(re int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++) {
                int k=(w[i]/p[j]<=Sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
                g[i]=g[i]-g[k]+j-1;
            }
	memset(f,0,sizeof(f));f[1]=1;
	for(re int i=2;i<=Sqr;i++){
		if(f[i]) continue;
		F[i][1]=1;
		for(re int j=i+i;j<=Sqr;j+=i) f[j]=1,F[j][1]=F[j][0],F[j][0]=i;
	}
	for(re int i=2;i<=Sqr;i++) 
		if(f[i]&&!F[i][1]) F[i][1]=F[i][0];
	for(re int i=2;i<=Sqr;i++) 
		if(f[i]&&(i/F[i][0])%F[i][0]==0) F[i][1]=F[i][0]; 
	for(re int i=1;i<=Sqr;i++)
		for(re int j=i;j<=Sqr;j+=i) 
			pre[j]+=mu[j/i]*D[F[i][1]];
	for(re int i=1;i<=Sqr;i++) pre[i]+=pre[i-1];
	LL now=0;
	for(re LL l=1,r;l<=n;l=r+1) {
		r=n/(n/l);
		now+=(n/l)*(n/l)*(calc(r)-calc(l-1));
	}
	printf("%u",(unsigned int)now);
	return 0;
}