LOJ-572 「LibreOJ Round #11」Misaka Network 与求和

Description

给定 \(N,k\)，求：

\[\sum_{i=1}^N\sum_{j=1}^N (f(\gcd(i,j)))^k\]

其中 \(f(x)\) 表示 \(x\) 质因子分解结果中次大的质因子，重复的质因数计算多次。

Constraints

\(1\le N,k\le 2\cdot 10^9\)

Solution

记 \(f_k(x)=(f(x))^k\)。推式子：

\[\begin{aligned}\sum_{i}^N\sum_j^N f_k(\gcd(i,j)) &= \sum_g^N f_k(g)\sum_{i}^{N/g}\sum_{j}^{N/g} [\gcd(i,j)=1] \\ &= \sum_g^N f_k(g)\sum_d^{N/g} \mu(d)\left\lfloor\frac{N}{dg}\right\rfloor^2 \\ &= \sum_g^N f_k(g)G\left(\left\lfloor\frac{N}{g}\right\rfloor\right) \\ G(N) &= \sum_{d}^N \mu(d)\lfloor N/d\rfloor^2 \end{aligned} \]

随便什么筛出 \(\mu\)，min-25 用可以用 UOJ Sanrd 的套路筛出 \(f_k\)​。

整除分块套整除分块复杂度 \(O\left( \int_1^\sqrt{N}\sqrt i\right)+O\left( \int_1^\sqrt{N}\sqrt \frac{N}{i} \right) = O(N^{3/4})\)。题中 \(N\le 2\cdot 10^9\)​​，算比较小的，可以直接过。

#include <algorithm>
#include <cmath>
#include <cstdio>
using uint = unsigned int;
using ull = unsigned long long;

inline uint fastpow(uint a, uint b) {
uint r = 1;
for (; b; b >>= 1, a *= a)
if (b & 1) r *= a;
return r;
}

const int N = 1e5;
uint n, K;

namespace min_25 {
uint n, m, tot, cnt, K;
uint D[N], p[N], pK[N];
bool vis[N];

inline uint id(uint x) {
return x <= m ? cnt - x + 1 : n / x;
  }

void sieve(uint n) {
vis[1] = 1;
for (uint i = 1; i <= n; i++) {
if (!vis[i]) p[++tot] = i;
for (uint j = 1; j <= tot && i * p[j] <= n; j++) {
vis[i * p[j]] = 1;
if (i % p[j] == 0) break;
      }
    }
for (uint i = 1; i <= tot; i++)
pK[i] = fastpow(p[i], K);
  }

uint s[N];
void init_s() {
for (uint i = 1; i <= cnt; i++)
s[i] = D[i] - 1;
for (uint j = 1; j <= tot; j++)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
s[i] -= s[id(D[i] / p[j])] - (j - 1);
  }

namespace fk {
uint g[N];
void init_g() {
for (uint i = 1; i <= cnt; i++)
g[i] = s[i];
for (uint j = tot; j; --j)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
for (uint pk = p[j]; (ull)pk * p[j] <= D[i]; pk *= p[j])
g[i] += g[id(D[i] / pk)] - s[id(D[i] / pk)]
+ pK[j] * (s[id(D[i] / pk)] - j + 1);
    }
inline uint get(uint x) {
if (!x) return 0u;
return g[id(x)];
    } 
  }

namespace mu {
uint g[N];
void init_g() {
for (uint i = 1; i <= cnt; i++)
g[i] = -s[i];
for (uint j = tot; j; --j)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
g[i] -= g[id(D[i] / p[j])] + j;
    }
inline uint get(uint x) {
if (!x) return 0u;
return g[id(x)] + 1;
    } 
  }

void load(uint _n, uint _K) {
n = _n, K = _K, m = sqrt(n);
for (uint x = 1; x <= n; x = n / (n / x) + 1)
D[++cnt] = n / x;
sieve(m);

init_s();
mu::init_g();
fk::init_g();
  }
}

signed main() {
scanf("%u%u", &n, &K);
min_25::load(n, K);

auto calc_S = [&](uint n) {
using min_25::mu::get;
uint res = 0;
for (uint l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
res += (get(r) - get(l - 1)) * (n / l) * (n / l);
    }
return res;
  };

auto calc_F = [&](uint n) {
using min_25::fk::get;
uint res = 0;
for (uint l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
res += calc_S(n / l) * (get(r) - get(l - 1));
    }
return res;
  };

printf("%u\n", calc_F(n));
return 0;
}

Solution Plus 1

做完发现有更优做法！考虑迪利克雷卷积：\(h=f_k*\mu\)。那么答案成了：

\[\sum_i^N h(i)\left\lfloor\frac{N}{i}\right\rfloor^2\]

式子非常简洁，尝试直接求 \(h\) 的前缀和。由于不是积性函数，而且要比 \(f_k\) 的魔改方法复杂得多，min-25 筛直接求 \(h\) 有些困难。而 \(h\) 和迪利克雷卷积关系较大，由人类智慧可得：

\[h*I = f_k*\mu *I=f_k*\varepsilon=f_k\]

凑出了卷积形式，考虑杜教筛。记 \(F, H\) 分别为 \(f_k,h\) 的前缀和。那么：

