n!+m!=x^p (n>m>1) の方程式を解く！

概要

以下の問題を解きます

この問題は1937年にErdosが解が有限個であることを示しています（ 論文 ）。ですが完全解決はまだされていませんでした。そこでChatGPTにこの問題を完全解決するように依頼したら一発で解いてきて、割と面白かったので正しさの確認がてら記事にすることにしました（ ChatGPTとの会話 ）

前提知識

以下、証明で使う前提知識を挙げておきます

有名なやつです。また、2010年のDusartの 論文 から以下の2つの命題を使います

証明

では証明に入っていきましょう。初めは自明です

正整数の組$(n,m,x,p)$を初めの方程式の解とする

正整数$N,M(N>M)$に対して$N!+M!=M!((M+1)(M+2)...(N-1)N+1))$なので、整数$M\geq1,d\geq0$に対して
$$R_{M,d}=M(M+1)...(M+d) $$
とします。($d=0$のときは$R_{M,d}=1$とします)　$N!+M!=M!(R_{M,N-M}+1)$です

ここからとんでも評価が続きます。便宜上正整数$N$に対して
$$D_N=2q_N-N-1 $$
とします。
補題4より明らかに$0< n-m\leq D_m$です。よって補題5とあわせて
$$P_m \leq R_{m,n-m}+1 \leq R_{m,D_m}+1 $$
が成り立ちます。この$P_m\leq R_{m,D_m}+1$が$m$にのみ依存している条件なので$m$で不等式評価ができます。
なお以降、正整数$a,b$に対して$v_a(b)$で$b$が$a$で割り切れる回数を表すとします。

さて、これで$m$の上限をひとまず与えられました。しかし、さすがにこの上限だと全探索は厳しそうです。よく見ると補題6の不等式評価は緩そうに見えるので、より小さい$m$でも$P_m>R_{m,D_m}+1$となってしまうのではないかと思えてきます。

ここでコンピューターの出番です。以下のコードで $234\leq m<793477$に対して$P_m>R_{m,D_m}+1$であることが確認できます。この不等式が$m$にしかよらないので、補題6くらいの大小評価でも$m$の全探索に持っていけるんですね～頭いい。

      #!/usr/bin/env python3
"""Exact check of P_m > R_{m,D_m}+1 for 234 <= m <= 793476.

Definitions used here:

    q_m     = the least prime q with q > m/2,
    D_m     = 2*q_m - m - 1,
    P_m     = product of primes r with m/2 < r <= m,
    R_m,d   = product_{i=1}^d (m+i), with the empty product equal to 1.

The script verifies, using exact integer arithmetic only, that

    P_m > R_{m,D_m} + 1

for every integer m in the interval 234 <= m <= 793476.
"""

from __future__ import annotations

from bisect import bisect_right
from math import prod

LOWER_M = 234
UPPER_M = 793_476


