普通版传送门
加强版传送门
首先我们先来认识一下题目中的\(f(x)\)函数是什么，我们知道，\(x\)有无平方因子，对应的是\(x\)的素因子分解后存在某个素因子指数位大于1，而我们知道\(\mu(x)=1或-1\)时，\(x\)的素因子位数均小于等与\(1\)，因此\(f(x)=\mu(x)^2\)

接下来就是喜闻乐见的推式子环节

\(\sum_{i=1}^n\sum_{j=1}^n(i+j)^k\mu(gcd(i,j))^2*gcd(i,j)=\)
\(\sum_{d=1}^nd\mu(d)^2\sum_{i=1}^n\sum_{j=1}^n(i+j)^k[gcd(i,j)==d]=\)
\(\sum_{d=1}^nd^{k+1}\mu(d)^2\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}(i+j)^k[gcd(i,j)==1]=\)
\(\sum_{d=1}^nd^{k+1}\mu(d)^2\sum_{x=1}^{\frac{n}{d}}\mu(x)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}(i+j)^k[x|gcd(i,j)]=\)
\(\sum_{d=1}^nd^{k+1}\mu(d)^2\sum_{x=1}^{\frac{n}{d}}x^k\mu(x)\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}(i+j)^k=\)
令\(T=dx\)
\(\sum_{T=1}^nT^k\sum_{d|T}d*\mu(d)^2\mu(T/d)\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}(i+j)^k\)

令\(f(T)=\sum_{d|T}d*\mu(d)^2\mu(\frac{T}{d})\)
\(S(n)=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k\)
则原式化为\(\sum_{T=1}^nT^kf(T)S(\frac{n}{T})\)
显然\(S(\frac{n}{T})\)可以整除分块优化，那麽重要的就是在\(O(n)\)时间段内求出其上的所有值。

1:\(T^k\)

可以用线性筛在O(n)时间内取得,令\(tk[n]=\sum_{i=1}^ni^k\)

2:\(f(T)\)

\(f(T)\)毫无疑问为积性函数那麽也可以考虑线性筛
对于一个质数\(p\),我们可以发现\(f(1)=1,f(p)=p-1,f(p^2)=-p,f(p^x)=0,(x>2)\)
因此，
当\((i/p)\%p!=0\)时，\(f(i*p)=f(i)*f(p)\)
当\(i\%p==0\)时，
        当\(((i/p)\%p!=0)f(i*p)=f(i)*f(p)\),否则\(f(i*p)=0\)
这样也就可以\(O(n)\)求出来了

3:\(S(n)\)

这个是最难的，因为\(k\)很大，所以我们不能用二项式定理展开，只能另寻他法。
我们不妨找找各项之间的关系
\(S(n)-S(n-1)=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k-\sum_{i=1}^{n-1}\sum_{j=1}^{n-1}(i+j)^k\)
                                \(=\sum_{i=1}^{n-1}\sum_{j=1}^n(i+j)^k-\sum_{i=1}^{n-1}\sum_{j=1}^{n-1}(i+j)^k+\sum_{j=1}^n(n+j)^k\)
                                \(=\sum_{j=1}^n(n+j)^k+\sum_{i=1}^{n-1}(n+i)^k\)
                                \(=tk[2*n]+tk[2*n-1]-tk[n]*2\)
那麽我们的\(S(n)\)也可以\(O(n)\)递推

综上，将\(\sum_{T=1}^nT^kf(T)\)的前缀和求出来后，我们就可以每次询问在\(O(\sqrt{n})\)的时间内得出

综上复杂度为\(O(n+T*\sqrt{n})\)
这里附上加强版的代码

ll T,n,k;

int p[maxn/10],cnt;
bitset<maxn>vis;

ll f[maxm];//中间部分
ll tk[maxm];//T^k的前缀和
ll S[maxn];

ll qpow(ll a,ll b){
    ll sum=1;
    while(b){if(b&1)sum=sum*a;a=a*a;b>>=1;}
    return sum;
}

void pre(ll m){
    f[1]=1;tk[1]=1;
    rep(i,2,m){
        if(!vis[i]){
            p[++cnt]=i;f[i]=i-1;
            tk[i]=qpow(i,k);
        }
        rep(j,1,cnt){
            if(1LL*p[j]*i>m) break;
            vis[p[j]*i]=1;
            tk[i*p[j]]=1LL*tk[i]*tk[p[j]];
            if(i%p[j]) 
                f[i*p[j]]=f[i]*f[p[j]];
            else{
                if((i/p[j])%p[j]==0) f[i*p[j]]=0;
                else f[i*p[j]]=f[i/p[j]]*-p[j];
                break;
            }
        }
    }
    rep(i,1,n) f[i]=f[i]*tk[i],f[i]=(f[i-1]+f[i]);
    S[1]=tk[2];
    rep(i,1,m)tk[i]=(tk[i]+tk[i-1]);
    rep(i,2,n) S[i]=S[i-1]+tk[i<<1]+tk[(i<<1)-1]-2*tk[i],S[i];
}

int main(){
    T=read(),n=read(),k=read();
    pre(n<<1);
    while(T--){
        n=read();
        ll ans=0;
        for(ll l=1,r;l<=n;l=r+1){
          r=n/(n/l);
          ans=(ans+(f[r]-f[l-1])*S[n/l]);
        }
        print((ans));pts;
    }
    return 0;
}