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