The Sum of Unitary Divisors

题目描述

A natural number $d$ is a unitary divisor of $n$ if $d$ is a divisor of $n$ and if $d$ and ${n \over d}$ are coprime.

Let $\sigma^{\ast}(n)$ be the sum of the unitary divisors of $n$. For example, $\sigma^{\ast}(1)=1, \ \sigma^{\ast}(2)=3$ and $\sigma^{\ast}(6)=12$.

Let $S(n)=\sum_{i=1}^n \sigma^{\ast}(i)$.

Given $n$, find $S(n) \bmod 2^{64}$.

题意概述

定义$d$为$n$的单位因子当且仅当$d \mid n$且$(d, {n \over d})=1$。

令$\sigma^{\ast}(n)$为$n$的所有单位因子之和，$S(n)=\sum_{i=1}^n \sigma^{\ast}(i)$。给定$n$，求$S(n) \bmod 2^{64}$。

数据范围：$1 \le n \le 5 \times 10^{13}$。

算法分析

推式子

$$
\begin{align}
S(n)&=\sum_{i=1}^n \sigma^{\ast}(i) \\
&=\sum_{i=1}^n i \sum_{j=1}^{\lfloor n/i \rfloor} [(i, j)=1]
\end{align}
$$

利用莫比乌斯反演

$$
\begin{align}
S(n)&=\sum_{d=1}^n \mu(d) \sum_{i=1}^{\lfloor n/d \rfloor} di \sum_{j=1}^{\lfloor n/d^2i \rfloor} 1 \\
&=\sum_{d=1}^n d \cdot \mu(d) \sum_{i=1}^{\lfloor n/d \rfloor} i \cdot \left\lfloor {n \over d^2i} \right\rfloor
\end{align}
$$

显然，当$d^2 \gt n$时不产生贡献。那么

$$
\begin{align}
S(n)=\sum_{d=1}^{\lfloor \sqrt{n} \rfloor} d \cdot \mu(d) \sum_{i=1}^{\lfloor n/d^2 \rfloor} i \cdot \left\lfloor {n \over d^2i} \right\rfloor
\end{align}
$$

问题变成了求$F(n)=\sum_{i=1}^n i \cdot \mu(i)$和$G(n)=\sum_{i=1}^n i \cdot \left\lfloor {n \over i} \right\rfloor$。

对于$F$，我们有$(I \cdot \mu) \ast I=\varepsilon$。所以

$$
\begin{align}
F(n)=1-\sum_{i=2}^n i \cdot F\left(\left\lfloor {n \over i} \right\rfloor\right)
\end{align}
$$

这样$F$和$G$都可以用分块的技巧求出。而求$S$只要枚举$d$，将$\left\lfloor {n \over d^2} \right\rfloor$相等的$d$一起计算。

代码

#include <cmath>
#include <cstdio>

#define int long long

static const int N = 10000005;
static const int M = 50005;
int prime[N], mui[N], g[M];
bool vis[N];

void init() {
  int top = 0; mui[1] = 1;
  for (int i = 2; i < N; ++ i) {
    if (! vis[i]) prime[top ++] = i, mui[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]) mui[k] = -mui[i]; else { mui[k] = 0; break; }
    }
  }
  for (int i = 2; i < N; ++ i) (mui[i] *= i) += mui[i - 1];
}

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

static const int HASH = 23456789;
int numr, h[HASH];
struct Record {
  int n, v, nxt;
} r[HASH];

void add_rec(int n, int v) {
  r[++ numr] = (Record) { n, v, h[n % HASH] }, h[n % HASH] = numr;
}

int get_f(int n) {
  if (n < N) return mui[n];
  for (int i = h[n % HASH]; i; i = r[i].nxt)
    if (r[i].n == n) return r[i].v;
  int ret = 1;
  for (int i = 2, j; i <= n; i = j + 1)
    j = n / (n / i), ret -= (get_sum(j) - get_sum(i - 1)) * get_f(n / i);
  return add_rec(n, ret), ret;
}

int get_g(int n) {
  if (n < M && g[n]) return g[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)) * (n / i);
  return ret;
}

int calc(int n) {
  int ret = 0;
  for (int i = 1, j; i * i <= n; i = j + 1)
    j = sqrt(n / (n / i / i) + 0.01), ret += (get_f(j) - get_f(i - 1)) * get_g(n / i / i);
  return ret;
}

signed main() {
  int T; scanf("%lld", &T), init();
  for (int i = 1; i < M; ++ i) g[i] = get_g(i);
  for (int n; T --;) scanf("%lld", &n), printf("%llu\n", calc(n));
  return 0;
}