def sieve(n: int) -> tuple[bytearray, list[int]]:
    """Return an is-prime table and the list of primes up to n."""
    if n < 1:
        return bytearray(b"\x00") * (n + 1), []

    isprime = bytearray(b"\x01") * (n + 1)
    isprime[0] = 0
    isprime[1] = 0

    p = 2
    while p * p <= n:
        if isprime[p]:
            start = p * p
            isprime[start : n + 1 : p] = b"\x00" * (((n - start) // p) + 1)
        p += 1

    primes = [i for i, flag in enumerate(isprime) if flag]
    return isprime, primes


ISPRIME, PRIMES = sieve(UPPER_M)


def q_m(m: int) -> int:
    """Return the least prime q with q > m/2."""
    # For both even and odd m, q > m/2 is equivalent to q > floor(m/2).
    idx = bisect_right(PRIMES, m // 2)
    return PRIMES[idx]


def D_m(m: int) -> int:
    """Return D_m = 2*q_m - m - 1."""
    return 2 * q_m(m) - m - 1


def R_m_D(m: int, D: int) -> int:
    """Return R_{m,D} = product_{i=1}^D (m+i)."""
    return prod(range(m + 1, m + D + 1))


def check_interval() -> tuple[int, int, int, int, int]:
    """Check the inequality throughout the interval.

    Returns:
        checked_count, min_margin_m, min_margin_bits, max_D_m, max_D
    """
    # This variable is maintained as P_m = product_{m/2 < r <= m, r prime} r.
    P = 1

    checked_count = 0
    min_margin_bits: int | None = None
    min_margin_m: int | None = None
    max_D = -1
    max_D_m: int | None = None

    for m in range(2, UPPER_M + 1):
        # When m increases by one, the new endpoint m may enter the interval.
        if ISPRIME[m]:
            P *= m

        # If m is even, the old lower endpoint m/2 is no longer allowed because
        # the interval is open: (m/2, m].  Remove it exactly when it is prime.
        if m % 2 == 0 and ISPRIME[m // 2]:
            P //= m // 2

        if m < LOWER_M:
            continue

        D = D_m(m)
        R_plus_1 = R_m_D(m, D) + 1

        if P <= R_plus_1:
            raise AssertionError(
                f"FAILED at m={m}: D_m={D}, P_m <= R_(m,D_m)+1"
            )

        checked_count += 1

        if D > max_D:
            max_D = D
            max_D_m = m

        # Positive bit-length margin is stronger than necessary, but useful as
        # a compact diagnostic of how far apart the two exact integers are.
        margin_bits = P.bit_length() - R_plus_1.bit_length()
        if min_margin_bits is None or margin_bits < min_margin_bits:
            min_margin_bits = margin_bits
            min_margin_m = m

    assert min_margin_m is not None
    assert min_margin_bits is not None
    assert max_D_m is not None
    return checked_count, min_margin_m, min_margin_bits, max_D_m, max_D


def main() -> None:
    checked, min_m, min_bits, maxD_m, maxD = check_interval()

    expected = UPPER_M - LOWER_M + 1
    if checked != expected:
        raise AssertionError(f"checked {checked} values, expected {expected}")

    print("finite inequality check: PASS")
    print(f"range checked: {LOWER_M} <= m <= {UPPER_M}")
    print(f"m checked: {checked}")
    print(f"smallest bit-length margin: {min_bits} at m={min_m}")
    print(f"largest D_m in this range: D_{maxD_m}={maxD}")


if __name__ == "__main__":
    main()
    

これで$m\leq233$になりました。あとは全探索すればいいですね？補題4があるので計算量もそんな大きくなりません。

      #!/usr/bin/env python3
"""Exact residual search for n! + m! = x^p with 2 <= m <= 233.

For each m, the Erdős--Obláth reduction restricts the search to

    m < n < 2*q_m,

where q_m is the least prime greater than m/2.  Equivalently, with
D_m = 2*q_m - m - 1, it is enough to check

    n = m+1, ..., m+D_m.

This script searches exactly for all representations

    n! + m! = x^p,    x > 1, p > 1,

in that residual range.  All arithmetic is integer arithmetic.
"""

from __future__ import annotations

from bisect import bisect_right
from math import factorial

MAX_M = 233

# A safe exponent bound for the residual range.  Bertrand gives q_m <= m, hence
# n < 2*q_m <= 2*m <= 466.  So factorial(2*MAX_M) is a harmless upper bound.
MAX_EXPONENT_NEEDED = factorial(2 * MAX_M).bit_length() + 1
SIEVE_LIMIT = max(MAX_M, MAX_EXPONENT_NEEDED)


def sieve(n: int) -> tuple[bytearray, list[int]]:
    """Return an is-prime table and the list of primes up to n."""
    if n < 1:
        return bytearray(b"\x00") * (n + 1), []

    isprime = bytearray(b"\x01") * (n + 1)
    isprime[0] = 0
    isprime[1] = 0

    p = 2
    while p * p <= n:
        if isprime[p]:
            start = p * p
            isprime[start : n + 1 : p] = b"\x00" * (((n - start) // p) + 1)
        p += 1

    primes = [i for i, flag in enumerate(isprime) if flag]
    return isprime, primes


ISPRIME, PRIMES = sieve(SIEVE_LIMIT)


def q_m(m: int) -> int:
    """Return the least prime q with q > m/2."""
    idx = bisect_right(PRIMES, m // 2)
    return PRIMES[idx]


def D_m(m: int) -> int:
    """Return D_m = 2*q_m - m - 1."""
    return 2 * q_m(m) - m - 1


def iroot(n: int, k: int) -> int:
    """Return floor(n**(1/k)), using only integer arithmetic."""
    if n < 0:
        raise ValueError("n must be non-negative")
    if k < 1:
        raise ValueError("k must be positive")
    if n < 2:
        return n

    lo = 1
    hi = 1 << ((n.bit_length() + k - 1) // k)
    while lo <= hi:
        mid = (lo + hi) // 2
        power = mid**k
        if power == n:
            return mid
        if power < n:
            lo = mid + 1
        else:
            hi = mid - 1
    return hi


def power_representations(N: int) -> list[tuple[int, int]]:
    """Return all pairs (x,p), p>1, x>1, with N = x**p.

    It is enough first to test prime exponents: if N is a perfect power with
    composite exponent p=ab, then N is also a perfect a-th power.  Once a prime
    exponent is found, we scan all exponents to list every representation.
    """
    if N <= 1:
        return []

    prime_hit = False
    for p in PRIMES:
        if p > N.bit_length():
            break
        x = iroot(N, p)
        if x > 1 and x**p == N:
            prime_hit = True
            break

    if not prime_hit:
        return []

    reps: list[tuple[int, int]] = []
    for p in range(2, N.bit_length() + 1):
        x = iroot(N, p)
        if x > 1 and x**p == N:
            reps.append((x, p))
    return reps


def residual_search() -> tuple[int, int, list[tuple[int, int, int, int, int]]]:
    """Search the residual range.

    Returns:
        pairs_checked, largest_n_checked, solutions

    Each solution is recorded as (n, m, x, p, value), where value = n! + m!.
    """
    solutions: list[tuple[int, int, int, int, int]] = []
    pairs_checked = 0
    largest_n_checked = 0

    factorials = [factorial(i) for i in range(2 * MAX_M + 1)]

    for m in range(2, MAX_M + 1):
        for n in range(m + 1, m + D_m(m) + 1):
            pairs_checked += 1
            largest_n_checked = max(largest_n_checked, n)

            value = factorials[n] + factorials[m]
            for x, p in power_representations(value):
                solutions.append((n, m, x, p, value))

    return pairs_checked, largest_n_checked, solutions


def main() -> None:
    pairs_checked, largest_n, solutions = residual_search()

    expected_solutions = [
        (3, 2, 2, 3, 8),
        (5, 4, 12, 2, 144),
    ]
    if solutions != expected_solutions:
        raise AssertionError(
            f"unexpected solution list: {solutions!r}; expected {expected_solutions!r}"
        )

    print("residual search for n! + m! = x^p: PASS")
    print(f"m range checked: 2 <= m <= {MAX_M}")
    print(f"residual (m,n) pairs checked: {pairs_checked}")
    print(f"largest n checked: {largest_n}")
    print("solutions:")
    for sol in solutions:
        print(sol)


if __name__ == "__main__":
    main()
    

これで解が$(3,2,2,3),(5,4,12,2)$だけであることが証明できました！
前提知識さえ認めてしまえば高校数学範囲でも理解できる証明でしたね

まとめ

ChatGPTはすごい！