洛谷 P6222

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

题意：
(T) 组数据，求 (sumlimits_{i=1}^nsumlimits_{j=1}^n(i+j)^kmu^2(gcd(i,j))gcd(i,j))
弱化版中 (T=1)，(n leq 5 imes 10^6)
强化版中 (T=10^4)，(n leq 10^7)

推式子：

[sumlimits_{i=1}^nsumlimits_{j=1}^n(i+j)^kmu^2(gcd(i,j))gcd(i,j) ]

[sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^n(i+j)^kmu^2(d) imes d[gcd(i,j)=d] ]

[sumlimits_{d=1}^nmu^2(d)d^{k+1} imessumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{n}{d}}(i+j)^k[gcd(i,j)=1] ]

[sumlimits_{d=1}^nmu^2(d)d^{k+1} imessumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{n}{d}}(i+j)^ksumlimits_{p|gcd(i,j)}mu(p) ]

[sumlimits_{d=1}^nmu^2(d)d^{k+1} imessumlimits_{p|gcd(i,j)}mu(p)p^k imessumlimits_{i=1}^{frac{n}{dp}}sumlimits_{j=1}^{frac{n}{dp}}(i+j)^k ]

令 (s(x)=sumlimits_{i=1}^xsumlimits_{j=1}^x(i+j)^k)

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

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

[sumlimits_{t=1}^ns(frac{n}{t})t^ksumlimits_{d|t}dmu^2(d)mu(frac{t}{d}) ]

令 (f(t)=sumlimits_{d|t}dmu^2(d)mu(frac{t}{d}))

[sumlimits_{t=1}^ns(frac{n}{t})f(t)t^k ]

预处理 (s(x)) 和 (g(x)=sumlimits_{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 imes 2^k+2 imes 3^k+3 imes 4^k+4 imes 5^k+3 imes 6^k+2 imes 7^k+1 imes 8^k)
观察每一项前面的系数，发现了什么。呈阶梯状分部！
预处理 (w(x)=sumlimits_{i=1}^xi^x imes (x-i+1))，稍微画个图就会发现 (s(x)=w(2x)-2w(x))
接下来是 (g(x))，要求出 (g(x))，肯定要先求出 (f(x))。
由于 (f(x)) 是 (xmu^2(x)) 与 (mu(x)) 的狄利克雷卷积，而两项都是积性函数，故 (f(x)) 也是积性函数。
因此可以用欧拉筛求出 (f(x))。假设 (x) 质因数分解里面有一项 (p^q)，我们考虑这一项对答案的贡献。
若 (q=1)，(1 imesmu^2(1) imesmu(p)+p imesmu^2(p) imesmu(1)=p-1)
若 (q=2)，(1 imesmu^2(1) imesmu(p^2)+p imesmu^2(p) imesmu(p)+p^2 imesmu^2(p^2) imesmu(1)=-p)
若 (qgeq 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;
}