「Luogu 3768」简单的数学题【杜教筛】

Time Limit: 6 Sec Memory Limit: 256 MB

Description

求下式的值 mod p\text{mod } p

i=1nj=1nijgcd(i,j)\sum_{i = 1}^n \sum_{j = 1}^n i \cdot j \cdot \gcd(i, j)

Input

一行两个整数 p,np, n

Output

一行一个整数,表示答案。

Sample Input

1
998244353 2000

Sample Output

1
883968974

Constraints

对于 100%100\% 的数据,满足 n1010,n \leqslant 10^{10}, 5×108p1.1×1095 \times 10^8 \leq p \leqslant 1.1 \times 10^9

Solution

S(n)=1+2++nS(n) = 1 + 2 + \cdots + n,则有:

i=1nj=1nijgcd(i,j)=d=1nd3i=1n/dj=1n/dij[gcd(i,j)=1]=d=1nd3i=1n/dj=1n/dijki,kjμ(k)=d=1nd3k=1nμ(k)k2S2(nkd)=T=1nS2(nT)dTd3(Td)2μ(Td)=T=1nS2(nT)T2dTdμ(Td)=T=1nS2(nT)T2φ(T)\begin{aligned} & \sum_{i = 1}^n \sum_{j = 1}^n ij \gcd(i, j) \\ = & \sum_{d = 1}^n d^3 \sum_{i = 1}^{\lfloor n / d \rfloor} \sum_{j = 1}^{\lfloor n / d \rfloor} ij [\gcd(i, j) = 1] \\ = & \sum_{d = 1}^n d^3 \sum_{i = 1}^{\lfloor n / d \rfloor} \sum_{j = 1}^{\lfloor n / d \rfloor} ij \sum_{k \mid i, k \mid j} \mu(k) \\ = & \sum_{d = 1}^n d^3 \sum_{k = 1}^n \mu(k) k^2 S^2 \left( \left \lfloor \frac{n}{kd} \right \rfloor \right) \\ = & \sum_{T = 1}^n S^2 \left( \left \lfloor \frac{n}{T} \right \rfloor \right) \sum_{d \mid T} d^3 \left( \frac{T}{d} \right)^2 \mu \left( \frac{T}{d} \right) \\ = & \sum_{T = 1}^n S^2 \left( \left \lfloor \frac{n}{T} \right \rfloor \right) T^2 \sum_{d \mid T} d \mu \left( \frac{T}{d} \right) \\ = & \sum_{T = 1}^n S^2 \left( \left \lfloor \frac{n}{T} \right \rfloor \right) T^2 \varphi(T) \end{aligned}

考虑如何处理 f(i)=i2φ(i)f(i) = i^2 \varphi(i),要找到能够方便计算下式的函数 gg

g(1)S(n)=i=1n(fg)(i)i=2ng(i)S(ni)g(1)S(n) = \sum_{i = 1}^n (f * g)(i) - \sum_{i = 2}^n g(i) S \left( \frac{n}{i} \right)

注意到 diφ(d)=i\sum_{d \mid i} \varphi(d) = i,考虑令 g(x)=x2g(x) = x^2,显然此时 (fg)(i)=i3(f * g)(i) = i^3

另外,可以用下面的公式求自然数三次方前缀和:

13+23++n3=(1+2++n)21^3 + 2^3 + \cdots + n^3 = (1 + 2 + \cdots + n)^2

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int maxn = 4641600;
ll n, P;
int tot, p[maxn], phi[maxn];
map<ll, int> ans;

ll G(ll x) {
ll y = x + 1;
x & 1 ? y >>= 1 : x >>= 1;
return x % P * (y % P) % P;
}

ll S(ll x) {
return G(x) * G(x) % P;
}

ll F(ll x) {
ll y = x + 1, z = x << 1 | 1;
x & 1 ? y >>= 1 : x >>= 1;
if (!(x % 3)) x /= 3;
else if (!(y % 3)) y /= 3;
else z /= 3;
return x % P * (y % P) % P * (z % P) % P;
}

void sieve() {
fill(p, p + maxn, 1), phi[1] = 1;
for (int i = 2; i < maxn; i++) {
if (p[i]) p[++tot] = i, phi[i] = i - 1;
for (int j = 1; j <= tot && i * p[j] < maxn; j++) {
p[i * p[j]] = 0;
if (!(i % p[j])) { phi[i * p[j]] = phi[i] * p[j]; break; }
else phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
for (int i = 1; i < maxn; i++) {
phi[i] = (phi[i - 1] + 1LL * phi[i] * i % P * i % P) % P;
}
}

ll calc(ll x) {
if (x < maxn) return phi[x];
if (ans.count(x)) return ans[x];
ll sum = S(x);
for (ll i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
sum = (sum - calc(x / i) * (F(j) - F(i - 1)) % P + P) % P;
}
return ans[x] = sum;
}

int main() {
scanf("%lld %lld", &P, &n), sieve();
ll sum = 0;
for (ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
sum = (sum + S(n / i) * (calc(j) - calc(i - 1) + P) % P) % P;
}
printf("%lld\n", sum);
return 0;
}