\[F(N)=\sum_{i=1}^N H(\lfloor N/i\rfloor) \Leftrightarrow H(N)=F(N)-\sum_{i=2}^N H(\lfloor N/i\rfloor)\]

杜教筛复杂度 \(O(N^{2/3})\)，总复杂度 \(O\left(\frac{N^{3/4}}{\log N}\right)\)，瓶颈在 min-25 求 \(F\)。

Solution Plus 2

经过 Mr_Spade 指导，可以不要杜教筛！

考虑目标在于快速求 \(H\)​​​。这里有一个船新科技：设阈值 \(n^{2/3}\)​，并以此分治：

• 对于 \(\le n^{2/3}\)​​ 的点值，线性筛出 \(f_k,\mu\)​ 的点值，暴力迪利克雷卷积求出 \(h\)​ 在 \(n^{2/3}\)​​ 内所有的点值，最后前缀和求出 \(H\)。复杂度 \(O(n^{2/3}\log n)\)​。
• 对于 \(>n^{2/3}\)​​ 的点值​，使用之前整除分块的做法算出 \(H\)​ 在较大整除点的点值。复杂度 \(O\left(\int_1^\sqrt[3]N \sqrt\frac{N}{i}\right)=O(n^{2/3})\)。

最后一次整除分块就做完了。复杂度仍然是 \(O\left(\frac{N^{3/4}}{\log N}\right)\)。常数比杜教筛小很多！

#include <algorithm>
#include <cmath>
#include <cstdio>
using uint = unsigned int;
using ull = unsigned long long;

inline uint fastpow(uint a, uint b) {
uint r = 1;
for (; b; b >>= 1, a *= a)
if (b & 1) r *= a;
return r;
}

const int N = 1e5;
const int L = 1e7;
uint n, K, gap;
uint sfk[L], smu[L], sh[L];

namespace min_25 {
uint n, m, tot, cnt, K;
uint D[N], p[L / 10], pK[L / 10];
bool vis[L];

inline uint id(uint x) {
return x <= m ? cnt - x + 1 : n / x;
  }

void sieve(uint n) {
vis[1] = 1;
sfk[1] = 0, smu[1] = 1;
for (uint i = 1; i <= n; i++) {
if (!vis[i]) {
p[++tot] = i;
sfk[i] = 1, smu[i] = -1;
pK[tot] = fastpow(i, K);
      }
for (uint j = 1; j <= tot && i * p[j] <= n; j++) {
vis[i * p[j]] = 1;
if (vis[i]) sfk[i * p[j]] = sfk[i];
else
sfk[i * p[j]] = pK[j];
if (i % p[j] == 0) break;
smu[i * p[j]] = -smu[i];
      }
    }
while (p[tot] > m) --tot;
  }

uint s[N];
void init_s() {
for (uint i = 1; i <= cnt; i++)
s[i] = D[i] - 1;
for (uint j = 1; j <= tot; j++)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
s[i] -= s[id(D[i] / p[j])] - (j - 1);
  }

namespace fk {
uint g[N];
void init_g() {
for (uint i = 1; i <= cnt; i++)
g[i] = s[i];
for (uint j = tot; j; --j)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
for (uint pk = p[j]; (ull)pk * p[j] <= D[i]; pk *= p[j])
g[i] += g[id(D[i] / pk)] - s[id(D[i] / pk)]
+ pK[j] * (s[id(D[i] / pk)] - j + 1);
    }
inline uint get(uint x) {
if (!x) return 0u;
return g[id(x)];
    } 
  }

namespace mu {
uint g[N];
void init_g() {
for (uint i = 1; i <= cnt; i++)
g[i] = -s[i];
for (uint j = tot; j; --j)
for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
g[i] -= g[id(D[i] / p[j])] + j;
    }
inline uint get(uint x) {
if (!x) return 0u;
return g[id(x)] + 1;
    } 
  }

void load(uint _n, uint _K) {
n = _n, K = _K, m = sqrt(n);
for (uint x = 1; x <= n; x = n / (n / x) + 1)
D[++cnt] = n / x;
sieve(std::max(gap, m));

init_s();
mu::init_g();
fk::init_g();
  }
}

signed main() {
scanf("%u%u", &n, &K);
gap = pow(n, 2.0 / 3);
min_25::load(n, K);

for (uint i = 1; i <= gap; i++)
for (uint j = 1; i * j <= gap; j++)
sh[i * j] += sfk[i] * smu[j];
for (uint i = 2; i <= gap; i++)
sh[i] += sh[i - 1];

auto calc_H = [&](uint x) {
auto mu = min_25::mu::get;
auto fk = min_25::fk::get;
if (x <= gap) return sh[x];
uint ret = 0;
for (uint l = 1, r; l <= x; l = r + 1) {
r = x / (x / l);
ret += (mu(r) - mu(l - 1)) * fk(x / l);
    }
return ret;
  };

uint ans = 0, lst = 0;
for (uint l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
uint cur = calc_H(r);
ans += (cur - lst) * (n / l) * (n / l);
lst = cur;
  }

printf("%u\n", ans);
return 0;
}