题面

题解

这题有毒……不知为啥的错误调了半天……

令\(f(i)={sgcd(i)}\)，那么容易看出\(f(i)\)就是\(i\)的次大质因子，用\(i\)除以它的最小质因子即可计算

于是题目所求即为

\[\sum_{i=1}^n\sum_{j=1}^n{f(\gcd(i,j))}^k\]

\[\sum_{d=1}^n {f(d)}^k\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]\]

\[\sum_{d=1}^n {f(d)}^k(2\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\varphi(i)-1)\]

后面那个根据\(\varphi\)的定义就可以推出来

因为后面括号中的式子只与\({\left\lfloor\frac{n}{d}\right\rfloor}\)的值有关，我们可以考虑整除分块，那么括号里的东西就是一个杜教筛

那么我们现在的问题就是对于前半部分怎么快速求和了。我们考虑\(Min\_25\)筛的过程，其中\(g(n,j)\)表示所有\(1\)到\(n\)中的质数或最小质因子大于\(P_j\)的所有数的\(k\)次方之和，即

\[g(n,j)=\sum_{i=1}^ni^k[i\in P\ or\ min_p(i)>P_j]\]

中间的转移中有这么一步

\[g(n,i)=g(n,j-1)-{P_j}^k(g({\left\lfloor\frac{n}{P_j}\right\rfloor},j-1)-\sum_{i=1}^{j-1}{P_i}^k)\]

考虑后面减去的东西，是所有最小质因子等于\(P_j\)且自身为合数的数的\(k\)次方之和，即

\[\sum_xx^k[x\notin P\ and\ min_p(x)=P_j]\]

然后我们惊喜地发现，后面减去的东西中，括号里的东西就是上面所有满足条件的\(x\)的\({f(x)}^k\)之和！

那么我们就可以直接\(Min\_25\)筛来计算\({f(x)}^k\)的前缀和了，那么一段区间只要差分一下就可以了

最后有个尴尬的问题，就是这个模数不是质数，所以我们在初始化的时候需要快速计算\(\sum_{i=1}^n i^k\)就不能拉格朗日插值了

那么只好用第二类斯特林数优化了

\[\begin{aligned} \sum\limits_{i=1}^{n}i^K&=\sum_{i=1}^n\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}i^\underline{j}\\ &=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\sum_{i=1}^ni^\underline{j}\\ &=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\frac{

{(n+1)}^\underline{j+1}}{j+1} \end{aligned}\]

然后没有然后了

//minamoto#include
    
     #define R register#define int unsigned int#define IT map
     
      ::iterator#define fp(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
       
        mp;bitset
        
         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\n",ans);    return 0;}