Project Euler 412

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

Project Euler 412

题目

Gnomon numbering

For integers \(m, n (0\le n<m)\), let \(L(m,n)\) be an \(m\times m\) grid with the top-right \(n\times n\) grid removed.

For example, \(L(5, 3)\) looks like this:

We want to number each cell of \(L(m,n)\) with consecutive integers \(1, 2, 3, \dots\) such that the number in every cell is smaller than the number below it and to the left of it.

For example, here are two valid numberings of \(L(5,3)\):

Let \(LC(m, n)\) be the number of valid numberings of \(L(m, n)\).

It can be verified that \(LC(3,0) = 42, LC(5,3) = 250250, LC(6,3) = 406029023400\) and \(LC(10,5) \bmod 76543217 = 61251715\).

Find \(LC(10000,5000) \bmod 76543217\).

解决方案

我们把格子看成偏序集：若格子 \(A\) 在格子 \(B\) 的上方或右方，并且可以通过若干次“向下/向左”移动到达 \(B\)，则由题目约束可知必须有 \(A<B\)。这种“在某个方向上严格递增”的填数问题，本质上是在数该偏序集的线性扩张。

为了与经典结论对接，做一个几何变换：将整个图形 顺时针旋转 \(90^\circ\)。旋转后：

  • “小于下方”变成“小于上方”，即向下递增。
  • “小于左方”变成“小于右方”，即向右递增。

于是旋转后的填数规则变为：用 \(1,2,\dots,k\) 填满图形（\(k\) 为格子总数），并且在每一行从左到右递增、每一列从上到下递增。这正是 标准杨表 的定义。

设\(a=m-n, b=n, m=a+b.\)旋转后图形的每一行长度为：

  • 前 \(a\) 行长度为 \(m=a+b\)（完整宽度）；
  • 后 \(b\) 行长度为 \(a\)（缺掉右侧 \(b\) 列）。

因此对应的分拆为

\[ \lambda=(\underbrace{m,m,\dots,m}_{a\text{ times}},\underbrace{a,a,\dots,a}_{b\text{ times}}). \]

格子总数为\(|\lambda|=am+ba=a(m+b)=m^2-n^2.\)记\(k=m^2-n^2.\)那么我们要的就是该形状的标准杨表数：\(LC(m,n)=f^\lambda.\)

对任意杨图 \(\lambda\)，标准杨表数满足著名的 hook-length 公式：

\[ f^\lambda=\frac{|\lambda|!}{\prod_{c\in\lambda} h(c)}, \]

其中 \(h(c)\) 是格子 \(c\) 的 hook length：它等于“该格子本身 + 同行右侧格子数 + 同列下方格子数”。

仍用 \(a=m-n,\ b=n\)。把杨图分成三块（按行列坐标理解）：

  • 左下大矩形：宽 \(a\)，高 \(a+b\)（所有行都有前 \(a\) 列）；
  • 右上附加矩形：宽 \(b\)，高 \(a\)（只有前 \(a\) 行有后 \(b\) 列）；
  • 它们的交叠部分是左上的 \(a\times a\) 正方形。

因此，我们可以直接地按“区域类型”写 hook-length。

左上 \(a\times a\) 正方形（行 \(1\ldots a\)，列 \(1\ldots a\)）

对其中任意格子，右侧还有 \(m-c\) 个格子，下方还有 \(m-r\) 个格子，因此\(h= (m-c)+(m-r)+1.\)

令 \(i=a-r\in[0,a-1]\)，\(j=a-c\in[0,a-1]\)，代入 \(m=a+b\) 得\(h = 2(a+b)-(a-i)-(a-j)+1 = 2b+i+j+1.\)

所以该块的 hook-length 乘积为

\[ H_1=\prod_{i=0}^{a-1}\prod_{j=0}^{a-1}(2b+i+j+1). \]

固定 \(i\)，对 \(j\) 的乘积是一个连续整数段：\(\displaystyle{\prod_{j=0}^{a-1}(2b+i+j+1)=\frac{(2b+i+a)!}{(2b+i)!}=\frac{(m+n+i)!}{(2n+i)!}.}\)

因此 \[ H_1=\prod_{i=0}^{a-1}\frac{(m+n+i)!}{(2n+i)!}. \]

