题面
题解
这题有毒……不知为啥的错误调了半天……
令(f(i)={sgcd(i)})，那么容易看出(f(i))就是(i)的次大质因子，用(i)除以它的最小质因子即可计算
于是题目所求即为
[sum_{i=1}^nsum_{j=1}^n{f(gcd(i,j))}^k
]
[sum_{d=1}^n {f(d)}^ksum_{i=1}^{leftlfloorfrac{n}{d}
ight
floor}sum_{j=1}^{leftlfloorfrac{n}{d}
ight
floor}[gcd(i,j)=1]
]
[sum_{d=1}^n {f(d)}^k(2sum_{i=1}^{leftlfloorfrac{n}{d}
ight
floor}varphi(i)-1)
]
后面那个根据(varphi)的定义就可以推出来
因为后面括号中的式子只与({leftlfloorfrac{n}{d} ight floor})的值有关，我们可以考虑整除分块，那么括号里的东西就是一个杜教筛
那么我们现在的问题就是对于前半部分怎么快速求和了。我们考虑(Min\_25)筛的过程，其中(g(n,j))表示所有(1)到(n)中的质数或最小质因子大于(P_j)的所有数的(k)次方之和，即
[g(n,j)=sum_{i=1}^ni^k[iin P or min_p(i)>P_j]
]
中间的转移中有这么一步
[g(n,i)=g(n,j-1)-{P_j}^k(g({leftlfloorfrac{n}{P_j}
ight
floor},j-1)-sum_{i=1}^{j-1}{P_i}^k)
]
考虑后面减去的东西，是所有最小质因子等于(P_j)且自身为合数的数的(k)次方之和，即
[sum_xx^k[x
otin P and min_p(x)=P_j]
]
然后我们惊喜地发现，后面减去的东西中，括号里的东西就是上面所有满足条件的(x)的({f(x)}^k)之和！
那么我们就可以直接(Min\_25)筛来计算({f(x)}^k)的前缀和了，那么一段区间只要差分一下就可以了
最后有个尴尬的问题，就是这个模数不是质数，所以我们在初始化的时候需要快速计算(sum_{i=1}^n i^k)就不能拉格朗日插值了
那么只好用第二类斯特林数优化了
[egin{aligned}
sumlimits_{i=1}^{n}i^K&=sum_{i=1}^nsum_{j=1}^Kegin{Bmatrix}K\jend{Bmatrix}i^underline{j}\
&=sum_{j=1}^Kegin{Bmatrix}K\jend{Bmatrix}sum_{i=1}^ni^underline{j}\
&=sum_{j=1}^Kegin{Bmatrix}K\jend{Bmatrix}frac{{(n+1)}^underline{j+1}}{j+1}
end{aligned}]
然后没有然后了
//minamoto
#include<bits/stdc++.h>
#define R register
#define int unsigned int
#define IT map<int,int>::iterator
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=1e5+5;
int w[N<<1],g[N<<1],h[N<<1],G[N<<1],id1[N],id2[N],p[N],pw[N],sp[N],phi[N],S[55][55];
map<int,int>mp;bitset<N>vis;int n,m,k,tot,sqr,id,ans,now,las;
int ksm(R int x,R int y){
R int res=1;
for(;y;y>>=1,x*=x)if(y&1)res*=x;
return res;
}
void init(int n){
phi[1]=1;
fp(i,2,n){
if(!vis[i])p[++tot]=i,pw[tot]=ksm(i,k),sp[tot]=sp[tot-1]+pw[tot],phi[i]=i-1;
for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){
vis[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
fp(i,1,n)phi[i]+=phi[i-1];
S[0][0]=1;
fp(i,1,k){
S[i][1]=1;
fp(j,2,i)S[i][j]=S[i-1][j]*j+S[i-1][j-1];
}
}
int Phi(int n){
if(n<=sqr)return phi[n];
IT it=mp.find(n);
if(it!=mp.end())return it->second;
int res=n*(n+1)>>1;
for(R int i=2,j;i<=n;i=j+1)
j=n/(n/i),res-=Phi(n/i)*(j-i+1);
return mp[n]=res;
}
int calc(int n){
int res=0,tmp;
fp(i,1,min(k,n)){
tmp=S[k][i];
fp(j,n-i+1,n+1)(j%(i+1)==0)?tmp*=j/(i+1):tmp*=j;
res+=tmp;
}return res;
}
signed main(){
// freopen("testdata.in","r",stdin);
scanf("%u%u",&n,&k),init(sqr=sqrt(n));
for(R int i=1,j;i<=n;i=j+1){
j=n/(n/i),w[++m]=n/i;
w[m]<=sqr?id1[w[m]]=m:id2[n/w[m]]=m;
g[m]=calc(w[m])-1,h[m]=w[m]-1;
}
fp(j,1,tot)for(R int i=1;1ll*p[j]*p[j]<=w[i];++i){
id=(w[i]/p[j]<=sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
g[i]-=pw[j]*(g[id]-sp[j-1]);
h[i]-=h[id]-j+1;
G[i]+=g[id]-sp[j-1];
}
for(R int i=2,j;i<=n;i=j+1){
j=n/(n/i),id=(j<=sqr)?id1[j]:id2[n/j];
now=G[id]+h[id],ans+=((Phi(n/i)<<1)-1)*(now-las),las=now;
}
printf("%u
",ans);
return 0;
}