Project Euler 382

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

Project Euler 382

题目

Generating polygons

A polygon is a flat shape consisting of straight line segments that are joined to form a closed chain or circuit. A polygon consists of at least three sides and does not self-intersect.

A set \(S\) of positive numbers is said to generate a polygon \(P\) if: - no two sides of \(P\) are the same length,

  • the length of every side of \(P\) is in \(S\), and

  • \(S\) contains no other value.

For example:

The set \(\{3, 4, 5\}\) generates a polygon with sides \(3, 4,\) and \(5\) (a triangle).

The set \(\{6, 9, 11, 24\}\) generates a polygon with sides \(6, 9, 11,\) and \(24\) (a quadrilateral).

The sets \(\{1, 2, 3\}\) and \(\{2, 3, 4, 9\}\) do not generate any polygon at all.

Consider the sequence s, defined as follows: - \(s_1 = 1, s_2 = 2, s_3 = 3\)

  • \(s_n = s_{n-1} + s_{n-3}\) for \(n > 3\).

Let \(U_n\) be the set \(\{s_1, s_2, \dots, s_n\}\). For example, \(U_{10} = \{1, 2, 3, 4, 6, 9, 13, 19, 28, 41\}\).

Let \(f(n)\) be the number of subsets of \(U_n\) which generate at least one polygon.

For example, \(f(5) = 7, f(10) = 501\) and \(f(25) = 18635853\).

Find the last \(9\) digits of \(f(10^{18})\).

解决方案

令\(N=10^{18}\)。将 \(S\) 的元素按从小到大记为 \(a_1<a_2<\dots<a_k\)（\(k=|S|\)）。经典多边形不等式表明：存在（非退化）简单多边形以这些数为边长，当且仅当最长边严格小于其余边之和，即\(a_k<a_1+\cdots+a_{k-1}.\)等价地写成\(\displaystyle{2\max(S)<\sum_{x\in S}x,}\)且显然还需要 \(|S|\ge 3\)。该条件可由三角不等式对折线的推广得到。在本题中，\(U_n\) 内元素互不相同，因此子集天然满足边长互异。

首先我们需要证明两个恒等式。

第一个是证明\(\displaystyle{\sum_{i=1}^{n}s_i=s_{n+3}-3 (n\ge 1).}\)使用归纳法进行证明。当\(n=1\) 时，\(s_4=4\)，右边 \(=4-3=1=s_1\)。

若对 \(n-1\) 成立，则\(\displaystyle{\sum_{i=1}^{n}s_i=\left(\sum_{i=1}^{n-1}s_i\right)+s_n=(s_{n+2}-3)+s_n.}\)由递推式得 \(s_{n+3}=s_{n+2}+s_n\)，因此\((s_{n+2}-3)+s_n=s_{n+3}-3.\)证毕。

由此立刻得到

\[ \sum_{i=1}^{n-3}s_i=s_n-3<s_n. \]

因此：若某个可行集合 \(S\subseteq U_n\) 的最大元素为 \(s_n\)，那么仅靠 \({s_1,\dots,s_{n-3}}\) 的全部和也不足以“压过” \(s_n\)，于是 \(S\) 必须包含 \(s_{n-1}\) 或 \(s_{n-2}\)（或两者都含）。

第二个是证明\(s_n=s_{n-2}+s_{n-3}+s_{n-4}.\)可见，由递推式\(s_n=s_{n-1}+s_{n-3}, s_{n-1}=s_{n-2}+s_{n-4}\)相加即得。

接下来考虑使用动态规划求解这道题目。令 \(f(n)\) 为 \(U_n\) 的可行子集数。

根据上面“必须含 \(s_{n-1}\) 或 \(s_{n-2}\)”的结论，把“最大边为 \(s_n\) 的可行子集”分成两类：

  • \(V(n)\)：包含 \(s_n\)，但 不包含 \(s_{n-1}\) 的可行子集数；
  • \(W(n)\)：同时包含 \(s_n\) 与 \(s_{n-1}\) 的可行子集数。

于是有\(f(n)=f(n-1)+V(n)+W(n).\)接下来推导 \(V(n),W(n)\) 的递推。

