题解 - [Luogu P3768] 简单的数学题

强烈谴责这种起无意义标题的做法

题目链接

原始题面

题目描述

由于出题人懒得写背景了,题目还是简单一点好

输入一个整数 \(n\) 和一个整数 \(p\), 你需要求出:

\[ \left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p \]

其中 \(\gcd(a,b)\) 表示 \(a\)\(b\) 的最大公约数

输入输出格式

输入格式

一行两个整数 \(p,n\)

输出格式

一行一个整数表示答案

输入输出样例

输入样例 #1

1
998244353 2000

输出样例 #1

1
883968974

说明

对于 20% 的数据,\(n \leq 1000\)

对于 30% 的数据,\(n \leq 5000\)

对于 60% 的数据,\(n \leq 10^6\), 时限 1s

对于另外 20% 的数据,\(n \leq 10^9\), 时限 3s

对于最后 20% 的数据,\(n \leq 10^{10}\), 时限 4s

对于 100% 的数据,\(5 \times 10^8 \leq p \leq 1.1 \times 10^9\)\(p\) 为质数

解题思路

推式子

\[ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^n ij (i,j)&\equiv\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} ij [(i,j)=1]\pmod p\\ &\equiv\sum_{d=1}^nd^3\sum_{e=1}^{\lfloor\frac{n}{d}\rfloor}e^2\mu(e)\sum_{i=1}^{\lfloor\frac{n}{de}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{de}\rfloor} ij\pmod p\\ &\overset{D=de}{\equiv}\sum_{D=1}^n\left({\lfloor\frac{n}{d}\rfloor\left(\lfloor\frac{n}{d}\rfloor+1\right)\over 2}\right)^2D^2\varphi(D)\pmod p\\ \end{aligned} \]

\(f(n)=n^2\varphi(n)\), 用杜教筛求 \(S(n):=\sum_{i=1}^nf(i)\)

\(g(n)=n^2\), 有

\[ \begin{aligned} S(n)=g(1)S(n)&=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\\ &=\left(\frac{n(n+1)}{2}\right)^2-\sum_{i=2}^ni^2S\left(\left\lfloor\frac{n}{d}\right\rfloor\right) \end{aligned} \]

时间复杂度

首先以 \(O(m)\) 的时间复杂度预处理出 \(S(1..m)\) (线性筛,\(f\) 是积性的所以可以做到)

则时间复杂度为

\[ O\left(O(m)+\int_1^{\sqrt n}O\left(\frac{x}{\sqrt m}\right)\mathrm{d}x\right)=O\left(m+\frac{x}{\sqrt m}\right) \]

选取合适的 \(m\) 后即为 \(O(n^\frac{2}{3})\)

代码参考

Show code

Luogu_P3768view raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
/*
* @Author: Tifa
* @Description: From <https://github.com/Tiphereth-A/CP-archives>
* !!! ATTENEION: All the context below is licensed under a
* GNU Affero General Public License, Version 3.
* See <https://www.gnu.org/licenses/agpl-3.0.txt>.
*/
#include <bits/stdc++.h>
using namespace std;
using data_type = long long;
const int N = 3e6 + 1, P = N / 10 + 1;
data_type sf[N];
unordered_map<data_type, data_type> sum_f;
data_type prime[P], prime2[P], prime3[P], cnt;
bool vis[N];
void init(const data_type &p, const data_type &n = N - 1) {
sf[1] = 1;
for (data_type i = 2; i <= n; ++i) {
if (!vis[i]) {
prime[++cnt] = i;
prime3[cnt] = (prime2[cnt] = i * i % p) * i % p;
sf[i] = (prime3[cnt] + p - prime2[cnt]) % p;
}
for (data_type j = 1; j <= cnt && i * prime[j] <= n; ++j) {
vis[i * prime[j]] = 1;
sf[i * prime[j]] = sf[i] * prime3[j] % p;
if (i % prime[j] == 0) break;
(sf[i * prime[j]] += p - sf[i] * prime2[j]) %= p;
}
}
for (data_type i = 2; i <= n; ++i) (sf[i] += sf[i - 1]) %= p;
}
constexpr data_type g(data_type n, const data_type &p) {
n %= p;
return n * n % p;
}
constexpr data_type inv3(const data_type &p) {
return (p + p * (p % 3 == 1) + 1) / 3;
};
constexpr data_type sum_g(data_type n, const data_type &p) {
n %= p;
return n * (n + 1) / 2 % p * (2 * n + 1) % p * inv3(p) % p;
}
constexpr data_type sum_conv_g_f(data_type n, const data_type &p) {
n %= p;
data_type _ = n * (n + 1) / 2 % p;
return _ * _ % p;
}
data_type get_sum_f_mul_g1(const data_type &n, const data_type &p) {
if (n < N) return sf[n];
if (sum_f[n]) return sum_f[n];
data_type ans = sum_conv_g_f(n, p);
for (data_type l = 2, r = 0; l <= n; l = r + 1) {
r = n / (n / l);
(ans += p - (sum_g(r, p) + p - sum_g(l - 1, p)) % p *
get_sum_f_mul_g1(n / l, p) % p) %= p;
}
return sum_f[n] = ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
data_type n, p;
cin >> p >> n;
init(p);
data_type ans = 0;
for (data_type l = 1, r = 0; l <= n; l = r + 1) {
r = n / (n / l);
(ans += (get_sum_f_mul_g1(r, p) + p - get_sum_f_mul_g1(l - 1, p)) % p *
sum_conv_g_f(n / l, p) % p) %= p;
}
cout << ans << '\n';
return 0;
}