Project Euler 671

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

Project Euler 671

题目

Colouring a Loop

A certain type of flexible tile comes in three different sizes - \(1\times1, 1\times2,\) and \(1\times3\) - and in \(k\) different colours. There is an unlimited number of tiles available in each combination of size and colour.

These are used to tile a closed loop of width \(2\) and length (circumference) \(n\), where \(n\) is a positive integer, subject to the following conditions:

  • The loop must be fully covered by non-overlapping tiles.

  • It is not permitted for four tiles to have their corners meeting at a single point.

  • Adjacent tiles must be of different colours.

For example, the following is an acceptable tiling of a \(2\times 23\) loop with \(k=4\) (blue, green, red and yellow):

but the following is not an acceptable tiling, because it violates the “no four corners meeting at a point” rule:

Let \(F_k(n)\) be the number of ways the \(2\times n\) loop can be tiled subject to these rules when \(k\) colours are available. (Not all \(k\) colours have to be used.) Where reflecting horizontally or vertically would give a different tiling, these tilings are to be counted separately.

For example, \(F_4(3) = 104\), \(F_5(7) = 3327300\), and \(F_6(101)\equiv 75309980 \pmod{1\,000\,004\,321}\). Find \(F_{10}(10\,004\,003\,002\,001) \bmod 1\,000\,004\,321\).

解决方案

令\(n=10004003002001\)。本题和题目 670 背景类似，这里只保留背景相似的分析结论。考虑从左到右逐列观察，令第 \(t\) 列上、下格颜色分别为 \(a_t,b_t\)。则有：同一行中，相邻列若颜色相同，则必属于同一块水平砖延伸；因此每一行的同色连续段长度最多为 \(3\)（对应水平砖长度最多为 \(3\)）。\(a_t=b_t\) 当且仅当第 \(t\) 列是一块竖直 \(2\times1\) 砖。

类似的，四角相遇只可能发生在某条列分界线 \(t\mid t+1\) 与上下两行分界线的交点处。因此，在 \(a_t\ne b_t\land a_{t+1}\ne b_{t+1}\) 时，必须满足 \(a_{t+1}=a_t\) 或 \(b_{t+1}=b_t\)（至少一行继续同色段）。

可见，颜色名称本身不重要，真正影响计数的是：颜色之间的相等/不等关系；新颜色需要避开多少种已出现颜色。对所有颜色做任意置换会把合法铺法双射到合法铺法，因此可以用规范代表色压缩状态。这一思想在题目 670 已经用过；本题因为是环形计数且要精确处理颜色选择数，需要更细的颜色标签信息。

在计算 \(G_k(n)\)（把循环平移视为不同，固定切口）时，采用如下规范化：用标签 \(1,2\) 表示两种被区分的代表色；用标签 \(0\) 表示其余 \(k-2\) 种颜色中的某一种。

沿环固定一条切口，按列推进。切口处只需记录当前被切到的砖形态。

若切口切到一块竖直 \(2\times1\) 砖，记为\(V(t), t\in\{0,1,2\},\)其中 \(t\) 为该竖砖颜色标签。共有 \(3\) 个竖直状态。

若切口处上下两格分别属于两块水平砖，则记为\(H(c,d,x,y),\)其中\(c,d\)表示切口右侧这一列上、下格的颜色标签；\(x,y\in\{0,1,2\}\)表示该格在所属水平砖中的位置减一（即第 \(x+1\) / \(y+1\) 格）。

因为上下两格相邻（共享边）必须异色，所以 \((c,d)\) 只可能属于以下 \(7\) 种规范化组合：\((0,0),(1,0),(0,1),(2,0),(0,2),(1,2),(2,1).\)其中 \((0,0)\) 表示上下都是其它色，且这两个其它色实际必须不同。因此水平状态数为：\(7\times 3\times 3=63.\)总状态数为：\(3+63=66.\)

令 \(M\) 为 \(66\times 66\) 转移矩阵。一次转移表示向右推进一列。若当前状态是 \(H(c,d,x,y)\)，则下一列上、下两行各自要么继续当前水平砖，要么开启新砖，因此上行位置参数 \(X\)、下行位置参数 \(Y\) 满足：\(X\in\{0,x+1\}, Y\in\{0,y+1\},\)并且必须满足长度上限（即 \(X,Y\le 2\)）。再由前文的等价条件，在水平/水平 \(\to\) 水平/水平的转移中，不允许上下两行同时在边界处换砖，因此禁止\(X=0\land Y=0.\)

