Project Euler 707

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

Project Euler 707

题目

Lights Out

Consider a \(w\times h\) grid. A cell is either ON or OFF. When a cell is selected, that cell and all cells connected to that cell by an edge are toggled on-off, off-on. See the diagram for the \(3\) cases of selecting a corner cell, an edge cell or central cell in a grid that has all cells on (white).

The goal is to get every cell to be off simultaneously. This is not possible for all starting states. A state is solvable if, by a process of selecting cells, the goal can be achieved.

Let \(F(w,h)\) be the number of solvable states for a \(w\times h\) grid. You are given \(F(1,2)=2\), \(F(3,3) = 512\), \(F(4,4) = 4096\) and \(F(7,11) \equiv 270016253 \pmod{1\,000\,000\,007}\).

Let \(f_1=f_2 = 1\) and \(f_n=f_{n-1}+f_{n-2}, n \ge 3\) be the Fibonacci sequence and define \[ S(w,n) = \sum_{k=1}^n F(w,f_k)\] You are given \(S(3,3) = 32\), \(S(4,5) = 1052960\) and \(S(5,7) \equiv 346547294 \pmod{1\,000\,000\,007}\).

Find \(S(199,199)\). Give your answer modulo \(1\,000\,000\,007\).

解决方案

令\(W=N=199.\)把每个格子的 ON/OFF 视为 \(\mathbb F_2\) 上的 0/1。对一个 \(w\times h\) 网格共有 \(wh\) 个格子，按键方案（每格按或不按）也对应一个 \(\mathbb F_2^{wh}\) 向量 \(x\)。

令 \(A\) 为 \(wh\times wh\) 的矩阵，表示一次按键集合对整盘状态的翻转线性算子：\(A\) 的第 \(j\) 列表示按第 \(j\) 个格子会翻转哪些格子（自己 + 四邻接）。初始状态为 \(c\in\mathbb F_2^{wh}\)，按键后变为 \(c + Ax\)。论文指出，把 \(c\) 变成全 0 等价于解线性方程 \(Ax=c\)。

因此，可解初始状态集合正是线性映射 \(A\) 的像空间 \(\mathrm{Im}(A)\)。在 \(\mathbb F_2\) 上，\(F(w,h)=|\mathrm{Im}(A)| = 2^{\mathrm{rank}(A)}.\)又因为 \(\mathrm{rank}(A)+\mathrm{nullity}(A)=wh\)，得\(\log_2 F(w,h)= wh - \mathrm{nullity}(A).\)

所以关键变成：如何求网格算子 \(A\) 的零空间维数 \(\mathrm{nullity}(A)\)。

论文对网格 \(G_{n,m}\)（\(n\) 列 \(m\) 行）引入 Fibonacci 多项式 \(f_k(x)\)（系数在 \(\mathbb F_2\) 上）：\(g_0(x)=1, g_1(x)=x, g_k(x)=x g_{k-1}(x)+g_{k-2}(x)(k\ge 2).\)并在定理11b给出：\(A^+_{n,m}\) 的零空间维数等于 \(\gcd(g_n(x+1), g_m(x))\) 的次数。

因此对本题的 \(w\times h\)（只要把列/行对换也无影响），可取\(\mathrm{nullity}(A) = \deg\bigl(\gcd(g_w(x+1), g_h(x))\bigr),\)从而

\[ F(w,h)=2^{wh-\deg(\gcd(g_w(x+1), g_h(x)))}. \]

我们把\(S(W,N)=\displaystyle \sum_{k=1}^N F(W,f_k)\)的第 \(k\) 项写成\(\displaystyle F(W,f_k)=2^{W\cdot f_k - d_k},d_k=\deg\bigl(\gcd(P(x), g_{f_k}(x))\bigr),\)其中固定\(P(x)=g_{W}(x+1).\)

利用欧几里得算法性质：\(\gcd(P,Q)=\gcd(P,Q\bmod P).\)所以求 \(d_k\) 时，只需知道\(r_k(x)= g_{f_k}(x)\bmod P(x),\)然后\(d_k=\deg(\gcd(P,r_k)).\)这样所有多项式运算都可在商环 \(\mathbb F_2[x]/(P)\) 中进行，次数永远 \(<\deg P\le W\)，不会随着 \(g_k\) 巨大而爆炸。

