洛谷 3768 - 简单的数学题
根据Crash的数字表格，很容易可以将式子化简为
\[\begin{aligned} Ans &= \sum\limits_{i = 1}^n \sum\limits_{j = 1} ij(i, j) \\ &= \sum\limits_{d = 1}^n d^3 \sum\limits_{k = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \mu(k) k^2 \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{kd}\right\rfloor} i \right)^2 \end{aligned}
\]
感觉 \(d, k\) 放在一起式子无法继续化简，主要是有 \(kd\) 存在，故令 \(T = kd\) ，则有
\[Ans = \sum\limits_{T = 1}^n \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{T}\right\rfloor} i \right)^2 T^2 \sum\limits_{d | T} d \mu(\frac{T}{d})
\]
那么考虑整除分块，现在需要处理的是后半部分
通过观察（看题解）可以发现， \(\sum\limits_{d | T} d \mu(\frac{T}{d})\) 可以看成 \(\mu * id\) ，故可替换成 \(\phi(T)\) ，则有
\[Ans = \sum\limits_{T = 1}^n \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{T}\right\rfloor} i \right)^2 T^2 \phi(T)
\]
现在考虑将 \(T^2 \phi(T)\) 部分用杜教筛解决
\[h(n) = \sum\limits_{d | n} d^2 \phi(d) g(\frac{n}{d})
\]
为了消除 \(d^2\) ，令 \(g(\frac{n}{d}) = n^2\) ，则有
\[h(n) = n^3
\]
故得
\[S(n) = \sum\limits_{i = 1}^n i^3 - \sum\limits_{d = 2}^n d^2 S(\left\lfloor\frac{n}{d}\right\rfloor)
\]
又（通过看题解）有一个知识点
\[\sum\limits_{i = 1}^n i^3 = \left( \sum\limits_{i = 1}^n i \right)
\]
那么就可以直接杜教筛了
至于复杂度，将最外围的整除分块与杜教筛看为一体，故复杂度为 \(O (n^{\frac{2}{3}})\)
代码
#include <iostream>
#include <cstdio>
#include <cstring>
#include <tr1/unordered_map>
using namespace std;
typedef long long LL;
const int MAXN = 6e06 + 10;
LL MOD, N;
LL power (LL x, LL p) {
LL cnt = 1;
while (p) {
if (p & 1)
cnt = cnt * x % MOD;
x = x * x % MOD;
p >>= 1;
}
return cnt;
}
LL inv2, inv6;
int prime[MAXN];
int vis[MAXN]= {0};
int pcnt = 0;
LL phi[MAXN]= {0};
LL sumphi[MAXN]= {0};
int MAX = 6e06;
void linear_sieve () {
phi[1] = 1;
for (int i = 2; i <= MAX; i ++) {
if (! vis[i]) {
prime[++ pcnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= pcnt && i * prime[j] <= MAX; j ++) {
vis[i * prime[j]] = 1;
if (! (i % prime[j])) {
phi[i * prime[j]] = phi[i] * 1ll * prime[j] % MOD;
break;
}
phi[i * prime[j]] = phi[i] * 1ll * (prime[j] - 1) % MOD;
}
}
for (int i = 1; i <= MAX; i ++)
sumphi[i] = (sumphi[i - 1] + 1ll * i % MOD * 1ll * i % MOD * phi[i] % MOD) % MOD;
}
tr1::unordered_map<LL, LL> maphi;
inline LL sqr (LL x) {
return x * x % MOD;
}
inline LL eqm (LL n) {
return (n + 1) % MOD * (n % MOD) % MOD * (2 * n % MOD + 1) % MOD * inv6 % MOD;
}
inline LL oseqm (LL n) {
return n % MOD * ((n + 1) % MOD) % MOD * inv2 % MOD;
}
LL phi_sieve (LL n) {
if (n <= MAX)
return sumphi[n];
if (maphi[n])
return maphi[n];
LL total = sqr (oseqm (n));
for (LL l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
total = (total - (eqm (r) - eqm (l - 1) + MOD) % MOD * phi_sieve (n / l) % MOD + MOD) % MOD;
}
return maphi[n] = total;
}
LL Solve () {
LL ans = 0;
for (LL l = 1, r; l <= N; l = r + 1) {
r = N / (N / l);
ans = (ans + sqr (oseqm (N / l)) * ((phi_sieve (r) - phi_sieve (l - 1) + MOD) % MOD) % MOD) % MOD;
}
return ans;
}
int main () {
scanf ("%lld%lld", & MOD, & N);
inv2 = power (2ll, MOD - 2), inv6 = power (6ll, MOD - 2);
MAX = (int) min (1ll * MAX, N), linear_sieve ();
LL ans = Solve ();
cout << ans << endl;
return 0;
}
/*
998244353 2000
*/
/*
1000000007 9786510294
*/
