题意
设 \(f(n)\) 为 \(i\) 的第二大因数（\(f(n)=\dfrac{n}{\min(n)}\)），求 \(\sum_{i=1}^n \sum_{j=1}^n f(\gcd(i,j))^k\)。\(n\leq 10^9\)，\(k\leq 50\)。模 \(2^{32}\)。
题解
简单莫反可得：
\[ans=\sum_{i=2}^{n}f(i)^k(-1+2\sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}\varphi(i)) \]整除分块显然，\(\varphi\) 的前缀和可以杜教筛求出，问题在于 \(f(i)^k\) 的前缀和。
注意到 \(f(i)^k\) 和 \(i^k\) 相差的是 \(\min(i)^k\)，而 min_25 筛便将 \(\min(i)\) 不同的数分开来。观察 min_25 筛第一部分的式子（节选）：
\[g(n,i-1)-f'(p_i) \color{red}\left(g(\lfloor {n\over p_i}\rfloor,i-1)-\sum_{j=1}^{i-1}f'(p_j) \right) \]红色部分便是 \(\min(j)=p_i\) 的合数 \(j\) 的 \(f(j)^k\)。
质数要单独拿出来统计。
注意到模数十分阴间，预处理自然数幂和时最好用斯特林数+下降幂。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int M=7e4+10,L=2e6+10,K=57;
const ll N=1e9+10;
#define uint unsigned int
uint n,k,sqr,lim;
uint pri[L+10],sp1[M],sp2[M],boo[L+10],ppow[M],cnt=0;
uint id1[M],id2[M],g1[M],g2[M],g[M],w[M];
uint fs[M];
uint phi[L+10],phi2[L+10];
int _(ll x){ return x<sqr?id1[x]:id2[n/x]; }
uint qpow(uint x,uint y){
uint ans=1;
while(y){
if(y&1)ans=ans*x;
x=x*x;
y>>=1;
}
return ans;
}
uint s2[K][K];
void init(int k){
s2[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=i;j++)
s2[i][j]=s2[i-1][j]*j+s2[i-1][j-1];
}
uint s(uint x){
// sum_{i=0}^{x-1}
uint ans=0;
for(int i=0;i<=k;i++){
uint p=1;
for(int j=0;j<=i;j++){
if((x-j)%(i+1)==0)p*=(x-j)/(i+1);
else p*=x-j;
}
ans+=p*s2[k][i];
}
return ans;
}
uint sphi(uint x){
if(x<=lim)return phi[x];
if(phi2[n/x])return phi2[n/x];
uint ans=x*1ll*(x+1)/2;
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*sphi(x/l);
}
// cerr<<"! "<<x<<" "<<ans<<endl;
return phi2[n/x]=ans;
}
int main(){
cin>>n>>k;
init(k);
sqr=sqrt(n+1);
for(int i=2;i<=sqr;i++){
if(!boo[i]){
pri[++cnt]=i;
ppow[cnt]=qpow(i,k);
sp1[cnt]=(sp1[cnt-1]+ppow[cnt]);
sp2[cnt]=(sp2[cnt-1]+1);
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=sqr;j++){
boo[i*pri[j]]=1;
if(i%pri[j]==0)break;
}
}
int cnt_=cnt;
phi[1]=1;
lim=pow(n+1,0.68);
cnt=0;
for(int i=2;i<=lim;i++){
if(!boo[i]){
pri[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=lim;j++){
boo[i*pri[j]]=1;
if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
else{
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
}
}
for(int i=1;i<=lim;i++)phi[i]+=phi[i-1];
int m=1;
for(ll l=1,r;l<=n;l=r+1,m++){
r=n/(n/l);
uint t=n/l;
w[m]=t;
g1[m]=s(t+1)-1;
g2[m]=t-1;
if(t<sqr)id1[t]=m;
else id2[n/t]=m;
}
--m;
for(int i=1;i<=cnt_;i++){
int k=1;
for(int j=1;j<=m&&w[j]>=pri[i]*1ll*pri[i];j++){
ll r=_(w[j]/pri[i]);
fs[j]+=(g1[r]-sp1[i-1]);
g1[j]-=ppow[i]*(g1[r]-sp1[i-1]);
g2[j]-=(g2[r]-sp2[i-1]);
}
}
for(int i=1;i<=m;i++)fs[i]+=g2[i];
uint ans=0;
for(int i=1;i<=m;i++)
ans+=(fs[i]-fs[i+1])*(sphi(n/w[i])*2-1);
cout<<ans;
}