原式 $= \displaystyle\sum_{i = 1}^n \sum_{j = 1}^n \lambda(\frac{ij}{\gcd(i, j)})$

$ = \displaystyle\sum_{d = 1}^n \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} [\gcd(i, j) = 1] \lambda(ijd)$

$ = \displaystyle\sum_{d = 1}^n \sum_{p = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(p) \sum_{i = 1}^{\lfloor \frac{n}{dp} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{dp} \rfloor} \lambda(ijdp^2)$

推到这一步看上去不是很好搞了，考虑探究 $\lambda$ 函数的性质。

• $\lambda$ 函数是完全积性函数。

证明：由于不去重，各个因数的贡献独立，于是命题显然成立。

• $\forall n \in Z^+$，$\lambda(n^2) = 1$。

证明：由于 $\lambda$ 函数的值只可能为 $\pm 1$，平方以后只能为 $1$。

原式 $= \displaystyle\sum_{d = 1}^n \sum_{p = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(p) \sum_{i = 1}^{\lfloor \frac{n}{dp} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{dp} \rfloor} \lambda(i) \lambda(j) \lambda(d) \lambda(p^2)$

$ = \displaystyle\sum_{d = 1}^n \lambda(d) \sum_{p = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(p) (\sum_{i = 1}^{\lfloor \frac{n}{dp} \rfloor} \lambda(i))^2$

问题转化为快速求出 $\lambda$ 函数的前缀和。由于多组询问，显然不能直接上 Min_25 筛，但容易发现 $\lambda$ 与 $\mu$ 有一定相似之处，且我们可以杜教筛求出 $\mu$ 函数的前缀和，不妨从这个角度考虑原问题。

不难发现 $\mu(n) = \mu^2(n) \lambda(n)$，为了在 $\lambda$ 和 $\mu$ 间建立卷积关系，不妨令 $f = \frac{\lambda}{\mu}$（狄利克雷卷积意义下的除法），不难发现 $f(n)$ 当且仅当 $n$ 为完全平方数时为 $1$，其他时候为 $0$。

于是我们设定杜教筛阈值 $S$，预处理到 $S$ 并每次杜教筛求出 $\mu$ 的前缀和，求答案时枚举范围内的完全平方数 / 二次数论分块并调用杜教筛即可。然而求解原式时的数论分块套数论分块还是要带来 $O(Tn^{\frac{3}{4}})$ 的时间复杂度，仍然不能通过。

设 $T = dp$，有：

原式 $= \displaystyle\sum_{T = 1}^n (\lambda * \mu)(T) (\sum_{i = 1}^{\lfloor \frac{n}{T} \rfloor} \lambda(i))^2$

我们已经可以快速求出 $\lambda$ 函数的前缀和，于是我们可以杜教筛求出 $\lambda * \mu$ 的前缀和。

设预处理阈值为 $S$，则时间复杂度为 $O(S + \frac{\sum n}{\sqrt{S}})$，当 $S = (\sum n)^{\frac{2}{3}}$ 时取最优时间复杂度为 $O((\sum n)^{\frac{2}{3}})$。

代码：

#include <stdio.h>
#include <math.h>

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef __uint128_t ulll;

typedef struct Division_tag {
	ulll a;
	Division_tag(){}
	Division_tag(ull mod){
		a = (((ulll)1 << 64) + mod - 1) / mod;
	}
} Division;

const int N = 2e7 + 7, M = 2e5 + 7;
ll cur_n, sqrt_n;
short mu[N], liouville[N], f[N];
int prime[N], id1[M], id2[M];
uint mu_sum[N], liouville_sum[N], f_sum[N], h1[M], h2[M], h3[M];
ll number[M];
bool p[N];
Division inv_number[M];

ull operator /(const ull &a, const Division &b){
	return a * b.a >> 64;
}

inline void init1(){
	int cnt = 0;
	p[0] = p[1] = true;
	mu[1] = 1;
	liouville[1] = 1;
	f[1] = 1;
	for (register int i = 2; i < N; i++){
		if (!p[i]){
			prime[++cnt] = i;
			mu[i] = -1;
			liouville[i] = -1;
			f[i] = -2;
		}
		for (register int j = 1; j <= cnt && i * prime[j] < N; j++){
			int t = i * prime[j];
			p[t] = true;
			liouville[t] = -liouville[i];
			if (i % prime[j] == 0){
				mu[t] = 0;
				f[t] = -f[i];
				break;
			}
			mu[t] = -mu[i];
			f[t] = -f[i] * 2;
		}
	}
	for (register int i = 1; i < N; i++){
		mu_sum[i] = mu_sum[i - 1] + mu[i];
		liouville_sum[i] = liouville_sum[i - 1] + liouville[i];
		f_sum[i] = f_sum[i - 1] + f[i];
	}
}

inline ll sqrt(ll n){
	ll ans = sqrt((double)n);
	while (ans * ans <= n) ans++;
	return ans - 1;
}

inline int init2(ll n){
	int id = 0;
	cur_n = n;
	sqrt_n = sqrt(n);
	for (register ll i = 1, j; i <= n; i = j + 1){
		ll tn = n / i;
		j = n / tn;
		id++;
		number[id] = tn;
		inv_number[id] = Division(tn);
		if (tn <= sqrt_n){
			id1[tn] = id;
		} else {
			id2[j] = id;
		}
	}
	return id;
}

inline int get_id(ll n){
	return n <= sqrt_n ? id1[n] : id2[cur_n / n];
}

inline uint sqr(uint n){
	return n * n;
}

int main(){
	int t;
	scanf("%d", &t);
	init1();
	for (register int i = 1; i <= t; i++){
		int id;
		uint ans = 0;
		ll n;
		scanf("%lld", &n);
		id = init2(n);
		for (register int j = id; j >= 1; j--){
			if (number[j] < N){
				h1[j] = mu_sum[number[j]];
				h2[j] = liouville_sum[number[j]];
				h3[j] = f_sum[number[j]];
			} else {
				ll t = sqrt(number[j]);
				h1[j] = 1;
				for (register ll k = 2, l; k <= number[j]; k = l + 1){
					int tn_id = get_id(number[j] / k);
					l = number[j] / inv_number[tn_id];
					h1[j] -= h1[tn_id] * (l - k + 1);
				}
				h2[j] = 0;
				for (register ll k = 1, l; k <= t; k = l + 1){
					int tn_id = get_id(number[j] / k / k);
					l = sqrt(number[j] / inv_number[tn_id]);
					h2[j] += h1[tn_id] * (l - k + 1);
				}
				h3[j] = h2[j];
				for (register ll k = 2, l; k <= number[j]; k = l + 1){
					int tn_id = get_id(number[j] / k);
					l = number[j] / inv_number[tn_id];
					h3[j] -= h3[tn_id] * (l - k + 1);
				}
			}
		}
		for (register ll j = 1, k; j <= n; j = k + 1){
			ll tn = n / j;
			k = n / tn;
			ans += sqr(h2[get_id(tn)]) * (h3[get_id(k)] - h3[get_id(j - 1)]);
		}
		printf("%u\n", ans);
	}
	return 0;
}