求
[sum_{i=1}^nsum_{j=1}^nsgcd(i,j)^k
]
首先这个次大公约数显然就是(gcd)除一下最小质因子了
于是
[sum_{i=1}^nsum_{j=1}^n(frac{(i,j)}{minp((i,j))})^k
]
显然可以枚举(gcd)，之后利用欧拉函数来求和
[sum_{d=1}^n(frac{d}{minp(d)})^k(2sum_{i=1}^{frac{n}{d}}varphi(i)-1)
]
后面那个东西显然可以直接杜教筛了，但是我们想要整除分块的话就需要前面这个函数的前缀和
考虑到最小质因子是一个跟埃氏筛非常有关系的东西，所以我们得知这题的正解应该是min_25
考虑我们用min_25去筛一个函数叫(f(i)=i^k)，我们只需要在min_25的过程中把每一个质数作为最小质因子的那些数的(k)次方和筛出来就好了
显然我们这样并没有计算(d)为质数的情况，于是我们顺便筛一个质数个数就好了
现在的问题变成了如何求自然数得幂次方和这样一个经典的问题，显然我们随便插个值就好了
但是就会发现模数不是质数，需要求逆元的方法就都凉了
有这样一个非常牛逼的做法是不需要逆元的，就是斯特林数
[sum_{i=1}^ni^k=sum_{i=1}^kS_2(k,i)frac{(n+1)^{underline {(j+1)}}}{j+1}
]
不想推了，就抄yyb的吧
至于上面那个柿子为明明有分母却不需要逆元是因为(n+1)
的(j+1)次下降幂一定会有一个数是(j+1)的倍数，只要我们乘这个数得时候乘上这个数除以(j+1)的商就好了
代码
#include<cmath>
#include<cstdio>
#define re register
#define uint unsigned int
const int maxn=1e5+5;
int n,k,U,Sqr,o,m;
uint phi[1500005],g[maxn],S[2][55],s[maxn],pre[maxn],sum[maxn],h[maxn];
int f[1500005],p[500005],w[maxn],id1[maxn],id2[maxn];
inline int chk(int x) {return x<=Sqr?id1[x]:id2[n/x];}
inline uint ksm(uint a,int b) {
uint S=1;
while(b) {if(b&1) S*=a;b>>=1;a*=a;}
return S;
}
uint Sphi(int x) {
if(x<=U) return phi[x];
if(sum[n/x]) return sum[n/x];
uint now=1ll*x*(x+1)/2ll;
for(re int l=2,r;l<=x;l=r+1) {
r=x/(x/l);
now-=Sphi(x/l)*(r-l+1);
}
return sum[n/x]=now;
}
inline uint calc(int n) {
uint cnt=0;
for(re int i=1;i<=k;i++) {
uint ans=1;
int flag=0;
for(re int t=i+1,j=n+1;t;--j,--t)
if(!flag&&j%(i+1)==0) flag=1,ans*=(j/(i+1));else ans*=j;
cnt+=S[o][i]*ans;
}
return cnt;
}
int main() {
scanf("%d%d",&n,&k);U=std::pow(n,0.666)+1;Sqr=std::sqrt(n)+1;
f[1]=1,phi[1]=1;o=0;S[0][1]=1;
for(re int i=2;i<=U;i++) {
if(!f[i]) p[++p[0]]=i,phi[i]=i-1,pre[p[0]]=pre[p[0]-1]+ksm(i,k);
for(re int j=1;j<=p[0]&&p[j]*i<=U;j++) {
f[p[j]*i]=1;
if(i%p[j]==0) {phi[p[j]*i]=p[j]*phi[i];break;}
phi[p[j]*i]=(p[j]-1)*phi[i];
}
}
for(re int i=2;i<=k;i++,o^=1)
for(re int j=1;j<=i;j++)
S[o^1][j]=S[o][j-1]+S[o][j]*(uint)j;
for(re int i=1;i<=U;i++) phi[i]+=phi[i-1];
for(re int 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]=calc(w[m])-1;h[m]=w[m]-1;
}
for(re int j=1;j<=p[0]&&1ll*p[j]*p[j]<=n;j++) {
uint tmp=ksm(p[j],k);
for(re int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++) {
int t=chk(w[i]/p[j]);
uint tot=tmp*(g[t]-pre[j-1]);
s[i]+=g[t]-pre[j-1];g[i]-=tot;
h[i]-=(h[t]-j+1);
}
}
for(re int i=1;i<=m;i++) s[i]+=h[i];
uint ans=0;
for(re int l=1,r;l<=n;l=r+1) {
r=n/(n/l);
ans+=(2ll*Sphi(n/l)-1)*(s[chk(r)]-s[chk(l-1)]);
}
printf("%u
",ans);
return 0;
}