LOV3 -> OV3

好友磊哥出了一道数学题，题面是这样的

我们需要高效计算：

$G(N) = \sum_{j=1}^N \sum_{i=1}^j \gcd(i,j) \mod 998244353$

其中 $N \leq 10^{11}$

已知 $G(10) = 122$

思路

我们不难写出暴力

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

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

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秒
完成