原版传送门 & 加强版传送门

题意：

\(T\) 组数据，求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^k\mu^2(\gcd(i,j))\gcd(i,j)\)

弱化版中 \(T=1\)，\(n \leq 5 \times 10^6\)

强化版中 \(T=10^4\)，\(n \leq 10^7\)

推式子：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^k\mu^2(\gcd(i,j))\gcd(i,j)
\]

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^k\mu^2(d)\times d[\gcd(i,j)=d]
\]

\[\sum\limits_{d=1}^n\mu^2(d)d^{k+1}\times\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}(i+j)^k[\gcd(i,j)=1]
\]

\[\sum\limits_{d=1}^n\mu^2(d)d^{k+1}\times\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}(i+j)^k\sum\limits_{p|\gcd(i,j)}\mu(p)
\]

\[\sum\limits_{d=1}^n\mu^2(d)d^{k+1}\times\sum\limits_{p|\gcd(i,j)}\mu(p)p^k\times\sum\limits_{i=1}^{\frac{n}{dp}}\sum\limits_{j=1}^{\frac{n}{dp}}(i+j)^k
\]

令 \(s(x)=\sum\limits_{i=1}^x\sum\limits_{j=1}^x(i+j)^k\)

\[\sum\limits_{dp \leq n}\mu^2(d)d^{k+1}\mu(p)p^ks(\frac{n}{dp})
\]

\[\sum\limits_{t=1}^ns(\frac{n}{t})\sum\limits_{d|t}\mu^2(d)d^{k+1}\mu(\frac{t}{d})(\frac{t}{d})^k
\]

\[\sum\limits_{t=1}^ns(\frac{n}{t})t^k\sum\limits_{d|t}d\mu^2(d)\mu(\frac{t}{d})
\]

令 \(f(t)=\sum\limits_{d|t}d\mu^2(d)\mu(\frac{t}{d})\)

\[\sum\limits_{t=1}^ns(\frac{n}{t})f(t)t^k
\]

预处理 \(s(x)\) 和 \(g(x)=\sum\limits_{i=1}^xf(i)i^k\)，就可以使用整除分块在 \(\sqrt{n}\) 的时间内求出。

接下来我们的问题就是如何求出 \(s(x)\) 和 \(g(x)\)。

首先预处理 \(i^k\) 肯定是需要的。不过一个个快速幂会超时，不过发现 \(i^k\) 是一个积性函数，欧拉筛解决，这部分时间复杂度 \(\pi(n)\log k\)。

对于 \(s(x)\)，直接求肯定不太容易，我们不妨转化为枚举 \(t \in [2,2x]\)，观察 \(t^k\) 被贡献了几次。

例如 \(x=4\)，\(s(4)=1 \times 2^k+2 \times 3^k+3 \times 4^k+4 \times 5^k+3 \times 6^k+2 \times 7^k+1 \times 8^k\)

观察每一项前面的系数，发现了什么。呈阶梯状分部！

预处理 \(w(x)=\sum\limits_{i=1}^xi^x \times (x-i+1)\)，稍微画个图就会发现 \(s(x)=w(2x)-2w(x)\)

接下来是 \(g(x)\)，要求出 \(g(x)\)，肯定要先求出 \(f(x)\)。

由于 \(f(x)\) 是 \(x\mu^2(x)\) 与 \(mu(x)\) 的狄利克雷卷积，而两项都是积性函数，故 \(f(x)\) 也是积性函数。

因此可以用欧拉筛求出 \(f(x)\)。假设 \(x\) 质因数分解里面有一项 \(p^q\)，我们考虑这一项对答案的贡献。

若 \(q=1\)，\(1\times\mu^2(1)\times\mu(p)+p\times\mu^2(p)\times\mu(1)=p-1\)

若 \(q=2\)，\(1\times\mu^2(1)\times\mu(p^2)+p\times\mu^2(p)\times\mu(p)+p^2\times\mu^2(p^2)\times\mu(1)=-p\)

若 \(q\geq 3\)，那么 \(d\) 与 \(\frac{p^q}{d}\) 中必定有一项的幂 \(\geq 2\)，故 \(f(p^q)=0\)。

\(s(x)\) 与 \(g(x)\) 都求出来了，本题也就迎刃而解了。

/*
Contest: -
Problem: P6222
Author: tzc_wk
Time: 2020.9.24
*/
#include <bits/stdc++.h>
using namespace std;
#define fi            first
#define se            second
#define pb            push_back
#define fz(i,a,b)    for(int i=a;i<=b;i++)
#define fd(i,a,b)    for(int i=a;i>=b;i--)
#define foreach(it,v) for(__typeof(v.begin()) it=v.begin();it!=v.end();it++)
#define all(a)        a.begin(),a.end()
#define fill0(a)    memset(a,0,sizeof(a))
#define fill1(a)    memset(a,-1,sizeof(a))
#define fillbig(a)    memset(a,0x3f,sizeof(a))
#define y1            y1010101010101
#define y0            y0101010101010
#define int unsigned int
typedef pair<int,int> pii;
typedef long long ll;
inline int read(){
    int x=0,neg=1;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') neg=-1;
        c=getchar();
    }
    while(isdigit(c)) x=x*10+c-'0',c=getchar();
    return x*neg;
}
int T=read(),N=read(),k=read();
inline int qpow(int x,int e){
    int ans=1;
    while(e){
        if(e&1) ans=ans*x;
        x=x*x;e>>=1;
    } return ans;
}
int pr[20000005],pcnt=0,f[20000005],p[20000005];
bool vis[20000005];
inline void prework(int n){
    f[1]=p[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){pr[++pcnt]=i;f[i]=i-1;p[i]=qpow(i,k);}
        for(int j=1;j<=pcnt&&pr[j]*i<=n;j++){
            p[pr[j]*i]=p[i]*p[pr[j]];vis[pr[j]*i]=1;
            if(i%pr[j]) f[pr[j]*i]=f[pr[j]]*f[i];
            else{
                int lft=i/pr[j];
                if(lft%pr[j]) f[pr[j]*i]=-pr[j]*f[lft];
                else f[pr[j]*i]=0;
                break;
            }
        }
    }
    fz(i,1,n) f[i]=f[i-1]+f[i]*p[i];
    fz(i,1,n) p[i]+=p[i-1];
    fz(i,1,n) p[i]+=p[i-1];
}
inline int sum(int x){
    return (p[x<<1]-(p[x]<<1));
}
signed main(){
    prework(N<<1);
    while(T--){
        int n=read(),ans=0;
        for(int l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            ans+=(f[r]-f[l-1])*sum(n/l);
        }
        cout<<ans<<endl;
    }
    return 0;
}