Project Euler 801

题目

\(x^y \equiv y^x\)

The positive integral solutions of the equation \(x^y=y^x\) are \((2,4)\), \((4,2)\) and \((k,k)\) for all \(k > 0\).

For a given positive integer \(n\), let \(f(n)\) be the number of integral values \(0 < x,y \leq n^2-n\) such that \[x^y\equiv y^x \pmod n\] For example, \(f(5)=104\) and \(f(97)=1614336\).

Let \(S(M,N)=\sum f(p)\) where the sum is taken over all primes \(p\) satisfying \(M\le p\le N\).

You are given \(S(1,10^2)=7381000\) and \(S(1,10^5) \equiv 701331986 \pmod{993353399}\).

Find \(S(10^{16}, 10^{16}+10^6)\). Give your answer modulo \(993353399\).

解决方案

在模 \(p\) 意义下，若 \(p\mid x\) 且 \(y>0\)，则 \(x^y\equiv 0\pmod p\)。记\(n=p-1\)。

• 若 \(p\mid x\) 且 \(p\mid y\)，则 \(x^y\equiv y^x\equiv 0\pmod p\) 恒成立。
• 若恰有一个能被 \(p\) 整除，例如 \(p\mid x\) 且 \(p\nmid y\)，则 \(x^y\equiv 0\) 而 \(y^x\not\equiv 0\)，不可能相等。