令矩阵元素 \(M[s][t]\) 表示从状态 \(s\) 转到状态 \(t\) 的真实颜色方案数。注意规范化标签 \(0\) 实际对应 \(k-2\) 种颜色，且新放置砖时要避开与其共享边相邻砖的颜色。于是权重只取决于：新颜色是否为标签 \(1/2\)（唯一选择）；若为标签 \(0\)，需要避开多少种已知颜色（可选数为 \(k-2-m\)）。

尤其在 \((0,0)\) 类型中，上下两格都属于其它色且必须彼此不同，因此会出现乘积项，例如：不额外禁用时，权值为\((k-2)(k-3)\)；如果需要额外避开某个其它色时，权值为 \((k-3)(k-4)\)。

我们先定义 \(G_k(n)\) 为固定切口，循环平移视为不同的合法铺法数。题目要求的 \(F_k(n)\) 是把旋转视为相同后的结果。对某个循环群 \(C_n\) 使用 Burnside 引理：

\[ F_k(n)=\frac{1}{n}\sum_{r=0}^{n-1}\mathrm{Fix}(r), \]

其中 \(\mathrm{Fix}(r)\) 表示旋转 \(r\) 列后不变的铺法数。

对按列状态模型，若旋转 \(r\) 固定某铺法，则其列模式周期为 \(\gcd(n,r)\)，因此可写成约数和形式：\(\displaystyle F_k(n)=\frac{1}{n}\sum_{d\mid n}\varphi(d)G_k\left(\frac{n}{d}\right).\)

本题中\(n\)是素数。于是除恒等旋转外，其余旋转若固定铺法，则该铺法必须每列相同。这在本题约束下不可能成立（会与相邻异色或四角约束冲突），因此非平凡旋转固定点数为 \(0\)，从而有：

\[F_k(n)=\dfrac{G_k(n)}{n}.\]

固定切口后，一个合法环铺法对应一条长度为 \(n\) 的状态闭环，因此对任意状态 \(s\)，回到自身的方案数是 \((M^n)[s][s]\)。要恢复真实计数 \(G_k(n)\)，只需对规范化起点加上颜色重数。

如果起点为竖直砖，取起点状态为 \(V(1)\)。代表切口处是一块竖直砖，其颜色被命名为标签 \(1\)。真实颜色有 \(k\) 种选择，因此贡献为\(k\cdot (M^n)[V(1),V(1)].\)

如果起点为水平/水平，取起点为 \(H(1,2,x,y)\)，表示切口处上下是两块不同水平砖，并把它们颜色命名为标签 \(1,2\)。真实有序选色数为 \(k(k-1)\)，且切口可能落在砖内任意位置，因此要对 \(x,y\in{0,1,2}\) 求和，贡献为\(\displaystyle k(k-1)\sum_{x=0}^{2}\sum_{y=0}^{2}(M^n)[H(1,2,x,y),H(1,2,x,y)].\)两类贡献相加即为 \(G_k(n)\)。

于是最终答案为：

\[ F_k(n)=\dfrac{1}{n} \left( k\cdot (M^n)[V(1),V(1)] + k(k-1)\sum_{x=0}^{2}\sum_{y=0}^{2}(M^n)[H(1,2,x,y),H(1,2,x,y)] \right) \]

代码

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
107
108
109
110
from itertools import product

K = 10
N = 10004003002001
MOD = 1000004321