在 \(V(n)\) 中，\(s_{n-1}\notin S\)。由于 \(\displaystyle{\sum_{i=1}^{n-3}s_i=s_n-3<s_n}\)，若不取 \(s_{n-2}\) 则其余所有可选元素之和仍不足以超过 \(s_n\)，因此必须有 \(s_{n-2}\in S\)。再利用恒等式\(s_n=s_{n-2}+s_{n-3}+s_{n-4}.\)令\(T=S\cap U_{n-3}\subseteq{s_1,\dots,s_{n-3}},\)则多边形条件（最大边 \(s_n\) 小于其余边之和）为\(\displaystyle{\sum_{x\in S\setminus{s_n}}x >s_n.}\)

而\(\displaystyle{\sum_{x\in S\setminus{s_n}}x=s_{n-2}+\sum_{t\in T}t,}\)于是\(\displaystyle{s_{n-2}+\sum_{t\in T}t > s_n\Longleftrightarrow s_{n-2}+\sum_{t\in T}t > s_{n-2}+s_{n-3}+s_{n-4}\Longleftrightarrow\sum_{t\in T}t > s_{n-3}+s_{n-4}.}\)对 \(T\) 是否包含 \(s_{n-3},s_{n-4}\) 分类。

\(s_{n-3},s_{n-4}\in T\)

设\(R=T\cap U_{n-5}\subseteq U_{n-5},\)则\(\displaystyle{\sum_{t\in T}t=s_{n-3}+s_{n-4}+\sum_{r\in R}r.}\)

要使 \(\displaystyle{\sum_{t\in T}t>s_{n-3}+s_{n-4}}\)，等价于 \(\displaystyle{\sum_{r\in R}r>0}\)，即 \(R\ne\varnothing\)。因此该类可选数为\(2^{n-5}-1.\)

\(s_{n-3}\in T,s_{n-4}\notin T\)

写\(T={s_{n-3}}\cup R, R\subseteq U_{n-5}.\)则条件\(\displaystyle{\sum_{t\in T}t>s_{n-3}+s_{n-4}}\)变为\(\displaystyle{s_{n-3}+\sum_{r\in R}r>s_{n-3}+s_{n-4}\Longleftrightarrow\sum_{r\in R}r>s_{n-4}.}\)

而 \(\displaystyle{\sum_{r\in R}r>s_{n-4}}\) 等价于集合 \({s_{n-4}}\cup R\subseteq U_{n-4}\) 以 \(s_{n-4}\) 为最大边满足多边形不等式（此时所有 \(r\in R\) 都小于 \(s_{n-4}\)，若不等式成立则自动至少选到两条“其他边”，从而边数 \(\ge 3\)）。

因此该类数量就是“最大边为 \(s_{n-4}\) 的可行子集数”，即\(V(n-4)+W(n-4).\)

\(s_{n-3}\notin T,s_{n-4}\in T\)

那么写\(T={s_{n-4}}\cup R, R\subseteq U_{n-5}.\)则\(\displaystyle{\sum_{t\in T}t=s_{n-4}+\sum_{r\in R}r.}\)

原条件\(\displaystyle{\sum_{t\in T}t>s_{n-3}+s_{n-4}}\)等价于\(\displaystyle{s_{n-4}+\sum_{r\in R}r>s_{n-3}+s_{n-4}\Longleftrightarrow\sum_{r\in R}r>s_{n-3}.}\)注意到 \(R\subseteq U_{n-5}\subseteq U_{n-4}\)，并且这里我们固定了 \(s_{n-4}\in T\)，等价于固定 \(s_{n-4}\in S\)，同时又固定 \(s_{n-3}\notin T\)（也就是 \(s_{n-3}\notin S\)）。

