GCD Extreme (Hard)

题目描述

Let $G(n)=\sum_{i=1}^n \sum_{j=i+1}^n (i, j)$.
For example, $G(1)=0, \ G(2)=(1, 2)=1, \ G(3)=(1, 2)+(1, 3)+(2, 3)=3$.
Given $N$, find $G(N)$ modulo $2^{64}$.

题意概述

令$G(n)=\sum_{i=1}^n \sum_{j=i+1}^n (i, j)$。给定$N$，求$G(N) \bmod 2^{64}$。

数据范围：$1 \le N \le 235711131719$。

算法分析

推式子

$$
\begin{align}
G(n)&=\sum_{i=1}^n \sum_{j=i+1}^n (i, j)=\sum_{i=2}^n \sum_{j=1}^{i-1} (i, j) \\
&=\sum_{d=1}^n d \sum_{i=2}^{\lfloor n/d \rfloor} \sum_{j=1}^{i-1} [(i, j)=1] \\
&=\sum_{d=1}^n d \sum_{i=2}^{\lfloor n/d \rfloor} \varphi(i)
\end{align}
$$

令$S(n)=\sum_{i=1}^n \varphi(i)$。由于$\varphi \ast 1=I$，因此

$$
\begin{align}
\sum_{i=1}^n (\varphi \ast 1)(i)&={n(n+1) \over 2} \\
&=\sum_{i=1}^n \sum_{j \mid i} \varphi(j) \\
&=\sum_{ij \le n} \varphi(j) \\
&=\sum_{i=1}^n \sum_{j=1}^{\lfloor n/i \rfloor} \varphi(j) \\
&=\sum_{i=1}^n S\left(\left\lfloor {n \over i} \right\rfloor\right)
\end{align}
$$

那么

$$
\begin{align}
S(n)&={n(n+1) \over 2}-\sum_{i=2}^n S\left(\left\lfloor {n \over i} \right\rfloor\right) \\
G(n)&=\sum_{d=1}^n d \cdot \left(S\left(\left\lfloor {n \over d} \right\rfloor\right)-1\right)
\end{align}
$$

接下来就可以利用分块的技巧计算。

代码

#include <map>
#include <cstdio>

#define int long long

static const int N = 10000000;
int prime[N], phi[N];
bool vis[N];

void init() {
  int top = 0; phi[1] = 1;
  for (int i = 2; i < N; ++ i) {
    if (! vis[i]) prime[top ++] = i, phi[i] = i - 1;
    for (int j = 0; j < top; ++ j) {
      int k = i * prime[j]; if (k >= N) break;
      vis[k] = true;
      if (i % prime[j]) phi[k] = phi[i] * (prime[j] - 1); else { phi[k] = phi[i] * prime[j]; break; }
    }
  }
  for (int i = 2; i < N; ++ i) phi[i] += phi[i - 1];
}

int get_sum(int n) {
  return n & 1 ? ((n + 1) >> 1) * n : (n >> 1) * (n + 1);
}

std :: map <int, int> mp;

int get_phi(int n) {
  if (n < N) return phi[n];
  if (mp.count(n)) return mp[n];
  int ret = get_sum(n);
  for (int i = 2, j; i <= n; i = j + 1)
    j = n / (n / i), ret -= (j - i + 1) * get_phi(n / i);
  return mp[n] = ret;
}

int calc(int n) {
  int ret = 0;
  for (int i = 1, j; i <= n; i = j + 1)
    j = n / (n / i), ret += (get_sum(j) - get_sum(i - 1)) * (get_phi(n / i) - 1);
  return ret;
}

signed main() {
  int T; scanf("%lld", &T), init();
  while (T --) {
    int n; scanf("%lld", &n), printf("%llu\n", calc(n));
  }
  return 0;
}