Project Euler 844

题目

\(k\)-Markov Numbers

Consider positive integer solutions to

\[a^2+b^2+c^2 = 3abc\]

For example, \((1,5,13)\) is a solution. We define a \(3\)-Markov number to be any part of a solution, so \(1\), \(5\) and \(13\) are all 3-Markov numbers. Adding distinct \(3\)-Markov numbers \(\le 10^3\) would give \(2797\).

Now we define a \(k\)-Markov number to be a positive integer that is part of a solution to:

\[\displaystyle \sum_{i=1}^{k}x_i^2=k\prod_{i=1}^{k}x_i,\quad x_i\text{ are positive integers}\]

Let \(M_k(N)\) be the sum of \(k\)-Markov numbers \(\le N\). Hence \(M_3(10^{3})=2797\), also \(M_8(10^8) = 131493335\).

Define \(\displaystyle S(K,N)=\sum_{k=3}^{K}M_k(N)\). You are given \(S(4, 10^2)=229\) and \(S(10, 10^8)=2383369980\).

Find \(S(10^{18}, 10^{18})\). Give your answer modulo \(1\,405\,695\,061\).

解决方案

令\(K=N=10^{18}\)。固定 \(k\)，并固定前 \(k-1\) 个正整数\(x_1,\dots,x_{k-1}.\)令\(\displaystyle{P=\prod_{i=1}^{k-1}x_i, Q=\sum_{i=1}^{k-1}x_i^2.}\)那么原方程等价于\(x_k^2-kPx_k+Q=0.\)因此若 \((x_1,\dots,x_{k-1},x_k)\) 是正整数解，则另一个根为\(x_k' = kP-x_k=\dfrac{Q}{x_k}.\)这给出一个解的自反变换：

\[ (x_1,\dots,x_{k-1},x_k)\longmapsto (x_1,\dots,x_{k-1},k\prod_{i=1}^{k-1}x_i-x_k). \]

它保持整数性与方程成立。由于方程对分量置换对称，任意排列后仍是解。

当 \(k=3\) 时，方程退化为经典的 Markov 方程；本文所用变换正是 Vieta 跳跃 的直接推广：固定其余坐标，把某一坐标视作二次方程的一个根，并以另一根替换，所得点仍满足原方程。再结合坐标置换可知，解集在这类变换下封闭。

我们把一个解\(x_1, x_2,\cdots, x_k.\)重排成满足\(x_1\le x_2\le\cdots\le x_k.\)的解。

在论文中，命题16(1)和推论17综合给出如下下降规约引理：由于\(x_k\)是最大分量，那么\(1\le x_k' < x_k.\)也就是说：对最大分量做一次 Vieta 替换，会得到另一个正整数解，且最大值严格下降。

接下来我们证明命题：若某个 \(m\le N\) 出现在某个正整数解里，则存在另一个同样包含 \(m\) 的正整数解，且所有分量都 \(\le N\)。

证明：考虑集合\(\mathcal S_m\)表示所有包含\(m\)的正整数解。定义\(\displaystyle\Phi(x)=\max_i \{x_i\}.\)在 \(\mathcal S_m\) 中取一个使 \(\Phi\) 最小的解，记作 \(x^\ast\)，并令 \(\Phi(x^\ast)=M_{\min}\)。

若 \(M_{\min}\le N\)，结论成立。反设 \(M_{\min}>N\)。由于 \(m\le N<M_{\min}\)，所以在 \(x^\ast\) 中，\(m\) 不可能是最大分量。 因此可以对最大分量应用归约引理：把最大分量 \(M_{\min}\) 替换为 \(M^\star\)，得到新解 \(y\)，满足\(\displaystyle\max_i \{y_i\} < M_{\min}.\)并且这个替换只改动最大分量，其他分量不变，所以 \(y\) 仍包含 \(m\)。

于是 \(y\in\mathcal S_m\) 且 \(\Phi(y)<\Phi(x^\ast)\)，这与 \(x^\ast\) 的极小性矛盾。故反设不成立，必有 \(M_{\min}\le N\)。命题得证。

因此，求所有 \(\le N\) 的 \(k\)-Markov 数时，确实只需要遍历\(\max_i x_i\le N\)的可达解，这不会漏掉任何 \(\le N\) 的数，因为任何这类数 \(m\) 若在某个大解里出现，都能通过上述下降过程压缩到一个整体不超过 \(N\) 的解里，且 \(m\) 保留。

由上文的自反变换可知，\((1,1,\dots,1)\) 是显然的正整数解，因为\(\displaystyle{\sum_{i=1}^{k}1^2=k, k\prod_{i=1}^{k}1=k.}\)于是可以从该起点出发，不断对某个坐标执行 Vieta 替换并允许置换，从而遍历可达解。

为了避免重复与排列冗余，可以把一个 \(k\)-元组当作多重集合处理：只记录每个取值出现了多少次。具体地，用\(\displaystyle\left\{v_1^{(c_1)},v_2^{(c_2)},\dots,v_t^{(c_t)}\right\}, v_1<\cdots<v_t,\sum_{i=1}^t c_i=k\)表示状态。这样置换不会产生新状态，去重也更容易。

对某个取值 \(v\) 做一次 Vieta 替换时，设全体乘积为\(\displaystyle M=\prod_{i=1}^{k}x_i,\)则其它坐标的乘积为 \(M/v\)，新值为\(v' = k\cdot(M/v)-v.\)若 \(v'\le 0\) 或 \(v'>N\)，则按上一节的原则剪枝；否则把多重集合中一个 \(v\) 替换成 \(v'\) 得到新状态，并继续搜索。

在搜索过程中维护一个整数集合，把所有出现过的分量值加入集合。搜索结束后，该集合恰为所有 \(\le N\) 的 \(k\)-Markov 数，求和即可得到 \(M_k(N)\)。

当 \(N\) 极大时，若对每个 \(k\) 都直接枚举解图，代价不可接受。为此先固定 \(N\)，定义\(A_k(N)\)表示不超过\(N\)的\(k\)-Markov数的集合。我们当前的目标是刻画 \(A_k(N)\) 随 \(k\) 变化的结构。

把第 \(t\) 步状态写成\(x^{(t)}=(x_1^{(t)},\dots,x_k^{(t)}), t=0,1,\dots,T,\)初态为\(x^{(0)}=(1,\dots,1).\) 记\(\displaystyle M^{(t)}:=\prod_{j=1}^k x_j^{(t)}.\) 一次 Vieta 替换是：选一个坐标下标 \(i_t\in\{1,\dots,k\}\)，并令

\[ x_j^{(t+1)}=\begin{cases} x_j^{(t)}, & j\neq i_t,\\ k\dfrac{M^{(t)}}{x_{i_t}^{(t)}}-x_{i_t}^{(t)}, & j=i_t. \end{cases} \]

于是跳跃路径就是下标序列\(\sigma=(i_0,i_1,\dots,i_{T-1}).\)给定 \(\sigma\) 与初态后，终态唯一确定。若关心终态第 \(r\) 个坐标，定义\(F_{\sigma,r}(k)=x_r^{(T)}.\)

接下来给出命题：对任意固定路径 \(\sigma\)，每一步各分量都是 \(k\) 的整系数多项式；即\(F_{\sigma,r}(k)\in\mathbb Z[k].\)

证明：初态各分量恒为常数 \(1\)，显然属于 \(\mathbb Z[k]\)。若第 \(t\) 步所有分量都在 \(\mathbb Z[k]\)，则\(\displaystyle M^{(t)}=\prod_{j=1}^k x_j^{(t)}\in\mathbb Z[k],\)从而\(\displaystyle k\frac{M^{(t)}}{x_{i_t}^{(t)}}-x_{i_t}^{(t)}\)也是 \(\mathbb Z[k]\) 中元素（这里 \(\displaystyle M^{(t)}/x_{i_t}^{(t)}=\prod_{j\ne i_t}x_j^{(t)}\) 仍为多项式）。因此第 \(t+1\) 步各分量仍为整系数多项式。归纳成立。

记\(r(k)=\lvert A_k(N)\rvert.\)将集合 \(A_k(N)\) 去重后按升序写成\(a_1(k)<a_2(k)<\cdots<a_{r(k)}(k).\)也就是说，\(a_r(k)\) 定义为 \(A_k(N)\) 的第 \(r\) 小元素（仅在 \(1\le r\le r(k)\) 时有定义）。

在足够大的 \(k\) 区间内，存在有限模板族\(\mathcal P_N=\{P_1,\dots,P_T\}\subset \mathbb Z[k],\)使得\(A_k(N)=\{P_j(k):1\le j\le T,P_j(k)\le N\}.\)

定义稳定起点：若整数 \(s\ge 3\) 满足上式对所有 \(k\ge s\) 都成立，则称 \(s\) 为稳定起点。记最小稳定起点为 \(k_0\)。计算上，\(k_0\) 通过自动验证确定：从小到大扫描候选 \(s\)，依次执行窗口稳定性检查、模板恢复（有限差分/插值）、短窗逐点验证与集合级强校验；首个通过全部检验的 \(s\) 即取为 \(k_0\)。

\[ S(K,N)=\sum_{k=3}^{k_0}M_k(N)+\sum_{k=k_0+1}^{K}M_k(N). \]

其中第一段直接枚举，第二段按上文做分段多项式区间和。

为进行大区间求和，定义事件集合\(\mathcal E=\{k:P_i(k)=P_j(k)\lor P_i(k)=N\}.\)事件点把区间 \([k_0+1,K]\) 划分为若干子区间 \(I_\ell=[L_\ell,R_\ell]\)。在每个 \(I_\ell\) 内，满足 \(P_j(k)\le N\) 的模板下标集合及去重规则均固定，记为 \(U_\ell\)，因此\(\displaystyle M_k(N)=\sum_{x\in A_k(N)}x=\sum_{j\in U_\ell}P_j(k) (k\in I_\ell).\)从而\(\displaystyle{\sum_{k=k_0+1}^{K}M_k(N)=\sum_{\ell}\sum_{j\in U_\ell}\sum_{k\in I_\ell}P_j(k).}\)

接下来只需计算内层和。对每个 \(\ell\) 和每个 \(j\in U_\ell\)，定义平移后的序列\(\displaystyle f_{\ell,j}(n)=P_j(L_\ell+n), n=0,1,\dots,m_\ell,m_\ell=R_\ell-L_\ell.\)这样就把模板多项式在区间上的求和转化为从 \(0\) 开始的前缀和：\(\displaystyle\sum_{k=L_\ell}^{R_\ell}P_j(k)=\sum_{n=0}^{m_\ell}f_{\ell,j}(n).\)

若用有限差分表示\(\displaystyle f_{\ell,j}(n)=\sum_{i\ge 0}\Delta^i f_{\ell,j}(0)\binom{n}{i},\)则由恒等式\(\displaystyle\sum_{t=0}^{n}\binom{t}{i}=\binom{n+1}{i+1}\)可得前缀和\(\displaystyle\sum_{t=0}^{n}f_{\ell,j}(t)=\sum_{i\ge 0}\Delta^i f_{\ell,j}(0)\binom{n+1}{i+1}.\)因此\(\displaystyle\sum_{k=L_\ell}^{R_\ell}P_j(k)=\sum_{i\ge 0}\Delta^i f_{\ell,j}(0)\binom{m_\ell+1}{i+1}.\)于是第二段总贡献写为

\[ \sum_{k=k_0+1}^{K}M_k(N) =\sum_{\ell}\sum_{j\in U_\ell}\sum_{k=L_\ell}^{R_\ell}P_j(k) =\sum_{\ell}\sum_{j\in U_\ell}\sum_{i\ge 0}\Delta^i f_{\ell,j}(0)\binom{m_\ell+1}{i+1}. \]

代码

from collections import deque

N = 10 ** 18
K = 10 ** 18
MOD = 1405695061

def markov_numbers(k, limit):
    start = ((1, k),)
    q = deque([start])
    seen = {start}
    nums = {1}

    while q:
        st = q.popleft()
        prod = 1
        for v, c in st:
            prod *= v ** c

        for idx, (v, c) in enumerate(st):
            other_prod = prod // v
            nv = k * other_prod - v
            if nv <= 0 or nv > limit or nv == v:
                continue

            new = list(st)
            if c == 1:
                new.pop(idx)
            else:
                new[idx] = (v, c - 1)

            inserted = False
            for j, (vv, cc) in enumerate(new):
                if vv == nv:
                    new[j] = (vv, cc + 1)
                    inserted = True
                    break
                if vv > nv:
                    new.insert(j, (nv, 1))
                    inserted = True
                    break
            if not inserted:
                new.append((nv, 1))

            new_st = tuple(new)
            if new_st not in seen:
                seen.add(new_st)
                q.append(new_st)
            nums.add(nv)

    return sorted(nums)


def forward_diffs(values):
    diffs = []
    cur = list(values)
    while True:
        diffs.append(cur[0])
        if len(cur) == 1:
            break
        cur = [cur[i + 1] - cur[i] for i in range(len(cur) - 1)]
    while diffs and diffs[-1] == 0:
        diffs.pop()
    return diffs


def eval_poly_binom(diffs, n):
    if n < 0:
        return 0
    res = 0
    c = 1
    for i, di in enumerate(diffs):
        if i == 0:
            c = 1
        else:
            if i > n:
                break
            c = c * (n - i + 1) // i
        res += di * c
    return res


def prefix_sum_poly_binom(diffs, n):
    if n < 0:
        return 0
    m = n + 1
    res = 0
    comb = 1
    for r, di in enumerate(diffs, 1):
        if r > m:
            break
        comb = comb * (m - r + 1) // r
        res += di * comb
    return res


def sum_range_k(diffs, k_l, k_r, shift):
    n_l = k_l - shift
    n_r = k_r - shift
    return prefix_sum_poly_binom(diffs, n_r) - prefix_sum_poly_binom(diffs, n_l - 1)


def find_cut(diffs, k_l, k_r, limit, shift):
    lo, hi = k_l, k_r
    ans = k_l - 1
    while lo <= hi:
        mid = (lo + hi) // 2
        if eval_poly_binom(diffs, mid - shift) <= limit:
            ans = mid
            lo = mid + 1
        else:
            hi = mid - 1
    return ans


def predicted_list(diffs_list, k, k0, limit):
    n = k - k0
    vals = []
    for diffs in diffs_list:
        v = eval_poly_binom(diffs, n)
        if v <= limit:
            vals.append(v)
    return sorted(set(vals))


def auto_k0(N, w_fit=20, v_confirm=10, v_validate=600, k_max=3000):
    cache = {}

    def get(k):
        if k not in cache:
            cache[k] = markov_numbers(k, N)
        return cache[k]

    for k in range(3, k_max + 1):
        get(k)
        if k < 3 + w_fit + v_confirm - 1:
            continue

        s = k - (w_fit + v_confirm) + 1
        R = len(get(s))
        ok = True

        for kk in range(s, k + 1):
            if len(get(kk)) != R:
                ok = False
                break
        if not ok:
            continue

        diffs_list = []
        for r in range(R):
            ys = [get(kk)[r] for kk in range(s, s + w_fit)]
            diffs_list.append(forward_diffs(ys))

        for kk in range(s + w_fit, k + 1):
            n = kk - s
            lst = get(kk)
            for r in range(R):
                if eval_poly_binom(diffs_list[r], n) != lst[r]:
                    ok = False
                    break
            if not ok:
                break
        if not ok:
            continue

        endv = min(k_max, s + v_validate)
        for kk in range(s, endv + 1):
            lst = get(kk)
            pred = predicted_list(diffs_list, kk, s, N)
            if lst != pred:
                ok = False
                break
        if not ok:
            continue

        return s, diffs_list, cache

k0, diffs_list, cache = auto_k0(N)
total = 0
for k in range(3, k0 + 1):
    if k not in cache:
        cache[k] = markov_numbers(k, N)
    total += sum(cache[k])

total += K - k0

for diffs in diffs_list[1:]:
    cut = find_cut(diffs, k0 + 1, K, N, k0)
    if cut >= k0 + 1:
        total += sum_range_k(diffs, k0 + 1, cut, k0)
total %= MOD
print(total)