Description
给定 (N,k)，求：
[sum_{i=1}^Nsum_{j=1}^N (f(gcd(i,j)))^k
]
其中 (f(x)) 表示 (x) 质因子分解结果中次大的质因子，重复的质因数计算多次。
Constraints
(1le N,kle 2cdot 10^9)
Solution
记 (f_k(x)=(f(x))^k)。推式子：
[egin{aligned}
sum_{i}^Nsum_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)leftlfloorfrac{N}{dg}
ight
floor^2 \
&= sum_g^N f_k(g)Gleft(leftlfloorfrac{N}{g}
ight
floor
ight) \
G(N) &= sum_{d}^N mu(d)lfloor N/d
floor^2
end{aligned}
]
随便什么筛出 (mu)，min-25 用可以用 UOJ Sanrd 的套路筛出 (f_k)。
整除分块套整除分块复杂度 (Oleft( int_1^sqrt{N}sqrt i ight)+Oleft( int_1^sqrt{N}sqrt frac{N}{i} ight) = O(N^{3/4}))。题中 (Nle 2cdot 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
", calc_F(n));
return 0;
}
Solution Plus 1
做完发现有更优做法！考虑迪利克雷卷积：(h=f_k*mu)。那么答案成了：
[sum_i^N h(i)leftlfloorfrac{N}{i}
ight
floor^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
floor) Leftrightarrow H(N)=F(N)-sum_{i=2}^N H(lfloor N/i
floor)
]
杜教筛复杂度 (O(N^{2/3}))，总复杂度 (Oleft(frac{N^{3/4}}{log N} ight))，瓶颈在 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) 在较大整除点的点值。复杂度 (Oleft(int_1^sqrt[3]N sqrtfrac{N}{i} ight)=O(n^{2/3}))。
最后一次整除分块就做完了。复杂度仍然是 (Oleft(frac{N^{3/4}}{log N} ight))。常数比杜教筛小很多！
#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
", ans);
return 0;
}