令\(M=\begin{pmatrix}x&1\\1&0\end{pmatrix}.\)可验证（对 \(n\ge 1\) 归纳）：\(M^n=\begin{pmatrix}g_n & g_{n-1}\\g_{n-1} & g_{n-2}\end{pmatrix},\)其中取 \(g_{-1}=0\).于是由 \(M^{m+n}=M^mM^n\) 的 \((1,1)\) 元在\(\mathbb F_2\)上得到\(g_{m+n}=g_m g_n + g_{m-1}g_{n-1}.\)

同理，\((1,1)\) 元取 \(M^{m+n+1}=M^{m+1}M^n\) 给出\(g_{m+n+1}=g_{m+1}g_n + g_m g_{n-1}.\)此外由递推式 \(g_{t+1}=x g_t+g_{t-1}\) 得\(g_{t-1}=g_{t+1}+x g_t\)（因为在 \(\mathbb F_2\) 上减法=加法）。

维护三元组：\(T_k=\bigl(g_{f_k-1},g_{f_k},g_{f_k+1}\bigr)\bmod P.\)已知 \(T_1=T_2=(g_0,g_1,g_2)=(1,x,x^2+1)\)。

设 \(m=g_{k-1},n=g_{k-2}\)，则\(g_{m+n}=g_m g_n + g_{m-1}g_{n-1},g_{m+n+1}=g_{m+1}g_n + g_m g_{n-1},g_{m+n-1}=g_{m+n+1}+x g_{m+n}.\)右侧都只用到 \(T_{k-1}\) 与 \(T_{k-2}\) 中的分量，加上一次乘以 \(x\)（移位）即可完成更新。

这样就能在次数始终 \(<W\) 的约束下，依次得到所有 \(r_k=f_{g_k}\bmod P\)，再用多项式 gcd 得到 \(d_k\)。最终回代到\(F\)中计算出\(S\)即可。

代码

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
from tools import phi
W = 199
N = 199
MOD = 1000000007
PHI = phi(MOD)

def deg(p):
    return p.bit_length() - 1

def polymod(p, mod):
    if p == 0:
        return 0
    md = deg(mod)
    while p and deg(p) >= md:
        p ^= mod << (deg(p) - md)
    return p

def polygcd(a, b):
    while b:
        a, b = b, polymod(a, b)
    return a

def polymul(a, b):
    r = 0
    while b:
        if b & 1:
            r ^= a
        b >>= 1
        a <<= 1
    return r

def polymul_mod(a, b, mod):
    return polymod(polymul(a, b), mod)

def fibpoly_xplus1(n):
    x1 = 3
    if n == 0:
        return 1
    if n == 1:
        return x1
    a, b = 1, x1
    for _ in range(2, n + 1):
        a, b = b, ((b << 1) ^ b) ^ a
    return b

def solve(w,n):
    P = fibpoly_xplus1(w)
    x = 2

    fibm = [0] * (n + 1)
    fibm[1] = 1
    fibm[2] = 1
    for k in range(3, n + 1):
        fibm[k] = (fibm[k - 1] + fibm[k - 2]) % PHI

    def xmul(p):
        return polymod(p << 1, P)

    f0 = 1
    f1 = x
    f2 = (x << 1) ^ 1

    t1 = [f0, f1, f2]
    t2 = [f0, f1, f2]

    def term(k, t):
        r = t[1]
        g = polygcd(P, r)
        d = deg(g)
        e = (w * fibm[k] - d) % PHI
        return pow(2, e, MOD)

    ans = (term(1, t1) + term(2, t2)) % MOD

    for k in range(3, n + 1):
        a0, a1, a2 = t2
        b0, b1, b2 = t1
        new1 = polymul_mod(b1, a1, P) ^ polymul_mod(b0, a0, P)
        new2 = polymul_mod(b2, a1, P) ^ polymul_mod(b1, a0, P)
        new0 = new2 ^ xmul(new1)
        t = [new0, new1, new2]
        ans = (ans + term(k, t)) % MOD
        t2, t1 = t1, t

    return ans

print(solve(W, N))