记(f(x)=)(x)的次大因数，那么(sgcd(i,j)=f(gcd(i,j)))。
下面来推式子：
[egin{aligned}
sum_{i=1}^nsum_{j=1}^nsgcd(i,j)^k&=sum_{i=1}^nsum_{j=1}^nf(gcd(i,j))^k&记d=gcd(i,j)\
&=sum_{d=1}^nf(d)^ksum_{frac i d=1}^{lfloor frac n d
floor}sum_{frac j d=1}^{lfloor frac n d
floor}[gcd(frac i d,frac j d)=1]&f(1)=0\
&=sum_{d=2}^nf(d)^ksum_{i=1}^{lfloor frac n d
floor}sum_{j=1}^{lfloor frac n d
floor}[gcd(i,j)=1]\
&=sum_{d=2}^nf(d)^k(2sum_{i=1}^{lfloor frac n d
floor}varphi(i)-1)
end{aligned}
]
最后一步的括号是用欧拉函数的定义直接替换出来的。
我们发现，可以按(lfloor frac n d
floor)的取值数论分段，因为括号显然只受(lfloor frac n d
floor)的取值影响，关键是如何求(f(x)^k)的前缀和？
记(S_x=sum_{i=2}^xf(x)^k)，
令min_25中的(g)数组以(s(x)=x^k)的定义计算(g)，考虑由(g_{i,j-1})推到(g_{i,j})的时候，减去的是什么？是(s(p_j)*(g_{lfloorfrac n {p_j}
floor,j-1}-g_{p_j-1,j-1}))，后面括号的部分是什么呢？恰好是那些最小质因子为(p_j)且除去(p_j)后剩余部分最小质因数不小于(p_j)的数的(k)次方之和，我们发现这些数的(f)值之和就是后面的括号。又因为每一个数至多被如此枚举到1次，所以对于每一个(i)，我们把括号的值累加起来，这就是那些([2,i])中非质数的(f)值之和；再加上([2,i])中的质数的(f)值之和，也就是质数个数，我们就求得了(S_i)。
min_25真牛逼。
所以就做完了，(g)的初始化需要用到自然数幂求和，考虑到(k)比较小，可以用第二类斯特林数求解。
[egin{aligned}
Sum_k(n)&=sum_{i=1}^ni^k\
&=sum_{j=1}^kegin{Bmatrix}k\jend{Bmatrix}frac{{(n+1)}^underline{j+1}}{j+1}
end{aligned}
]
Code
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <map>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef map<ll,ui> mlu;
const int SQRTN=32000,LEN=1000001;
bool vis[LEN];
ui p[LEN],pcnt,pphi[LEN],s[51][51],pk[LEN];
ll n,k,sqrtn,m;
ll a[SQRTN*2],cnt,pos1[SQRTN],pos2[SQRTN];
ui g[SQRTN*2],sum[SQRTN*2],g1[SQRTN*2];
mlu record;
ui ksm(ui x,int y){
ui res=1;
for(;y;x=x*x,y>>=1)
if(y&1) res=res*x;
return res;
}
void prework(){
pphi[1]=1;
for(int i=2;i<LEN;i++){
if(!vis[i]){
p[++pcnt]=i;
pphi[i]=i-1;
}
for(int j=1;j<=pcnt&&i*p[j]<LEN;j++){
int x=i*p[j];
vis[x]=true;
if(i%p[j]==0){
pphi[x]=pphi[i]*p[j];
break;
}
pphi[x]=pphi[i]*(p[j]-1);
}
}
for(int i=2;i<LEN;i++) pphi[i]+=pphi[i-1];
s[0][0]=s[0][1]=1;
for(int i=1;i<=50;i++){
s[i][1]=1;
for(int j=2;j<=i;j++)
s[i][j]=s[i-1][j]*(ui)j+s[i-1][j-1];
}
}
inline int gp(ll x){
return x<=sqrtn?pos1[x]:pos2[n/x];
}
void Diz(){
for(ll i=1,j;i<=n;i=j+1){
j=n/(n/i);
a[++cnt]=n/i;
}
reverse(a+1,a+1+cnt);
for(int i=1;i<=cnt;i++)
if(a[i]<=sqrtn) pos1[a[i]]=i;
else pos2[n/a[i]]=i;
}
ui sumup(ll l,ll r){
ll a=l+r,b=r-l+1;
if(a&1) b/=2; else a/=2;
return (ui)a*(ui)b;
}
ui getPhi(ll n){
if(n<LEN) return pphi[n];
mlu::iterator pos=record.find(n);
if(pos!=record.end()) return pos->second;
ui res=sumup(1,n);
for(ll i=2,j;i<=n;i=j+1){
j=n/(n/i);
res-=getPhi(n/i)*(j-i+1);
}
record[n]=res;
return res;
}
ui getSumk(ll n){
ui res=0,t;
for(ll j=1;j<=k;j++){
t=1;
int md=(n-j+1)%(j+1);
for(ll i=n-j+1;i<=n+1;i++,md++)
if(!md||md==(j+1)) t*=(ui)i/(j+1);
else t*=(ui)i;
res+=t*s[k][j];
}
return res;
}
void calc_g(){
for(int i=2;i<=cnt;i++) g[i]=getSumk(a[i])-1,g1[i]=a[i]-1,sum[i]=0;
for(int j=1;j<=m;j++)
for(int i=cnt;i>=2&&a[i]>=p[j]*p[j];i--){
g1[i]-=g1[gp(a[i]/p[j])]-g1[gp(p[j]-1)];
ui delta=(g[gp(a[i]/p[j])]-g[gp(p[j]-1)]);
g[i]-=pk[j]*delta;
sum[i]+=delta;
}
}
ui calc(ll n){
return sum[gp(n)]+g1[gp(n)];
}
int main
prework();
scanf("%lld%lld",&n,&k);
sqrtn=(ll)sqrt(n);
m=upper_bound(p+1,p+1+pcnt,sqrtn)-p-1;
for(int i=1;i<=m;i++) pk[i]=ksm(p[i],k);
Diz();
calc_g();
ui ans=0,last=0,tmp;
for(ll i=2,j;i<=n;i=j+1){
j=n/(n/i);
tmp=calc(j);
ans+=(getPhi(n/i)*2-1)*(tmp-last);
last=tmp;
}
printf("%lu
",ans);
return 0;
}