Project Euler 843

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

Project Euler 843

题目

Periodic Circles

This problem involves an iterative procedure that begins with a circle of \(n\ge 3\) integers. At each step every number is simultaneously replaced with the absolute difference of its two neighbours.

For any initial values, the procedure eventually becomes periodic.

Let \(S(N)\) be the sum of all possible periods for \(3\le n \leq N\). For example, \(S(6) = 6\), because the possible periods for \(3\le n \leq 6\) are \(1, 2, 3\). Specifically, \(n=3\) and \(n=4\) can each have period \(1\) only, while \(n=5\) can have period \(1\) or \(3\), and \(n=6\) can have period \(1\) or \(2\).

You are also given \(S(30) = 20381\).

Find \(S(100)\).

解决方案

令\(N=100\)。记一步变换为 \(T\)。对任意位置 \(i\)，有\(|a_{i-1}-a_{i+1}|\le \max(a_{i-1},a_{i+1})\le \max_j \{a_j\}.\)因此全局最大值 \(M_t=\max_i a_i^{(t)}\) 单调不增且为非负整数，故最终稳定为某个 \(M\)。状态空间离散，因而序列最终进入一个循环（可能先经过前导段）。

进入循环段后，最大值恒为 \(M\)。如果某一步出现了介于 \(0\) 和 \(M\) 之间的值，它会通过差分在邻域扩散，迫使某处最大值下降，和最大值已稳定为矛盾。更形式化地：在循环段内若存在 \(0<a_i<M\)，则其两侧在下一步产生的差分也会落入 \((0,M)\) 或 \([0,M)\)，从而使严格小于 \(M\)的连续区间长度增长；在环上增长到覆盖全体后就会导致所有值都 \(<M\)，矛盾。因此循环段里只能出现 \(0\) 或 \(M\)。

当 \(n\) 为偶数时还有一个结构性现象：在顶点集 \(V=\{0,1,\dots,n-1\}\) 上定义无向图 \(G_2=(V,E)\)，其中\(E=\{(i,(i+2)\bmod n) : i\in V\},\)也就是说用步长 \(2\)（模 \(n\)）把位置相连；当 \(n\) 为偶数时该图分成两个长度为 \(n/2\) 的环（偶下标与奇下标各自成环），它们在反向约束下可对应不同的最大值 \(M_{\text{even}},M_{\text{odd}}\)。于是循环段里偶位只可能是 \({0,M_{\text{even}}}\)，奇位只可能是 \({0,M_{\text{odd}}}\)。

关键结论： 循环节长度只取决于哪里为 0、哪里为最大值的 0/1 形状；把所有最大值同时除掉（偶/奇两环分别除以各自最大值）不会改变 period。因此讨论所有可能 period 时，足够只研究初值在 \({0,1}\) 上的情形。

现在令 \(x_i\in\{0,1\}\)。此时\(|x_{i-1}-x_{i+1}| = x_{i-1}\oplus x_{i+1},\)而在 \(\mathbb F_2\) 上 \(\oplus\) 就是加法，于是更新变为线性变换\(x_i' \equiv x_{i-1}+x_{i+1}\pmod 2.\) 写成向量 \(x' = M_n x\)，其中 \(M_n\) 是对称循环矩阵：每行在相邻两列为 1，其余为 0。比如 \(n=5\)： \[ M_5=\begin{pmatrix} 0&1&0&0&1\\ 1&0&1&0&0\\ 0&1&0&1&0\\ 0&0&1&0&1\\ 1&0&0&1&0 \end{pmatrix}. \]

因此“某个 period \(p\) 可实现” \(\Longleftrightarrow\) 存在某个非零 \(x\in\mathbb F_2^n\) 的轨道循环节为 \(p\)。

我们先确定一个所有 period 的公共上界 \(L(n)\)。

当\(n\)是奇数时， 考虑把 \(M_n\) 放到一个能包含 \(n\) 次单位根的扩域中。循环矩阵可被傅里叶基对角化：对每个满足 \(\zeta^n=1\) 的 \(\zeta\)，存在特征值\(\lambda(\zeta)=\zeta+\zeta^{-1}.\)在特征方向上，迭代就是乘 \(\lambda(\zeta)\)。

