题目

有生之年我竟然能\(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;
}