好友磊哥出了一道数学题，题面是这样的
我们需要高效计算：
$$
G(N) = \sum_{j=1}^N \sum_{i=1}^j \gcd(i,j) \mod 998244353
$$
其中 $N \leq 10^{11}$
已知 $G(10) = 122$
思路
我们不难写出暴力
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| def gcd(i, j):
while j:
i, j = j, i % j
return i
def calculate_G(N, mod):
result = 0
for j in range(1, N + 1):
for i in range(1, j + 1):
result += gcd(i, j)
result %= mod
return result
N = 10
mod = 998244353
result = calculate_G(N, mod)
print(result)
|
但这还是太慢了，问问磊哥怎么说
Solution
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
| import sys
from math import isqrt
from functools import lru_cache
from tqdm import tqdm
MOD = 998244353
INV2 = (MOD + 1) // 2
def compute_G_optimized(N):
# ================== 杜教筛参数优化 ==================
L = int(N ** (2/3)) + 1 # 调整预处理阈值为 N^(2/3)
if L > 10**7: # 内存保护限制
L = 10**7
# ================== 步骤1：线性筛预处理欧拉函数 ==================
print(f"预计算欧拉函数 (范围: {L})...")
phi = list(range(L + 1))
is_prime = [True] * (L + 1)
is_prime[0] = is_prime[1] = False
# 优化筛法实现
with tqdm(total=L-1, desc="线性筛进度") as pbar:
for p in range(2, L + 1):
if is_prime[p]:
phi[p] = p - 1
for multiple in range(p*2, L+1, p):
is_prime[multiple] = False
if (multiple // p) % p == 0:
phi[multiple] = phi[multiple // p] * p
else:
phi[multiple] = phi[multiple // p] * (p - 1)
pbar.update(1)
else:
for q in range(2, p):
if q * p > L:
break
if is_prime[q] and p % q == 0:
break
# 计算前缀和
Phi_prefix = [0] * (L + 1)
with tqdm(total=L, desc="前缀和计算") as pbar:
acc = 0
for i in range(1, L + 1):
acc += phi[i]
Phi_prefix[i] = acc % MOD
pbar.update(1)
# ================== 步骤2：记忆化递归优化 ==================
memo_Phi = {}
call_stack = []
@lru_cache(maxsize=None)
def compute_Phi(m):
if m <= L:
return Phi_prefix[m]
# 杜教筛公式: Φ(m) = m(m+1)/2 - ΣΦ(m//k) [k=2 to m]
if m in memo_Phi:
return memo_Phi[m]
res = (m % MOD) * ((m + 1) % MOD) % MOD * INV2 % MOD
# 分块优化
k = 2
while k <= m:
q = m // k
k_next = m // q + 1
res = (res - (k_next - k) % MOD * compute_Phi(q)) % MOD
k = k_next
memo_Phi[m] = res
return res
# ================== 步骤3：分块并行预加载 ==================
# ================== 智能预加载优化 ==================
print("预加载高频Φ值...")
max_reuse_q = isqrt(N)
precompute_start = L + 1
precompute_end = max_reuse_q
if precompute_end > precompute_start:
precompute_count = precompute_end - precompute_start
with tqdm(total=precompute_count, desc="预加载高频Φ值", unit="q") as pbar:
for q in range(precompute_end, precompute_start-1, -1):
compute_Phi(q)
pbar.update(1)
else:
print(f"所有高频q值({max_reuse_q})已在预处理范围{L}内")
# ================== 步骤4：分块加速计算 ==================
print("主计算阶段...")
G = 0
d = 1
total_blocks = 2 * isqrt(N)
with tqdm(total=total_blocks, desc="处理分块") as pbar:
while d <= N:
q = N // d
d_last = N // q
# 计算块贡献
block_size = d_last - d + 1
sum_d = (d + d_last) * block_size // 2
contribution = (sum_d % MOD) * compute_Phi(q) % MOD
G = (G + contribution) % MOD
# 进度更新
pbar.update(block_size * total_blocks // N)
d = d_last + 1
return G
if __name__ == "__main__":
N = int(sys.stdin.readline())
import time
time_s = time.time()
print(f"G({N}) = {compute_G_optimized(N)}")
time_e = time.time()
print(f"计算时间: {time_e - time_s:.2f}秒")
# print("内存使用: {:.2f}MB".format(sys.getsizeof(compute_G_optimized) / (1024 * 1024)))
print("完成")
|
Output
1
2
3
4
5
6
7
8
9
10
11
| 100000000000
预计算欧拉函数 (范围: 10000000)...
线性筛进度: 7%|█████████▉ | 664579/9999999 [00:10<02:26, 63782.90it/s]
前缀和计算: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10000000/10000000 [00:02<00:00, 3396765.53it/s]
智能预加载高频Φ值...
所有高频q值(316227)已在预处理范围10000000内
主计算阶段...
处理分块: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋| 631290/632454 [00:32<00:00, 19583.69it/s]
G(100000000000) = 551614306
计算时间: 46.13秒
完成
|