Project Euler 638

题目

Weighted lattice paths

Let \(P_{a,b}\) denote a path in a \(a\times b\) lattice grid with following properties:

• The path begins at \((0,0)\) and ends at \((a,b)\).
• The path consists only of unit moves upwards or to the right; that is the coordinates are increasing with every move.

Denote \(A(P_{a,b})\) to be the area under the path. For the example of a \(P_{4,3}\) path given below, the area equals 6.

Define \(G(P_{a,b},k)=k^{A(P_{a,b})}\). Let \(C(a,b,k)\) equal the sum of \(G(P_{a,b},k)\) over all valid paths in a \(a\times b\) lattice grid.

You are given that

• \(C(2,2,1)=6\)
• \(C(2,2,2)=35\)
• \(C(10,10,1)=184\,756\)
• \(C(15,10,3) \equiv 880\,419\,838 \mod 1\,000\,000\,007\)
• \(C(10\,000,10\,000,4) \equiv 395\,913\,804 \mod 1\,000\,000\,007\)

Calculate \(\displaystyle\sum_{k=1}^7 C(10^k+k, 10^k+k,k)\) . Give your answer modulo \(1\,000\,000\,007\)

解决方案

把路径写成长度 \(a+b\) 的字符串，包含 \(a\) 个 \(R\) 与 \(b\) 个 \(U\)。

当路径在高度 \(y\) 处走一步向右时，这一段横边覆盖的面积增量恰好是 \(y\)（对应这一列下方有 \(y\) 个单位格）。而此时的 \(y\) 等于这一步之前出现的 \(U\) 的个数。

因此 \(A(P_{a,b})\) 等于字符串中所有满足一个 \(U\) 出现在一个 \(R\) 之前的有序对 \((U,R)\) 的数量。

考虑所有从 \((0,0)\) 到 \((a,b)\) 的路径的最后一步：

• 若最后一步是 \(U\)，来自某条到 \((a,b-1)\) 的路径，面积不变，贡献为 \(C(a,b-1,k)\)。
• 若最后一步是 \(R\)，来自某条到 \((a-1,b)\) 的路径；此时最后这一步在高度 \(b\)，面积增加 \(b\)。权重多乘 \(k^b\)，贡献为 \(k^b C(a-1,b,k)\)。

于是对 \(a,b\ge 1\)：

\[ C(a,b,k)=C(a,b-1,k)+k^bC(a-1,b,k), \]

边界条件为\(C(a,0,k)=C(0,b,k)=1.\)

对固定的 \(a\)，定义关于 \(b\) 的生成函数\(\displaystyle F_a(x)=\sum_{b=0}^{\infty} C(a,b,k)x^b.\)把递推乘 \(x^b\) 并对 \(b\ge 1\) 求和：\(\displaystyle\sum_{b\ge1}C(a,b,k)x^b=\sum_{b\ge1}C(a,b-1,k)x^b+\sum_{b\ge1}k^b C(a-1,b,k)x^b.\)左边是 \(F_a(x)-1\)。第一项为 \(xF_a(x)\)。第二项为 \(F_{a-1}(kx)-1\)。因此\(F_a(x)-1=xF_a(x)+F_{a-1}(kx)-1\Longrightarrow(1-x)F_a(x)=F_{a-1}(kx).\)

又 \(\displaystyle F_0(x)=\sum_{b\ge0}1\cdot x^b=\dfrac1{1-x}\)。递推展开得到 \[ F_a(x)=\frac1{(1-x)(1-kx)(1-k^2x)\cdots(1-k^a x)}. \]

于是\(\displaystyle C(a,b,k)=[x^b] \prod_{i=0}^{a}\frac1{1-k^i x}.\)

记\(\displaystyle C_a(x)=\prod_{i=0}^{a}\frac1{1-k^i x}.\)注意到\(\displaystyle C_a(kx)=\prod_{i=0}^{a}\frac1{1-k^{i+1}x}=\frac1{1-k^{a+1}x}\prod_{i=1}^{a}\frac1{1-k^{i}x}.\)因此有恒等式\((1-x)C_a(x)=(1-k^{a+1}x)C_a(kx).\)

对两边取 \([x^b](b\ge 1)\)：\([x^b]C_a(x)-[x^{b-1}]C_a(x)= k^b[x^b]C_a(x)-k^{a+b}[x^{b-1}]C_a(x).\)令 \([x^b]C_a(x)=C(a,b,k)\)，整理得\((1-k^b)C(a,b,k)=(1-k^{a+b})C(a,b-1,k),\)即\(C(a,b,k)=\dfrac{1-k^{a+b}}{1-k^b}C(a,b-1,k),C(a,0,k)=1.\) 连乘得到 \[ C(a,b,k)=\prod_{i=1}^{b}\frac{1-k^{a+i}}{1-k^i}. \]

这正是高斯二项式的经典形式：\(\displaystyle C(a,b,k)=\binom{a+b}{b}_{k}.\)

当 \(k\to 1\) 时，\(\dfrac{1-k^{t}}{1-k}\to t,\)于是 \(q\)-二项式退化为普通二项式：\(C(a,b,1)=\dbinom{a+b}{b}.\)

因此，计算时用迭代维护幂即可：顺序计算 \(k^1,k^2,\dots,k^{2a}\)；前 \(a\) 项累乘分母，后 \(a\) 项累乘分子；最终乘上分母的逆元。

代码

K = 7
MOD = 1000000007

def nCr_mod(n, r):
    if r < 0 or r > n:
        return 0
    r = r if r < n - r else n - r
    num = 1
    den = 1
    for i in range(1, r + 1):
        num = num * (n - r + i) % MOD
        den = den * i % MOD
    return num * pow(den, MOD - 2, MOD) % MOD

def C_diag(a, k):
    if k == 1:
        return nCr_mod(2 * a, a)
    den = 1
    num = 1
    p = 1
    for _ in range(a):
        p = p * k % MOD
        den = den * (1 - p) % MOD
    for _ in range(a):
        p = p * k % MOD
        num = num * (1 - p) % MOD
    return num * pow(den, MOD - 2, MOD) % MOD

def solve(K):
    ans = 0
    for k in range(1, K + 1):
        a = 10 ** k + k
        ans = (ans + C_diag(a, k)) % MOD
    return ans

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