Project Euler 408

发表于 分类于 阅读次数： 本文字数： 2.9k 阅读时长 ≈ 3 分钟

Project Euler 408

题目

Admissible paths through a grid

Let’s call a lattice point \((x, y)\) inadmissible if \(x, y\) and \(x+y\) are all positive perfect squares.

For example, \((9, 16)\) is inadmissible, while \((0, 4), (3, 1)\) and \((9, 4)\) are not.

Consider a path from point \((x_1, y_1)\) to point \((x_2, y_2)\) using only unit steps north or east.

Let’s call such a path admissible if none of its intermediate points are inadmissible.

Let \(P(n)\) be the number of admissible paths from \((0, 0)\) to \((n, n)\).

It can be verified that \(P(5) = 252, P(16) = 596994440\) and \(P(1000) \bmod 1000000007 = 341920854\).

Find \(P(10000000) \bmod 1000000007\).

解决方案

若 \((x,y)\) 是障碍点，则存在正整数 \(u,v,w\) 使得\(x=u^2, y=v^2, x+y=w^2.\)代入得到\(u^2+v^2=w^2,\) 也就是说 \((u,v,w)\) 是一个毕达哥拉斯三元组。反过来，任取正整数解 \(u^2+v^2=w^2\)，都给出一个障碍点 \((u^2,v^2)\)。

因此，所有障碍点恰好对应于所有三元组 \((u,v,w)\) 的两条直角边平方：\((u,v,w)\Rightarrow(u^2,v^2)\text{与}(v^2,u^2).\)

又因为我们只关心从 \((0,0)\) 到 \((n,n)\) 的路径，中间点满足 \(0\le x,y\le n\)，而障碍点还要求 \(x,y>0\)，所以我们只需要枚举所有满足\(1\le u^2\le n, 1\le v^2\le n\)的三元组，即\(1\le u,v\le \lfloor\sqrt n\rfloor.\)

令 \(D(A,B)\) 表示从点 \(A=(x_1,y_1)\) 到点 \(B=(x_2,y_2)\)（满足 \(x_1\le x_2,y_1\le y_2\)）的单调路径条数，允许经过障碍点。每条路径由\(\Delta x=x_2-x_1\)次向东，\(\Delta y=y_2-y_1\)次向北组成，总步数 \(\Delta x+\Delta y\)，故

\[ D(A,B)=\binom{\Delta x+\Delta y}{\Delta x}=\binom{\Delta x+\Delta y}{\Delta y}. \]

目标 \(P(n)\) 本质上是“从 \((0,0)\) 到 \((n,n)\) 的所有路径”减去“至少经过一个障碍点的路径”。

把所有障碍点收集为集合 \(\mathcal{S}\)，再加入终点 \(T=(n,n)\)。将 \(\mathcal{S}\cup{T}\) 按字典序排序：先按 \(x\) 递增，若 \(x\) 相同则按 \(y\) 递增，得到序列 \(p_1,p_2,\dots,p_m\)，其中 \(p_m=T\)。

记 \(p_i=(x_i,y_i)\)。对 \(1\le i<m\)（即 \(p_i\in\mathcal{S}\)），定义 \(F_i\) 为从 \((0,0)\) 到 \(p_i\) 的单调路径条数，满足路径的所有中间点都不是障碍点（等价地：路径在到达 \(p_i\) 之前不经过任何其它障碍点）。换句话说，\(F_i\) 统计的是“第一次遇到障碍点恰好发生在 \(p_i\)”的路径数。对于终点 \(p_m=T\)，令 \(F_m\) 表示从 \((0,0)\) 到 \(T\) 的 admissible 路径条数，则所求答案为 \(P(n)=F_m\)。

对任意 \(i\)，一条从 \((0,0)\) 到 \(p_i\) 的单调路径，如果会经过某个障碍点 \(p_j\)，那么必然满足\(x_j\le x_i, y_j\le y_i.\)而由于我们按 \(x\) 递增排序，\(j<i\) 自动保证 \(x_j\le x_i\)，因此只需额外检查 \(y_j\le y_i\)。

从“总路径数”中扣掉“先到某个更早障碍点，再从那里到 \(p_i\) 的路径”即可得到递推。