现在考虑集合\(S'={s_{n-3}}\cup R\subseteq U_{n-3}.\)其最大元素为 \(s_{n-3}\)，并且\(\displaystyle{\sum_{x\in S'\setminus\{s_{n-3}\}}x=\sum_{r\in R}r.}\)因此\(\displaystyle{\sum_{r\in R}r>s_{n-3}\Longleftrightarrow\sum_{x\in S'\setminus\{s_{n-3}\}}x>s_{n-3},}\)也就是说 \(S'\) 以 \(s_{n-3}\) 为最长边满足多边形不等式。

更进一步，由于我们在此分支中要求 \(s_{n-4}\in T\)，而 \(s_{n-4}\notin S'\)（因为 \(S'\) 的元素来自 \({s_{n-3}}\cup U_{n-5}\)），这正对应最大边为 \(s_{n-3}\) 且不包含 \(s_{n-4}\)的计数，即\(V(n-3).\)

三类互斥且覆盖所有可行情况，于是

\[ V(n)=(2^{n-5}-1)+\bigl(V(n-4)+W(n-4)\bigr)+V(n-3). \]

接下来开始推导\(W(n)\)。在 \(W(n)\) 中，\(s_n,s_{n-1}\in S\)。由递推式\(s_n=s_{n-1}+s_{n-3},\)令\(T=S\cap U_{n-2}\subseteq{s_1,\dots,s_{n-2}},\)即 \(T\) 表示除去已固定的 \(s_n,s_{n-1}\) 后，其余被选入的元素集合。多边形条件为\(\displaystyle{\sum_{x\in S\setminus\{s_n\}}x>s_n.}\)

而\(\displaystyle{\sum_{x\in S\setminus{s_n}}x=s_{n-1}+\sum_{t\in T}t,}\)所以\(\displaystyle{s_{n-1}+\sum_{t\in T}t>s_{n-1}+s_{n-3}\Longleftrightarrow\sum_{t\in T}t>s_{n-3}.}\)对 \(T\) 是否包含 \(s_{n-2},s_{n-3}\) 分类：

\(s_{n-2}\in T\)

因为 \(s_{n-2}>s_{n-3}\)，所以\(\displaystyle{\sum_{t\in T}t\ge s_{n-2}>s_{n-3},}\)条件自动满足。此时除去必须包含的 \(s_{n-2}\) 外，其余元素可在 \(U_{n-3}\) 中任取子集，故该类数量为\(2^{n-3}.\)

\(s_{n-2}\notin T,s_{n-3}\in T\)

那么写\(T={s_{n-3}}\cup R, R\subseteq U_{n-4}.\)则\(\displaystyle{\sum_{t\in T}t=s_{n-3}+\sum_{r\in R}r,}\)条件\(\displaystyle{\sum_{t\in T}t>s_{n-3}}\)等价于\(\displaystyle{\sum_{r\in R}r>0,}\)也就是 \(R\ne\varnothing\)。因此该类数量为\(2^{n-4}-1.\)

\(s_{n-2}\notin T,s_{n-3}\notin T\)

此时\(T\subseteq U_{n-4},\)并要求\(\displaystyle{\sum_{t\in T}t>s_{n-3}.}\)将 \(S'={s_{n-3}}\cup T\subseteq U_{n-3}\)，则 \(S'\) 的最大元素为 \(s_{n-3}\)，且\(\displaystyle{\sum_{x\in S'\setminus{s_{n-3}}}x=\sum_{t\in T}t.}\)

因此上述条件等价于 \(S'\) 以 \(s_{n-3}\) 为最大边满足多边形不等式。由于这里并未限制 \(s_{n-4}\) 是否被选入（因为 \(T\subseteq U_{n-4}\) 可以包含也可以不包含 \(s_{n-4}\)），所以该类计数正是“最大边为 \(s_{n-3}\) 的可行子集总数”，即\(V(n-3)+W(n-3).\)

综上得到 \[ W(n)=2^{n-3}+(2^{n-4}-1)+V(n-3)+W(n-3). \]

汇总得到： \[ \begin{aligned} V(n)&=(2^{n-5}-1)+V(n-3)+V(n-4)+W(n-4),\\ W(n)&=(2^{n-3}+2^{n-4}-1)+V(n-3)+W(n-3),\\ f(n)&=f(n-1)+V(n)+W(n). \end{aligned} \]

由于需要 \(f(N)\bmod 10^9\)，必须把含 \(2^k\) 与常数项的“仿射递推”线性化。做法是把状态向量扩充，加入一个“幂变量”和常数 \(1\)。

把递推整体向前挪一位（这样可以避免模 \(10^9\) 下做除法）：

令\(p_n=2^{n-4}\)则\(2^{n-4}=p_n, 2^{n-3}=2p_n, 2^{n-2}=4p_n.\)从而对 \(n\ge 4\) 有

\[ \begin{aligned} V(n+1)&=p_n-1+V(n-2)+V(n-3)+W(n-3),\\ W(n+1)&=6p_n-1+V(n-2)+W(n-2),\\ f(n+1)&=f(n)+V(n+1)+W(n+1),\\ p_{n+1}&=2p_n. \end{aligned} \]

取状态向量 \[ X_n=[f(n),V(n),V(n-1),V(n-2),V(n-3),W(n),W(n-1),W(n-2),W(n-3),p_n,1]^T, \]

则存在 \(11\times 11\) 的矩阵 \(A\)，其中\(A\)是：

\[ A= \begin{pmatrix} 1&0&0&2&1&0&0&1&1&7&-2\\ 0&0&0&1&1&0&0&0&1&1&-1\\ 0&1&0&0&0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0&0&0&0\\ 0&0&0&1&0&0&0&0&0&0&0\\ 0&0&0&1&0&0&0&1&0&6&-1\\ 0&0&0&0&0&1&0&0&0&0&0\\ 0&0&0&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0&0&0&2&0\\ 0&0&0&0&0&0&0&0&0&0&1 \end{pmatrix}. \]

使得\(X_{n+1}= A X_n.\)于是\(X_N= A^{,N-4}X_4,\)而答案就是 \(X_N\) 的第一分量 \(X_{N,0}\)。

代码

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
# Project Euler 382: Generating polygons
# Simplified: treat vectors as (n x 1) matrices; use only mat_mul + mat_pow.

MOD = 10**9
N = 10 ** 18

def mat_mul(A, B, mod=MOD):
    """Multiply A (r x k) and B (k x c) modulo mod."""
    r, k, c = len(A), len(A[0]), len(B[0])
    C = [[0] * c for _ in range(r)]
    for i in range(r):
        Ai = A[i]
        for t in range(k):
            a = Ai[t]
            if a == 0:
                continue
            Bt = B[t]
            for j in range(c):
                C[i][j] = (C[i][j] + a * Bt[j]) % mod
    return C


def mat_pow(A, e, mod=MOD):
    """Fast exponentiation of square matrix A^e modulo mod."""
    n = len(A)
    R = [[0] * n for _ in range(n)]
    for i in range(n):
        R[i][i] = 1
    while e > 0:
        if e & 1:
            R = mat_mul(R, A, mod)
        A = mat_mul(A, A, mod)
        e >>= 1
    return R


def solve(N: int) -> int:
    if N < 3:
        return 0
    if N == 3:
        return 0
    if N == 4:
        return 2

    # X_n = [f(n), V(n), V(n-1), V(n-2), V(n-3), W(n), W(n-1), W(n-2), W(n-3), p_n, 1]^T
    # p_n = 2^(n-4) mod MOD

    m = 11
    A = [[0] * m for _ in range(m)]

    # f(n+1) = f(n) + V(n+1) + W(n+1)
    # V(n+1) = p_n - 1 + V(n-2) + V(n-3) + W(n-3)
    # W(n+1) = 6 p_n - 1 + V(n-2) + W(n-2)
    # => f(n+1) = f(n) + (2*V(n-2) + V(n-3) + W(n-2) + W(n-3) + 7*p_n - 2)
    A[0][0] = 1
    A[0][3] = 2
    A[0][4] = 1
    A[0][7] = 1
    A[0][8] = 1
    A[0][9] = 7
    A[0][10] = MOD - 2

    # V(n+1)
    A[1][3] = 1
    A[1][4] = 1
    A[1][8] = 1
    A[1][9] = 1
    A[1][10] = MOD - 1

    # shift V: V(n) -> V((n+1)-1), etc.
    A[2][1] = 1
    A[3][2] = 1
    A[4][3] = 1

    # W(n+1)
    A[5][3] = 1
    A[5][7] = 1
    A[5][9] = 6
    A[5][10] = MOD - 1

    # shift W
    A[6][5] = 1
    A[7][6] = 1
    A[8][7] = 1

    # p_{n+1} = 2 p_n
    A[9][9] = 2

    # constant 1
    A[10][10] = 1

    # X_4: f(4)=2, V(4)=0, W(4)=2, p_4=1
    X4 = [
        [2],  # f(4)
        [0], [0], [0], [0],  # V(4..1)
        [2], [0], [0], [0],  # W(4..1)
        [1],  # p_4
        [1],  # constant
    ]

    P = mat_pow(A, N - 4, MOD)
    XN = mat_mul(P, X4, MOD)  # (11x11) * (11x1) -> (11x1)
    return XN[0][0] % MOD


print(solve(N))