def solve(k, n, mod=MOD):
    def state_to_int(b, t):
        if b:
            return t
        c, d, x, y = t
        mp = {(0, 0): 0, (1, 0): 1, (0, 1): 2, (2, 0): 3, (0, 2): 4, (1, 2): 5, (2, 1): 6}
        j = mp[(c, d)]
        return 3 + 9 * j + 3 * x + y

    def int_to_state(i):
        if i < 3:
            return True, i
        j, m = divmod(i - 3, 9)
        mp = [(0, 0), (1, 0), (0, 1), (2, 0), (0, 2), (1, 2), (2, 1)]
        c, d = mp[j]
        x, y = divmod(m, 3)
        return False, (c, d, x, y)

    def build_matrix():
        size = 66
        M = [[0] * size for _ in range(size)]
        for i in range(size):
            b, t = int_to_state(i)
            if b:
                for c in (0, 1, 2):
                    if t != 0 and t == c:
                        continue
                    z = (k - 2 - int(t == 0)) if c == 0 else 1
                    j = state_to_int(True, c)
                    M[i][j] = (M[i][j] + z) % mod
                for c, d in product((0, 1, 2), (0, 1, 2)):
                    if c != 0 and c == d:
                        continue
                    if t != 0 and (c == t or d == t):
                        continue
                    if (c, d) == (0, 0):
                        z = (k - 3) * (k - 4) if t == 0 else (k - 2) * (k - 3)
                    elif c == 0 or d == 0:
                        z = (k - 3) if t == 0 else (k - 2)
                    else:
                        z = 1
                    j = state_to_int(False, (c, d, 0, 0))
                    M[i][j] = (M[i][j] + z) % mod
            else:
                c, d, x, y = t
                for C in (0, 1, 2):
                    if C == 0:
                        z = k - 2 - int(c == 0) - int(d == 0)
                    else:
                        z = int(C != c and C != d)
                    j = state_to_int(True, C)
                    M[i][j] = (M[i][j] + z) % mod
                for C, D, X, Y in product((0, 1, 2), (0, 1, 2), (0, x + 1), (0, y + 1)):
                    if X == 0 and Y == 0:
                        continue
                    if X > 2 or Y > 2:
                        continue
                    if C != 0 and C == D:
                        continue
                    if (X > 0 and C != c) or (Y > 0 and D != d):
                        continue
                    if (X == 0 and C == c and C != 0) or (Y == 0 and D == d and D != 0):
                        continue
                    if (X == 0 and C == 0) or (Y == 0 and D == 0):
                        z = k - 2 - int(c == 0) - int(d == 0)
                    else:
                        z = 1
                    j = state_to_int(False, (C, D, X, Y))
                    M[i][j] = (M[i][j] + z) % mod
        return M

    def mat_mul(A, B, mod):
        return [[sum(A[i][k] * B[k][j] for k in range(len(A[0]))) % mod for j in range(len(B[0]))] for i in
                range(len(A))]

    def mat_pow(M, e):
        n0 = len(M)
        R = [[0] * n0 for _ in range(n0)]
        for i in range(n0):
            R[i][i] = 1
        A = M
        while e:
            if e & 1:
                R = mat_mul(R, A, mod)
            e >>= 1
            if e:
                A = mat_mul(A, A, mod)
        return R

    M = build_matrix()
    Mn = mat_pow(M, n)

    res = 0
    i = state_to_int(True, 1)
    res = (res + k * Mn[i][i]) % mod
    for x in range(3):
        for y in range(3):
            i = state_to_int(False, (1, 2, x, y))
            res = (res + k * (k - 1) * Mn[i][i]) % mod

    inv_n = pow(n, mod - 2, mod)
    return res * inv_n % mod

print(solve(K, N))
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
from itertools import product
import sympy

K = 10
N = 10004003002001
MOD = 1000004321

def factorization(n):
    ls = list(sympy.factorint(n).items())
    ls.sort()
    return ls

