1575 Gcd and Lcm

1575 Gcd and Lcm

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^i\sum _{k=1}^{i}lcm(gcd(i,j),gcd(i,k)) \\ 设f(n)=\sum _{i=1}^n\sum _{j=1}^{n}lcm(gcd(i,n),gcd(j,n)) \\ f(p)=3p^2-3p+1 \\ f(p^k)=(2k+1)(p^{2k}-p^{2k-1})+p^{k-1} \end{gathered}$$

/*
  Author : lifehappy
*/
#include <bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10;

typedef unsigned int uint;

const uint inv3 = 2863311531;

uint prime[N], id1[N], id2[N], m, cnt, T;

uint a[N], g1[N], sum1[N], g2[N], sum2[N], g3[N], sum3[N], n;

bool st[N];

int ID(int x) {
    return x <= T ? id1[x] : id2[n / x];
}

void init() {
    T = sqrt(n + 0.5);
    m = 0;
    for(uint l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        a[++m] = n / l;
        if(a[m] <= T) id1[a[m]] = m;
        else id2[n / a[m]] = m;
        g1[m] = (uint)(1ull * a[m] * (a[m] + 1) * (2 * a[m] + 1) / 2 * inv3) - 1;
        g2[m] = (uint)(1ull * a[m] * (a[m] + 1) / 2) - 1;
        g3[m] = a[m] - 1;
    }
    for(int j = 1; j <= cnt && prime[j] <= T; j++) {
        for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
            g1[i] -= prime[j] * prime[j] * (g1[ID(a[i] / prime[j])] - sum1[j - 1]);
            g2[i] -= prime[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]);
            g3[i] -= g3[ID(a[i] / prime[j])] - sum3[j - 1];
        }
    }
}

uint quick_pow(uint a, int n) {
    uint ans = 1;
    while(n) {
        if(n & 1) ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

uint f(uint p, uint k) {
    return uint((2 * k + 1) * (quick_pow(p, 2 * k) - quick_pow(p, 2 * k - 1)) + quick_pow(p, k - 1));
}

uint solve(int n, int m) {
    if(n < prime[m]) return 0;
    uint ans = (3 * g1[ID(n)] - 3 * g2[ID(n)] + g3[ID(n)]) - (3 * sum1[m - 1] - 3 * sum2[m - 1] + sum3[m - 1]);
    for(uint j = m; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {
        for(uint i = prime[j], k = 1; 1ll * i * prime[j] <= n; i *= prime[j], k++) {
            ans += f(prime[j], k) * solve(n / i, j + 1) + f(prime[j], k + 1);
        }
    }
    return ans;
}

int main() {
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    for(uint i = 2; i < N; i++) {
        if(!st[i]) {
            prime[++cnt] = i;
            sum1[cnt] = sum1[cnt - 1] + i * i;
            sum2[cnt] = sum2[cnt - 1] + i;
            sum3[cnt] = sum3[cnt - 1] + 1;
        }
        for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {
            st[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                break;
            }
        }
    }
    int T;
    cin >> T;
    while(T--) {
        cin >> n;
        init();
        cout << solve(n, 1) + 1 << "\n";
    }
	return 0;
}