右上 \(a\times b\) 矩形（行 \(1\ldots a\)，列 \(a+1\ldots a+b\)）、左下 \(b\times a\) 矩形（行 \(a+1\ldots a+b\)，列 \(1\ldots a\)）

该矩形中任意格子右侧与下方都只在这块矩形范围内延伸。令 \(t\) 表示该块内的列号（\(t=1\ldots b\)），则右侧个数为 \(b-t\)；下方（在该块内）个数为 \(a-r\)，于是\(h=(b-t)+(a-r)+1.\)

令 \(i=a-r\in[0,a-1]\)，\(j=b-t\in[0,b-1]\)，得到\(h=i+j+1.\)所以

\[ H_2=\prod_{i=0}^{a-1}\prod_{j=0}^{b-1}(i+j+1). \]

类似的，固定 \(i\)：\(\displaystyle{\prod_{j=0}^{b-1}(i+j+1)=\frac{(i+b)!}{i!}=\frac{(n+i)!}{i!}.}\)

因此 \[ H_2=\prod_{i=0}^{a-1}\frac{(n+i)!}{i!}. \]

完全对称可得左下 \(b\times a\) 矩形（行 \(a+1\ldots a+b\)，列 \(1\ldots a\)）块 hook-length 乘积也等于 \(H_2\)。

于是

\[ H = H_1\cdot H_2^2 = \left(\prod_{i=0}^{a-1}\frac{(m+n+i)!}{(2n+i)!}\right) \left(\prod_{i=0}^{a-1}\frac{(n+i)!}{i!}\right)^2. \]

因此

\[ LC(m,n)=\frac{(m^2-n^2)!}{ \left(\prod_{i=0}^{m-n-1}\frac{(m+n+i)!}{(2n+i)!}\right) \left(\prod_{i=0}^{m-n-1}\frac{(n+i)!}{i!}\right)^2 }. \]

代码

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
62
63
64
65
66
67
68
69
M = 10000
N = 5000
mod = 76543217

def factorial_mod_prime(n: int, p: int) -> int:
    """
    Compute n! mod p for 0 <= n < p, where p is prime.
    Uses Wilson complement when n is close to p.
    """
    if n <= 1:
        return 1
    # choose the shorter product
    if n <= (p - 1 - n):
        r = 1
        for x in range(2, n + 1):
            r = (r * x) % p
        return r
    else:
        prod = 1
        for x in range(n + 1, p):
            prod = (prod * x) % p
        # n! ≡ -(prod)^{-1} (mod p)
        return (-pow(prod, p - 2, p)) % p


def lc_mod(m: int, n: int, p: int = 76543217) -> int:
    """
    LC(m,n) mod p for the L-shape (m x m minus top-right n x n),
    with the monotone condition of the problem.
    Assumes p is prime and k = m^2 - n^2 < p (true for the target case).
    """
    assert 0 <= n < m
    a = m - n
    k = m * m - n * n
    assert k < p  # holds for m=10000, n=5000, p=76543217

    # k! mod p via Wilson complement if needed
    fact_k = factorial_mod_prime(k, p)

    # Precompute factorials up to 2m-1 (small: 19999 when m=10000)
    max_t = 2 * m - 1
    fac = [1] * (max_t + 1)
    for i in range(1, max_t + 1):
        fac[i] = (fac[i - 1] * i) % p

    invfac = [1] * (max_t + 1)
    invfac[max_t] = pow(fac[max_t], p - 2, p)
    for i in range(max_t, 0, -1):
        invfac[i - 1] = (invfac[i] * i) % p

    # Hook product H = prod_i ( (m+n+i)!/(2n+i)! ) * ( (n+i)!/i! )^2
    H = 1
    mn = m + n
    two_n = 2 * n
    for i in range(a):
        H = (H * fac[mn + i]) % p
        H = (H * invfac[two_n + i]) % p

        t = (fac[n + i] * invfac[i]) % p
        H = (H * t) % p
        H = (H * t) % p

    # LC = k! * H^{-1} mod p
    return (fact_k * pow(H, p - 2, p)) % p


if __name__ == "__main__":
    print(lc_mod(M,N,mod))