def solve_any_n(k, n, mod=MOD):
    def state_to_int(b, t):
        if b:
            return t
        c, d, x, y = t
        mp = {(0, 0): 0, (1, 0): 1, (0, 1): 2, (2, 0): 3, (0, 2): 4, (1, 2): 5, (2, 1): 6}
        j = mp[(c, d)]
        return 3 + 9 * j + 3 * x + y

    def int_to_state(i):
        if i < 3:
            return True, i
        j, m = divmod(i - 3, 9)
        mp = [(0, 0), (1, 0), (0, 1), (2, 0), (0, 2), (1, 2), (2, 1)]
        c, d = mp[j]
        x, y = divmod(m, 3)
        return False, (c, d, x, y)

    def build_matrix(mod_now):
        size = 66
        M = [[0] * size for _ in range(size)]
        for i in range(size):
            b, t = int_to_state(i)
            if b:
                for c in (0, 1, 2):
                    if t != 0 and t == c:
                        continue
                    z = (k - 2 - int(t == 0)) if c == 0 else 1
                    j = state_to_int(True, c)
                    M[i][j] = (M[i][j] + z) % mod_now
                for c, d in product((0, 1, 2), (0, 1, 2)):
                    if c != 0 and c == d:
                        continue
                    if t != 0 and (c == t or d == t):
                        continue
                    if (c, d) == (0, 0):
                        z = (k - 3) * (k - 4) if t == 0 else (k - 2) * (k - 3)
                    elif c == 0 or d == 0:
                        z = (k - 3) if t == 0 else (k - 2)
                    else:
                        z = 1
                    j = state_to_int(False, (c, d, 0, 0))
                    M[i][j] = (M[i][j] + z) % mod_now
            else:
                c, d, x, y = t
                for C in (0, 1, 2):
                    if C == 0:
                        z = k - 2 - int(c == 0) - int(d == 0)
                    else:
                        z = int(C != c and C != d)
                    j = state_to_int(True, C)
                    M[i][j] = (M[i][j] + z) % mod_now
                for C, D, X, Y in product((0, 1, 2), (0, 1, 2), (0, x + 1), (0, y + 1)):
                    if X == 0 and Y == 0:
                        continue
                    if X > 2 or Y > 2:
                        continue
                    if C != 0 and C == D:
                        continue
                    if (X > 0 and C != c) or (Y > 0 and D != d):
                        continue
                    if (X == 0 and C == c and C != 0) or (Y == 0 and D == d and D != 0):
                        continue
                    if (X == 0 and C == 0) or (Y == 0 and D == 0):
                        z = k - 2 - int(c == 0) - int(d == 0)
                    else:
                        z = 1
                    j = state_to_int(False, (C, D, X, Y))
                    M[i][j] = (M[i][j] + z) % mod_now
        return M

    def mat_mul(A, B, mod):
        return [[sum(A[i][k] * B[k][j] for k in range(len(A[0]))) % mod for j in range(len(B[0]))] for i in
                range(len(A))]

    def precompute_powers(M, emax, mod_now):
        pw = [M]
        while (1 << len(pw)) <= emax:
            pw.append(mat_mul(pw[-1], pw[-1], mod_now))
        return pw

    def mat_pow_from_powers(pw, e, mod_now):
        n0 = len(pw[0])
        R = [[0] * n0 for _ in range(n0)]
        for i in range(n0):
            R[i][i] = 1
        b = 0
        while e:
            if e & 1:
                R = mat_mul(R, pw[b], mod_now)
            e >>= 1
            b += 1
        return R

    def rooted_count_from_power(Mn, mod_now):
        res = 0
        i = state_to_int(True, 1)
        res = (res + k * Mn[i][i]) % mod_now
        for x in range(3):
            for y in range(3):
                i = state_to_int(False, (1, 2, x, y))
                res = (res + k * (k - 1) * Mn[i][i]) % mod_now
        return res

    def gen_divisors_with_phi(fac):
        ans = []
        def dfs(idx, d, ph):
            if idx == len(fac):
                ans.append((d, ph))
                return
            p, e = fac[idx]
            dfs(idx + 1, d, ph)
            dd = d
            pp = ph
            for t in range(1, e + 1):
                dd *= p
                if t == 1:
                    pp2 = pp * (p - 1)
                else:
                    pp2 = pp * (p - 1) * (p ** (t - 1))
                dfs(idx + 1, dd, pp2)
        dfs(0, 1, 1)
        return ans

    fac = factorization(n)
    divs_phi = gen_divisors_with_phi(fac)
    need_ms = sorted({n // d for d, _ in divs_phi})
    mod_for_dp = mod if n % mod != 0 else mod * n

    M = build_matrix(mod_for_dp)
    pw = precompute_powers(M, max(need_ms), mod_for_dp)

    g_cache = {}
    for m in need_ms:
        Mn = mat_pow_from_powers(pw, m, mod_for_dp)
        g_cache[m] = rooted_count_from_power(Mn, mod_for_dp)

    num = 0
    for d, ph in divs_phi:
        m = n // d
        num = (num + (ph % mod_for_dp) * g_cache[m]) % mod_for_dp

    if n % mod != 0:
        return num * pow(n % mod, mod - 2, mod) % mod

    return (num // n) % mod

print(solve_any_n(K, N))