记从原点到 \(p_i\) 的总路径数为 \(D((0,0),p_i)=\dbinom{x_i+y_i}{x_i}\)。若路径第一次撞到的障碍点是 \(p_j\)（\(j<i\) 且 \(y_j\le y_i\)），那么它的条数为\(F_j\cdot D(p_j,p_i)=F_j\dbinom{(x_i-x_j)+(y_i-y_j)}{x_i-x_j}.\)将所有可能的 \(j\) 累加并扣除。

于是 DP 方程可以写成：

\[ F_i= \left\{\begin{aligned} &\binom{x_1+y_1}{x_1} & & \text{if} \quad i=1\\ &\binom{x_i+y_i}{x_i}-\sum\limits_{\substack{1\le j<i\\y_j\le y_i}} F_j\binom{(x_i-x_j)+(y_i-y_j)}{x_i-x_j} & & \text{if}\quad i>1 \end{aligned}\right. \]

障碍点来自所有 \(u^2+v^2=w^2\) 的解。经典做法是：

  1. 先生成所有原始毕达哥拉斯三元组 \((a,b,c)\)（\(\gcd(a,b,c)=1\)），
  2. 再用倍率 \(t\ge 1\) 扩展为 \((ta,tb,tc)\)，
  3. 对每个 \((u,v,w)\) 记录障碍点 \((u^2,v^2)\) 与对称点 \((v^2,u^2)\)。

本源勾股树给出：从 \((3,4,5)\) 出发，对每个原始三元组 \((a,b,c)\)，以下三种变换会生成新的原始三元组，并且所有原始三元组都能且仅能通过这种方式生成：

\[ \begin{aligned} (a',b',c')&=(a-2b+2c,2a-b+2c,2a-2b+3c),\\ (a',b',c')&=(a+2b+2c,2a+b+2c,2a+2b+3c),\\ (a',b',c')&=(-a+2b+2c,-2a+b+2c,-2a+2b+3c). \end{aligned} \]

对于本题，只需要 \(u,v\le \lfloor\sqrt n\rfloor\)。因此枚举到某个原始三元组 \((a,b,c)\) 时，令 \(L=\max(a,b)\)，所有可用倍率满足\(tL\le \lfloor\sqrt n\rfloor.\)

由此，最终计算出\(P(n)=F_m\)。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from tools import *

N = 10_000_000
MOD = 1000000007

def solve(n: int) -> int:
    r = int_sqrt(n)  # u,v <= r
    co = Combinatorics(n + n, MOD)
    # 1) Generate inadmissible points via Berggren tree of primitive triples.
    pts = set()
    stack = [(3, 4, 5)]
    while stack:
        a, b, c = stack.pop()
        if a > r or b > r:
            continue

        L = a if a > b else b
        tmax = r // L
        for t in range(1, tmax + 1):
            u = t * a
            v = t * b
            x = u * u
            y = v * v
            if x <= n and y <= n:
                pts.add((x, y))
                pts.add((y, x))

        # children primitive triples
        stack.append((a - 2 * b + 2 * c, 2 * a - b + 2 * c, 2 * a - 2 * b + 3 * c))
        stack.append((a + 2 * b + 2 * c, 2 * a + b + 2 * c, 2 * a + 2 * b + 3 * c))
        stack.append((-a + 2 * b + 2 * c, -2 * a + b + 2 * c, -2 * a + 2 * b + 3 * c))

    # Append destination as the last DP target
    pts.add((n, n))
    points = sorted(pts)  # sort by x then y

    # 3) Inclusion-exclusion DP on points.
    m = len(points)
    xs = [p[0] for p in points]
    ys = [p[1] for p in points]
    dp = [0] * m

    for i in range(m):
        xi = xs[i]
        yi = ys[i]
        val = co.C(xi + yi, xi)
        for j in range(i):
            if ys[j] <= yi:
                dx = xi - xs[j]
                dy = yi - ys[j]
                val -= dp[j] * co.C(dx + dy, dx)
                val %= MOD
        dp[i] = val

    # (n,n) is included, and it is the maximum point, so it must be last
    return dp[-1]


if __name__ == "__main__":
    print(solve(N))