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

Counting Divisors (Square)

题目描述

Let $\sigma_0(n)$ be the number of positive divisors of $n$.
For example, $\sigma_0(1)=1, \sigma_0(2)=2$ and $\sigma_0(6)=4$.
Let $S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$.
Given $N$, find $S_2(N)$.

题意概述

令$\sigma_0(n)$表示$n$的正因数个数,$S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$。给定$N$,求$S_2(N)$。
数据范围:$1 \le N \le 10^{12}$。

算法分析

为了方便,我们令$f(n)=\sigma_0(n^2)$,$\omega(n)$表示$n$的不同质因子个数。
设$n=\prod_{k=1}^{\omega(n)} p_k^{a_k}$,其中$p_k$表示不同质数。
先来推$f(n)$。对于$n$的每个因数$d$,将$d$的某些质因子$p_k$的次数加上$a_k$得到$d’$,有$d’ \nmid n, d’ \mid n^2$。这样的方案数有$(2^{\omega(d)}-1)$种,再算上$d$本身,可得
$$
\begin{align}
f(n)=\sum_{d \mid n} 2^{\omega(d)}
\end{align}
$$
化简$2^{\omega(d)}$。它等价于$d$的每个质因数要么选一个,要么不选。那么它就等于$d$的因数中square-free的数的个数。因此
$$
\begin{align}
2^{\omega(d)}=\sum_{i \mid d} \mu^2(i)
\end{align}
$$
回到原式
$$
\begin{align}
S_2(n)=\sum_{i=1}^n \sum_{j \mid i} \sum_{k \mid j} \mu^2(k)
\end{align}
$$
还是毫无头绪。我们继续化简$f(n)$。考虑它的卷积形式$f=\mu^2 \ast 1 \ast 1=\mu^2 \ast \sigma_0$。所以
$$
\begin{align}
S_2(n)&=\sum_{i=1}^n \sum_{j \mid i} \mu^2(j) \sigma_0\left({i \over j}\right) \\
&=\sum_{ij \le n} \mu^2(i) \sigma_0(j)
\end{align}
$$
那么只要用反演求出$\mu^2$的前缀和、用分块求出$\sigma_0$的前缀和,就可以利用分块的技巧了。具体如下
$$
\begin{align}
\sum_{i=1}^n \mu^2(i)&=\sum_{i=1}^n [i\text{最大平方因子}=1] \\
&=\sum_{i=1}^n \sum_{j^2 \mid i} \mu(j) \\
&=\sum_{i=1}^{\lfloor \sqrt{n} \rfloor} \mu(i) \left\lfloor {n \over i^2} \right\rfloor \\
\sum_{i=1}^n \sigma_0(i)&=\sum_{i=1}^n \left\lfloor {n \over i} \right\rfloor
\end{align}
$$

代码

#include <cmath>
#include <cstdio>
#include <cstring>

void read(long long &t) {
  char c; while ((c = getchar()) < '0' || c > '9') ; t = c - '0';
  while ((c = getchar()) >= '0' && c <= '9') (t *= 10) += c - '0';
}

int N;
char mui[100000000];
int d[100000000], sum[100000000], prime[10000000];
bool vis[100000000];
long long n, sigma[100000000];

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

inline long long get_mui(long long n) {
  if (n < N) return sum[n];
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += mui[i] * (n / i / i);
  return ret;
}

inline long long get_sigma(long long n) {
  if (n < N) return sigma[n];
  long long ret = 0;
  for (register long long i = 1, j; i <= n; i = j + 1)
    j = n / (n / i), ret += (j - i + 1) * (n / i);
  return ret;
}

long long calc(long long n) {
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += get_sigma(n / i);
  for (register long long i = q + 1, j, cur, last = sum[q]; i <= n; i = j + 1, last = cur)
    j = n / (n / i), cur = get_mui(j), ret += (cur - last) * get_sigma(n / i);
  return ret;
}

signed main() {
  long long T; read(T);
  if (T > 800) N = 20000; else N = 100000000;
  init();
  while (T --) read(n), printf("%lld\n", calc(n));
  return 0;
}

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();
  for (int n; T --;) scanf("%lld", &n), printf("%llu\n", calc(n));
  return 0;
}