625. 最大公约数求和(Gcd Sum)

定义

$$ G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j) $$
已知 \(G(10)=122\)。求 \(G(10^{11})\)，并将结果对 \(998244353\) 取模。

分析：本题的直接枚举会超时，关键在于把双重求和改写成“欧拉函数前缀和 + 整除分块”模型。核心难点是快速计算 \(\Phi(n)=\sum_{k=1}^{n}\varphi(k)\)。

一、数学背景

先固定外层的 \(j\)，考虑

$$ F(j)=\sum_{i=1}^{j}\gcd(i,j) $$

可以用经典分解式得到

$$ F(j)=\sum_{d\mid j} d\cdot \varphi(j/d) $$

于是原式变为

$$ G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\cdot\varphi(j/d) $$

令 \(j=dm\)，得到

$$ G(N)=\sum_{m=1}^{N}\varphi(m)\sum_{d=1}^{\lfloor N/m\rfloor} d =\sum_{m=1}^{N}\varphi(m)\cdot T(\lfloor N/m\rfloor) $$

其中

$$ T(x)=\frac{x(x+1)}{2} $$

这样问题转成： 1. 频繁访问 \(\Phi(x)=\sum_{k\le x}\varphi(k)\)； 2. 对 \(\lfloor N/m\rfloor\) 做整除分块。

二、算法设计

1) 整除分块

对任意固定的商 \(q=\lfloor N/l\rfloor\)，满足该商不变的区间是

$$ r=\left\lfloor\frac{N}{q}\right\rfloor $$

于是区间 \([l,r]\) 的贡献为

$$ (\Phi(r)-\Phi(l-1))\cdot T(q) $$

总区间数约为 \(O(\sqrt N)\)。

2) 递推计算 \(\Phi(n)\)

利用恒等式

$$ \sum_{k=1}^{n}\varphi(k)\left\lfloor\frac{n}{k}\right\rfloor=\frac{n(n+1)}{2} $$

整理可得

$$ \Phi(n)=\frac{n(n+1)}{2}-\sum_{l=2}^{n}(r-l+1)\cdot \Phi(\lfloor n/l\rfloor) $$

同样用整除分块处理上式，并用记忆化缓存重复子问题。

3) 小范围筛预处理

对较小的 \(n\)（文中实现取 5_000_000）先用线性筛求 \(\varphi\) 与其前缀和，作为递推的“地基”，大幅减少递归深度与重复开销。

三、复杂度分析

四、代码实现与说明

"""Project Euler 625 - Gcd Sum."""

from __future__ import annotations

from functools import lru_cache
from time import perf_counter

MOD = 998244353
TARGET_N = 10**11
# Benchmarked on this machine to keep the full solve under 60 seconds.
SIEVE_LIMIT = 5_000_000


def build_phi_prefix(limit: int, mod: int) -> list[int]:
    """Build prefix sums of Euler's totient up to `limit` modulo `mod`."""
    phi = [0] * (limit + 1)
    is_composite = bytearray(limit + 1)
    primes: list[int] = []
    phi[1] = 1

    for i in range(2, limit + 1):
        if not is_composite[i]:
            primes.append(i)
            phi[i] = i - 1
        for p in primes:
            v = i * p
            if v > limit:
                break
            is_composite[v] = 1
            if i % p == 0:
                phi[v] = phi[i] * p
                break
            phi[v] = phi[i] * (p - 1)

    prefix = [0] * (limit + 1)
    prefix[1] = 1
    for i in range(2, limit + 1):
        prefix[i] = (prefix[i - 1] + phi[i]) % mod
    return prefix


def build_solver(mod: int = MOD, sieve_limit: int = SIEVE_LIMIT):
    """Prepare reusable closures for Phi-prefix and G(n) queries."""
    prefix_small = build_phi_prefix(sieve_limit, mod)
    inv2 = pow(2, mod - 2, mod)

    @lru_cache(maxsize=None)
    def phi_prefix_sum(n: int) -> int:
        if n <= sieve_limit:
            return prefix_small[n]

        total = (n % mod) * ((n + 1) % mod) % mod
        total = (total * inv2) % mod
        l = 2
        while l <= n:
            q = n // l
            r = n // q
            count = r - l + 1
            total = (total - (count % mod) * phi_prefix_sum(q)) % mod
            l = r + 1
        return total

    def gcd_sum(n: int) -> int:
        ans = 0
        l = 1
        while l <= n:
            q = n // l
            r = n // q
            phi_interval = (phi_prefix_sum(r) - phi_prefix_sum(l - 1)) % mod
            tri = (q % mod) * ((q + 1) % mod) % mod
            tri = (tri * inv2) % mod
            ans = (ans + phi_interval * tri) % mod
            l = r + 1
        return ans

    return gcd_sum


def main() -> None:
    build_start = perf_counter()
    gcd_sum = build_solver()
    build_elapsed = perf_counter() - build_start
    print(f"build_elapsed_seconds = {build_elapsed:.3f}")

    sample = gcd_sum(10)
    print(f"G(10) mod {MOD} = {sample}")
    assert sample == 122, f"Sample check failed: got {sample}, expected 122"

    t0 = perf_counter()
    result = gcd_sum(TARGET_N)
    elapsed = perf_counter() - t0
    print(f"G(10^11) mod {MOD} = {result}")
    print(f"solve_elapsed_seconds = {elapsed:.3f}")


if __name__ == "__main__":
    main()

代码说明（按执行顺序）：

常量区

build_phi_prefix

build_solver

phi_prefix_sum

gcd_sum

main

 

本站所有文章如非特别声明皆为原创，转载请与作者本人联系并注明出处(原文链接)。 如对文章有任何疑问或需指出错误，可以此页面提出问题，我会尽量及时回复。

←返回主页