由于我们在特征 \(2\) 的域上工作，Frobenius 映射 \(x\mapsto x^2\) 是环同态，因此对任意 \(k\ge 0\) 有 \((a+b)^{2^k}=a^{2^k}+b^{2^k}\)。于是\(\lambda(\zeta)^{2^k}=(\zeta+\zeta^{-1})^{2^k}=\zeta^{2^k}+\zeta^{-2^k}.\)若 \(2^k\equiv 1\pmod n\)，则 \(\zeta^{2^k}=\zeta\)，从而 \(\lambda(\zeta)^{2^k}=\lambda(\zeta)\)； 若 \(2^k\equiv -1\pmod n\)，则 \(\zeta^{2^k}=\zeta^{-1}\)，也同样得到 \(\lambda(\zeta)^{2^k}=\lambda(\zeta)\)。

因此对所有非零特征值都有 \(\lambda^{2^k-1}=1\)，从而线性变换在“非零特征空间”上的所有轨道循环节都整除 \(2^k-1\)。令 \[ k(n)=\min\{k\ge 1:2^k\equiv \pm 1\pmod n\}, L(n)=2^{k(n)}-1. \] 则所有可能 period 必为 \(L(n)\) 的因子。

当 \(n\) 是偶数时， 令 \(n=2m\)。把序列按下标奇偶拆成两条长度为 \(m\) 的子序列\(e_k=a_{2k\bmod m}, o_k=a_{2k+1\bmod m}\)，并定义长度为 \(m\) 的相邻差分Ducci算子\((\Delta b)_k=\lvert b_k-b_{(k+1)\bmod m}\rvert.\)由更新式 \(a_i'=\lvert a_{i-1}-a_{i+1}\rvert\) 直接计算可得\(o_k'=\lvert a_{2k}-a_{2k+2}\rvert=\lvert e_k-e_{k+1}\rvert=(\Delta e)_k,\)以及\(e_k'=\lvert a_{2k-1}-a_{2k+1}\rvert=\lvert o_{k-1}-o_k\rvert=(\Delta o)_{k-1}.\)

再做一步更新：\(o_k''=(\Delta e')_k, e_k''=(\Delta o')_{k-1}.\)把上一步的表达式代回去（注意一切下标均按 \(m\) 取模）：\(o_k''=(\Delta e')_k=\lvert e_k'-e_{k+1}'\rvert=|( \Delta o)_{k-1}-(\Delta o)_k|,e_k''=(\Delta o')_{k-1}=\lvert o_{k-1}'-o_k'\rvert=|(\Delta e)_{k-1}-(\Delta e)_k|.\)在 \(\{0,1\}\)（等价于 \(\mathbb F_2\)）情形下有 \(\lvert x-y\rvert=x\oplus y\)，于是差分再差分可化简为隔一个取邻居：\((\Delta(\Delta b))_k=(b_k\oplus b_{k+1})\oplus(b_{k+1}\oplus b_{k+2})=b_k\oplus b_{k+2}.\)因此两步迭代满足

\[ o_k''=o_{k-1}\oplus o_{k+1}, e_k''=e_{k-1}\oplus e_{k+1}, \]

这说明：做两步之后，偶位子序列 \((e_k)\) 自身按同样规则演化，奇位子序列 \((o_k)\) 也自身按同样规则演化，且二者互不影响。于是长度为 \(2m\) 的系统在“两步算子”意义下分解成两个长度为 \(m\) 的同类系统，从而所有可能周期都整除 \(2L(m)\)。

另一方面，取初态使得奇位全为 \(0\)，偶位取一个在长度 \(m\) 上具有最大周期 \(L(m)\) 的状态 \((e_k)\)。则每一步后奇位仍由偶位确定、偶位仍由奇位确定，因此轨道必然存在；并且由上面的两步关系知偶位子序列每两步按长度 \(m\) 的规则前进一次，因此其最小回到初态所需步数恰为 \(2L(m)\)。综上得到严格递推：

  • 若 \(n\) 是 2 的幂，则会在某步“熄灭”为全 0，只有 period \(1\)；
  • 否则 \(L(2m)=2L(m)\)。

于是对 \(n\le N\) 可直接用递推求得所有 \(L(n)\)，并保证每个真实 period 都整除 \(L(n)\)。

对某个候选周期 \(d\mid L(n)\)，令\(V_d=\{x\in\mathbb F_2^n: M_n^d x=x\}=\ker(M_n^d+I).\)其中\(V_d\) 是“周期整除 \(d\)”的全部状态集合。如果对两个因子 \(d_1<d_2\) 有 \(V_{d_1}=V_{d_2}\)，那么任何周期整除 \(d_2\) 的状态其实周期也整除 \(d_1\)，因此 \(d_2\) 不可能作为最小周期出现。反之，只要 \(V_d\) 严格新增了约束结构（与更小因子不同），就会对应出现一个新的可实现最小周期。

矩阵 \(M_n^d+I\) 是循环矩阵，其行空间对循环位移不变。把它做高斯消元后，行空间的一个自然基可以选为某一行模式的所有循环位移，其余行全为 0。于是 \(V_d\)（也就是核）等价于满足这一行模式在所有位移下给出的线性约束。

因此，用消元后的最后一条非零行作为一个签名模式即可区分不同的 \(V_d\)。若两个 \(d\) 得到相同签名，则它们的 \(V_d\) 相同，较大的那个 \(d\) 不会产生新周期。

因此，最终我们对每个 \(n=3,4,\dots,N\)，计算上界 \(L(n)\)，并枚举 \(L(n)\) 的所有因子 \(d\)（从小到大），计算签名 \((n,d)\)：快速幂求 \(M_n^d\)，再做 \(M_n^d+I\)；用位运算高斯消元得到规范形；并取最后一条非零行作为签名。若该签名首次出现，则把 \(d\) 记为该 \(n\) 的可实现最小周期。

最后把所有 \(n\) 的周期集合取并集求和即得 \(S(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
from functools import lru_cache

from tools import divisors
N = 100

def get_step_matrix(n: int):
    M = [(1 << (n - 2)) | 1]
    row = (1 << (n - 1)) | (1 << (n - 3))
    M.append(row)
    for _ in range(2, n - 1):
        row >>= 1
        M.append(row)
    M.append((1 << (n - 1)) | 2)
    return M


def matrix_mult(A, B):
    n = len(A)
    msb = 1 << (n - 1)
    C = [0] * n
    row0 = A[0]
    new_row = 0
    for j in range(n):
        new_row <<= 1
        t = row0 & B[j]
        if t and (t.bit_count() & 1):
            new_row |= 1
    C[0] = new_row
    for i in range(1, n):
        new_row = (new_row >> 1) | (msb * (new_row & 1))
        C[i] = new_row
    return C


def matrix_pow(M, exp: int):
    n = len(M)
    A = [1 << (n - i - 1) for i in range(n)]
    while exp:
        if exp & 1:
            A = matrix_mult(M, A)
        exp >>= 1
        if exp:
            M = matrix_mult(M, M)
    return A


def matrix_reduce(M):
    M = M.copy()
    n = len(M)
    bit = 1 << (n - 1)
    i = 0
    while bit and i < n:
        M.sort(reverse=True)
        if M[i] == 0:
            break
        row = M[i]
        for j in range(i + 1, n):
            if (M[j] & bit) == 0:
                break
            M[j] ^= row
        i += 1
        bit >>= 1
    return M


def get_divisor_pattern(n: int, d: int, stepM):
    Mpow = matrix_pow(stepM, d)
    MpowI = Mpow.copy()
    bit = 1 << n
    for i in range(n):
        bit >>= 1
        MpowI[i] ^= bit
    Mred = matrix_reduce(MpowI)
    for i in range(n - 1, -1, -1):
        if Mred[i] != 0:
            return Mred[i]
    return 0


def find_max_periods(N: int):
    res = [0] * (N + 1)
    for n in range(3, N + 1):
        if n % 2 == 0:
            if n & (n - 1) == 0:
                res[n] = 1
            else:
                res[n] = 2 * res[n // 2]
        else:
            pow2 = 1
            k = 0
            while True:
                k += 1
                pow2 = (pow2 * 2) % n
                if pow2 == 1 or pow2 == n - 1:
                    break
            res[n] = (1 << k) - 1
    return res


@lru_cache(None)
def all_divisors(m: int):
    return tuple(divisors(m))


def periods_for_n(n: int, max_period: int):
    stepM = get_step_matrix(n)
    seen = set()
    per = set()
    for d in all_divisors(max_period):
        pat = get_divisor_pattern(n, d, stepM)
        if pat not in seen:
            seen.add(pat)
            per.add(d)
    return per


def S(N: int) -> int:
    max_periods = find_max_periods(N)
    allp = set()
    for n in range(3, N + 1):
        allp |= periods_for_n(n, max_periods[n])
    return sum(allp)


print(S(N))