区间 \([1,pn]\) 中恰好有 \(n\) 个 \(p\) 的倍数（分别是 \(p,2p,\dots,np\)），因此贡献为\(\#\{x:1\le x\le pn,p\mid x\}=n,\Rightarrow\#\{(x,y):p\mid x,p\mid y\}=n^2.\)

于是\(f(p)=n^2+\#\{(x,y):1\le x,y\le pn,p\nmid x,p\nmid y,x^y\equiv y^x\pmod p\}.\)

下面只讨论 \(p\nmid x,y\) 的部分。

因为 \(p\) 是素数，乘法群 \((\mathbb Z/p\mathbb Z)^\times\) 是阶为 \(n\) 的循环群，取其原根 \(g\)。那么对任意 \(p\nmid x\)，存在唯一 \(i\in\{0,1,\dots,n-1\}\) 使得\(x\equiv g^i\pmod p.\)同理对 \(y\)，存在 \(j\) 使得 \(y\equiv g^j\pmod p\)。

又因为 \(g^n\equiv 1\pmod p\)，对任意整数指数 \(t\) 有\((g^i)^t\equiv g^{i(t\bmod n)}\pmod p.\)

因此在 \(p\nmid x,y\) 时，有

\[ \begin{aligned} &x^y\equiv y^x\pmod p\\ \iff&(g^i)^y\equiv (g^j)^x\pmod p\\ \iff&g^{i(y\bmod n)}\equiv g^{j(x\bmod n)}\pmod p\\ \iff&i(y\bmod n)\equiv j(x\bmod n)\pmod n. \end{aligned} \]

关键在于：现在比较的是模 \(n\) 的线性同余。

固定一个非零剩余类 \(a\in\{1,2,\dots,n\}\)，考虑所有满足\(x\equiv a\pmod p, 1\le x\le pn\)的 \(x\)。它们恰好是\(x=a+tp, t=0,1,\dots,n-1,\)共 \(n\) 个。由于 \(p=n+1\equiv 1\pmod n\)，于是\(x\bmod n=(a+tp)\bmod n\equiv (a+t)\bmod n.\)当 \(t\) 取遍 \(0,1,\dots,n-1\) 时，\((a+t)\bmod n\) 会把 \({0,1,\dots,n-1}\) 恰好覆盖一遍。

也就是说：当固定 \(x\) 在模 \(p\) 下的非零同余类后，区间 \([1,pn]\) 中所有满足该同余类的整数，其在模 \(n\) 下的取值会恰好覆盖 \(\mathbb Z_n\) 的每一个元素一次（即构成一个全排列）。同理这对 \(y\) 也成立。

因此，对于固定的离散对数 \((i,j)\)，我们只需要统计\(T(i,j)=\#\{(e,f)\in \mathbb Z_n^2:i f\equiv j e\pmod n\},\)

其中 \(e=x\bmod n\)、\(f=y\bmod n\)。由于对每个 \(i\) 都有唯一 \(x\bmod p\) 与之对应，所以 \(p\nmid x,y\) 的总贡献就是\(\displaystyle{\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} T(i,j).}\)

接下来计算 \(T(i,j)\)。为统计同余方程\(if \equiv je \pmod n\)在 \((e,f)\in \mathbb Z_n^2\) 中的解数，定义群同态\(\varphi:\mathbb Z_n^2\to \mathbb Z_n, \varphi(e,f)=if-je\pmod n.\)

则满足 \(if\equiv je\pmod n\) 的 \((e,f)\) 恰好等价于 \(\varphi(e,f)\equiv 0\)，因此\(T(i,j)=\#\{(e,f)\in\mathbb Z_n^2:\varphi(e,f)=0\}=\#\ker\varphi.\)这里\(\ker\varphi={(e,f)\in\mathbb Z_n^2:\ if-je\equiv 0\pmod n}\)就是该同余方程的解集。

另一方面，\(\varphi\) 的像为\(\operatorname{Im}\varphi=\{\varphi(e,f): (e,f)\in\mathbb Z_n^2\}=\{if-je \bmod n:\ e,f\in\mathbb Z_n\}.\)由于 \(e,f\) 可以自由取遍 \(\mathbb Z_n\)，上式中的所有值正好是由 \(i\) 与 \(j\) 在加法群 \(\mathbb Z_n\) 中所能生成的子群，即\(\operatorname{Im}\varphi=\langle i,j\rangle=\{ai+bj \bmod n:\ a,b\in\mathbb Z_n\}.\)

在循环群 \(\mathbb Z_n\) 中，任意子群都形如 \(d\mathbb Z_n\)，其中 \(d\mid n\)，并且由 \(i,j\) 生成的子群的最小正生成元为\(d=\gcd(i,j,n)\Rightarrow\operatorname{Im}\varphi=d\mathbb Z_n=\{0,d,2d,\dots,n-d\}.\)

于是像的大小为\(|\operatorname{Im}\varphi|=\left|d\mathbb Z_n\right|=\frac{n}{d}.\)利用有限群同态的基本计数关系\(|\mathbb Z_n^2|=|\ker\varphi|\cdot|\operatorname{Im}\varphi|,\)即可得到\(|\ker\varphi|=\dfrac{|\mathbb Z_n^2|}{|\operatorname{Im}\varphi|}=\dfrac{n^2}{n/d}=nd.\)因此

\[ T(i,j)=nd=n\gcd(i,j,n). \]

从而 \(p\nmid x,y\) 的贡献为\(\displaystyle{\sum_{i,j}T(i,j)=n\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\gcd(i,j,n),}\)最终有

\[ f(p)=n^2+n\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\gcd(i,j,n),n=p-1. \]

令 \(d\mid n\)，考虑满足 \(\gcd(i,j,n)=d\) 的 \((i,j)\) 对的个数。写 \(i=du,j=dv\)，其中 \(u,v\) 在模 \(m=n/d\) 的意义下取值，则条件等价于\(\gcd(u,v,m)=1.\)而这正是 Jordan’s totient function 的定义计数：\(J_2(m)=\#\{(u,v)\in\mathbb Z_m^2:\gcd(u,v,m)=1\}.\)

因此\(\#\{(i,j):\gcd(i,j,n)=d\}=J_2\left(\dfrac{n}{d}\right).\)于是\(\displaystyle{\sum_{i,j}\gcd(i,j,n)=\sum_{d\mid n} d\cdot J_2\left(\frac{n}{d}\right)=\sum_{t\mid n} \frac{n}{t}J_2(t),}\)其中做了换元 \(t=n/d\)。

代回\(f(p)\)，得到

\[ f(p)=n^2+n\sum_{t\mid n}\frac{n}{t}J_2(t)=n\left(n+\sum_{t\mid n}\frac{n}{t}J_2(t)\right),n=p-1. \]

定义\(\displaystyle{A(n)=\sum_{t\mid n}\frac{n}{t}J_2(t),}\)则\(f(p)=n(n+A(n)), n=p-1.\)由于 \(J_2\) 是乘法函数，\(A(n)\) 也是乘法函数，因此只需在素幂 \(n=q^e\) 上求闭式。根据\(J_2\)的定义，对 \(e\ge 1\)，有

\[ J_2(q^k)= \begin{cases} 1, & k=0,\\ q^{2k}-q^{2k-2}, & k\ge 1. \end{cases} \]

于是\(\displaystyle{A(q^e)=\sum_{k=0}^e q^{e-k}J_2(q^k)=q^e+\sum_{k=1}^e q^{e-k}(q^{2k}-q^{2k-2})=q^e+\sum_{k=1}^e (q^{e+k}-q^{e+k-2}).}\)把两段等比求和整理后可化简为

\[ A(q^e)=q^{e-1}(q^{e+1}+q^e-1). \]

因此若\(n\)有分解\(\displaystyle{n=\prod_{i=1}^{t} q_i^{e_i},}\)则\(\displaystyle{A(n)=\prod_{i=1}^{t} A(q_i^{e_i})=\prod_{i=1}^{t} q_i^{e_i-1}(q_i^{e_i+1}+q_i^{e_i}-1).}\)最终有

\[ f(p)=(p-1)((p-1)+A(p-1)). \]

代码

from tools import is_prime, Factorizer, factorization

L = 10 ** 16
R = 10 ** 16 + 10 ** 6
MOD = 993353399

def primes_in_interval(L: int, R: int):
    W = R - L + 1
    seg = bytearray(b"\x01") * W
    if L == 0:
        seg[0] = 0
        if W > 1:
            seg[1] = 0
    if L == 1:
        seg[0] = 0
    primes = Factorizer(10 ** 6 * 3).primes
    for q in primes:
        start = (L + q - 1) // q * q
        for m in range(start, R + 1, q):
            seg[m - L] = 0
    res = []
    for i in range(W):
        if seg[i]:
            x = L + i
            if is_prime(x):
                res.append(x)
    return res

def A_prime_power(q: int, e: int) -> int:
    qm = q % MOD
    t = pow(qm, e - 1, MOD) * (pow(qm, e + 1, MOD) + pow(qm, e, MOD) - 1) % MOD
    return t

def A_from_factor(n: int) -> int:
    fac = factorization(n)
    res = 1
    for q, e in fac:
        res = (res * A_prime_power(q, e)) % MOD
    return res

def f_of_prime(p: int) -> int:
    n = p - 1
    a = A_from_factor(n)
    nm = n % MOD
    return nm * (nm + a) % MOD

def solve(L, R):
    ps = primes_in_interval(L, R)
    ans = 0
    for p in ps:
        ans += f_of_prime(p)
        ans %= MOD
    return ans

print(solve(L, R))