1575 Gcd and Lcm

1575 Gcd and Lcm

∑i=1n∑j=1i∑k=1ilcm(gcd(i,j),gcd(i,k))设f(n)=∑i=1n∑j=1nlcm(gcd(i,n),gcd(j,n))f(p)=3p2−3p+1f(pk)=(2k+1)(p2k−p2k−1)+pk−1\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}\\ i=1∑n​j=1∑i​k=1∑i​lcm(gcd(i,j),gcd(i,k))设f(n)=i=1∑n​j=1∑n​lcm(gcd(i,n),gcd(j,n))f(p)=3p2−3p+1f(pk)=(2k+1)(p2k−p2k−1)+pk−1

/*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;}