Project Euler 361

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

Project Euler 361

题目

Subsequence of Thue-Morse sequence

The Thue-Morse sequence \(\{T_n\}\) is a binary sequence satisfying:

  • \(T_0 = 0\)
  • \(T_{2n} = T_n\)
  • \(T_{2n+1} = 1 - T_n\)

The first several terms of \({T_n}\) are given as follows:

\(01101001{\color{red}{10010}}1101001011001101001\dots\).

We define \(\{A_n\}\) as the sorted sequence of integers such that the binary expression of each element appears as a subsequence in \(\{T_n\}\).

For example, the decimal number \(18\) is expressed as \(10010\) in binary. \(10010\) appears in \(\{T_n\}\) (\(T_8\) to \(T_{12}\)), so \(18\) is an element of \(\{A_n\}\).

The decimal number \(14\) is expressed as \(1110\) in binary. \(1110\) never appears in \(\{T_n\}\), so \(14\) is not an element of \(\{A_n\}\).

The first several terms of \(A_n\) are given as follows:

\(n\) \(0\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\) \(7\) \(8\) \(9\) \(10\) \(11\) \(12\) \(\dots\)
\(A_n\) \(0\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\) \(9\) \(10\) \(11\) \(12\) \(13\) \(18\) \(\dots\)

We can also verify that \(A_{100} = 3251\) and \(A_{1000} = 80852364498\).

Find the last \(9\) digits of \(\sum_{k=1}^{18}A_{10^k}\).

解决方案

Thue–Morse 序列 \((T_n)\) 定义为\(T_n = s_2(n)\bmod 2,\)其中 \(s_2(n)\) 是 \(n\) 的二进制 1 的个数。题目中 \(A\) 是所有满足「其二进制表示作为 连续子串 出现在 \((T_n)\) 中」的非负整数按从小到大排序得到的序列，且 \(A_0=0\)。

对 \(x>0\)，其标准二进制表示首位一定是 1，因此：

  • 对每个长度 \(k\ge 1\)，落在区间 \([2^{k-1},2^k)\) 的数恰好对应 所有以 1 开头、长度为 \(k\) 的二进制串。
  • 同一长度 \(k\) 内，数值大小与二进制串的字典序一致。
  • 不同长度按长度分块：所有 \(k\) 位数都小于任何 \((k+1)\) 位数。

所以问题转化为：对每个 \(k\)，令 \(L(k)\) 为 Thue–Morse 中出现过的、以 1 开头的长度 \(k\) 子串的不同个数。 令\(\displaystyle{S(k)=\sum_{i=1}^k L(i)}\)，则 \(A_n(n\ge1)\) 在长度分块里：找最小 \(k\) 使 \(S(k)\ge n\)，并在长度 \(k\) 的那一块里取第\(o=n-S(k-1)-1\)个（从 0 开始）字典序子串，再把该子串当作二进制数。

接下来要解决三件事：

  1. 快速算 \(L(k)\) 与 \(S(k)\)，以便从 \(n\) 找到长度 \(k\) 与 \(o\)。
  2. 给定 \((k,o)\)，快速找到该子串在 Thue–Morse 中的一个起点位置 \(pos\)。
  3. 给定 \((pos,k)\)，快速求该长度 \(k\) 子串对应的二进制整数 \(\bmod 10^9\)，而不显式生成 \(k\) 位。

可见，Thue–Morse 也是如下 2-一致代换的不动点：\(\mu(0)=01, \mu(1)=10,\)注意这里的\(\mu\)只是一个函数名，用来表示一个替换规则。并把它作用在任意二进制串上：把串中每个字符分别替换后再拼接（例如 \(\mu(010)=\mu(0)\mu(1)\mu(0)=011001\)）。从初始串 \(0\) 出发不断重复这一替换，会得到一串不断增长的前缀：

\[ 0, 01, 0110, 01101001, \dots \]

这些前缀在长度上依次翻倍，且每一步都是上一步的延长；因此它们确定了一条唯一的无限二进制序列。该无限序列的第 \(n\) 位恰好等于 Thue–Morse 序列的 \(T_n\)，也就是说 Thue–Morse 就是由这条替换规则生成的。

考虑长度 \(k>4\) 的一个出现过的子串 \(w=T[pos:pos+k)\)。由于代换长度为 2，子串 \(w\) 必然来自两种“对齐方式”之一：

  • 偶对齐（pos 为偶数）：\(pos=2q\)。此时 \(w\) 落在 \(\mu(T[q:])\) 的块边界上。把 \(w\) 每两位分组 \((w_0w_1)(w_2w_3)\cdots\)，就能通过规则\(\mu\)唯一还原出一个较短串 \(u\)，使得 \(w\) 是 \(\mu(u)\) 的前 \(k\) 位。此时 \(|u|=\lceil k/2\rceil\)，且首位保持不变（\(\mu(1)\) 以 1 开头）。
  • 奇对齐（pos 为奇数）：\(pos=2q+1\)。此时 \(w\) 从某个 \(\mu(\cdot)\) 块的第二位开始，相当于对 \(\mu(v)\) 做“左移一位再截取”。此时 \(|v|=\lfloor k/2\rfloor+1\)，并且首位翻转：因为 \(\mu(0)=01\) 的第二位是 1，而 \(\mu(1)=10\) 的第二位是 0。

并且当 \(k>4\) 时，这种可对齐方式是唯一的（因为通过为\(k\)的串生成新的长度为\(2k,2k-1,2k-2\)这4类串是唯一的，短的交替串是唯一会产生二义性的情况，长度 \(\le4\) 可以直接硬处理/预处理）。

于是对“以 1 开头的长度 \(k\) 子串集合”，它由两部分构成：

  • 来自偶对齐：对应“以 1 开头、长度 \(\lceil k/2\rceil\)”的子串集合，数量 \(L(\lceil k/2\rceil)\)；
  • 来自奇对齐：对应“以 0 开头、长度 \(\lfloor k/2\rfloor+1\)”的子串集合。由于 Thue–Morse 的因子集合对按位取反封闭（出现一个子串则其按位取反也出现），所以“以 0 开头”的数量与“以 1 开头”相同，数量为 \(L(\lfloor k/2\rfloor+1)\)。

从而得到（对 \(k>4\)）：

  • 若 \(k=2n:L(2n)=L(n)+L(n+1)\)
  • 若 \(k=2n+1:L(2n+1)=2L(n+1)\)

初值：\(L(1)=1,L(2)=2,L(3)=3,L(4)=5.\)

因此，我们要用二分从 \(n\) 找 \(k\)。因此也需要 \(S(k)\) 的 \(O(\log k)\) 算法。对 \(k\ge 7\)，用上面的 \(L\) 递推把 \(S(2n)\)、\(S(2n+1)\) 写成 \(S(n),S(n+1)\) 的线性组合（推导时需要把小的特殊项常数单独扣出来，最终得到对 \(n\ge3\) 的统一形式）：

\[ \begin{aligned} S(2n) &= 3S(n) + S(n+1) - 4,\\ S(2n+1) &= S(n) + 3S(n+1) - 4, \end{aligned} \]

其中，\(n\ge 3\)。并且由\(L(1\dots6)\)给出基值求得：\(S(0)=0,S(1)=1,S(2)=3,S(3)=6,S(4)=11,S(5)=17,S(6)=25.\)

这样 \(S(k)\) 也能在 \(O(\log k)\) 时间得到，于是可以对 \(k\) 做指数扩张 + 二分，找到最小 \(k\) 使 \(S(k)\ge n\)，并得到 \(o\)。

定义 \(p(k,o,\texttt{bit})\)：表示 以 \(\texttt{bit}\) 开头 的所有长度 \(k\) 因子按字典序排序后，第 \(o\) 个因子的一个出现位置起点。

对 \(k\le4\) 直接预处理表即可（长度很小，枚举 Thue–Morse 前几十位就足够覆盖所有长度 \(\le4\) 因子）。

对 \(k>4\)，令\(a=\left\lceil\dfrac{k}{2}\right\rceil,b=\left\lfloor\dfrac{k}{2}\right\rfloor+1.\)

那么偶对齐那部分有 \(L(a)\) 个，奇对齐那部分有 \(L(b)\) 个。关键点是 字典序里两部分的先后顺序取决于首位 \(\texttt{bit}\)：

  • 首位为 1 时：偶对齐那组排在前（直观上它来自 \(\mu(1)=10\) 的块边界结构），奇对齐排在后；
  • 首位为 0 时：顺序相反（可由按位取反导致字典序整体反向来推出）。

同时，位置在代换层级上满足：

  • 偶对齐：子串起点是偶数，映射到上一层起点乘 2；
  • 奇对齐：子串起点是奇数，映射到上一层起点乘 2 再加 1，且首位翻转。

因此递推为（\(k>4\)）：

  • 若 \(\texttt{bit}=1\)：

    • 若 \(o < L(a):p(k,o,1)=2\cdot p(a,o,1)\)
    • 否则：\(p(k,o,1)=2\cdot p(b,o-L(a),0)+1\)

  • 若 \(\texttt{bit}=0\)：

    • 若 \(o < L(b):p(k,o,0)=2\cdot p(b,o,1)+1\)
    • 否则：\(p(k,o,0)=2\cdot p(a,o-L(b),0)\)

最终我们需要的是 \(pos = p(k,o,1)\)。

接下来需要计算二进制子串 \(T[pos:pos+k)\) 所表示的整数。

对每个 \(h\ge 0\)，考虑 Thue–Morse 序列中从位置 \(0\) 开始、长度为 \(2^h\) 的前缀块，记为\(Z_h, |Z_h|=2^h.\)由于 Thue–Morse 具有“自相似 + 取反”的结构，长度 \(2^h\) 的对齐块只可能是两种之一：

  • \(Z_h\) 本身；
  • \(\overline{Z_h}\)：对 \(Z_h\) 逐位取反得到的块。

并且有如下对齐性质：若 \(pos\) 是 \(2^h\) 的倍数，则

\[ T[pos:pos+2^h)= \begin{cases} Z_h,& T_{pos/2^h}=0,\\ \overline{Z_h},& T_{pos/2^h}=1. \end{cases} \]

设模数 \(M=10^9\)。预处理两类量：

  • 幂\(P_h = 2^{2^h}\bmod M;\)
  • 块值\(V_h(0)\equiv \text{value}(Z_h)\pmod M, V_h(1)\equiv \text{value}(\overline{Z_h})\pmod M,\) 其中 \(\text{value}(\cdot)\) 表示把一个二进制串当作二进制整数来解释。

块的构造满足\(Z_{h+1}=Z_h\Vert \overline{Z_h},\)

因此得到递推： \[ \begin{aligned} P_{h+1} &\equiv P_h^2 \pmod M,\\ V_{h+1}(0) &\equiv V_h(0)\cdot P_h + V_h(1)\pmod M,\\ V_{h+1}(1) &\equiv V_h(1)\cdot P_h + V_h(0)\pmod M. \end{aligned} \]

接下来，将区间 \([pos,pos+k)\) 做标准的 dyadic 分解：从左到右，每次取一个起点对齐到 \(2^h\)（即当前 \(cur\) 可被 \(2^h\) 整除），且\(2^h\) 不超过剩余长度的最大块。

若取到的块长度为 \(2^h\)，起点为 \(cur\)，则它对应的块类型由\(b=T_{cur/2^h}\)决定，于是该块的模值就是 \(V_h(b)\)。

用“二进制拼接”的方式滚动更新答案：\(ans \leftarrow ans\cdot P_h + V_h(b)\pmod M.\)

dyadic 分解的块数为 \(O(\log k)\)，因此整个取值过程也是 \(O(\log k)\)，无需显式生成长度为 \(k\) 的子串。

代码

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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
#!/usr/bin/env python3
from functools import lru_cache

MOD = 10**9


# -------------------------
# 1) Thue–Morse bit
# -------------------------
def tm_bit(i: int) -> int:
    # T_i = popcount(i) mod 2
    return i.bit_count() & 1


# -------------------------
# 2) L(k): number of distinct length-k factors starting with '1'
# -------------------------
@lru_cache(None)
def L(k: int) -> int:
    if k == 1:
        return 1
    if k == 2:
        return 2
    if k == 3:
        return 3
    if k == 4:
        return 5
    if k & 1 == 0:
        n = k // 2
        return L(n) + L(n + 1)
    else:
        n = (k - 1) // 2
        return 2 * L(n + 1)


# -------------------------
# 3) S(k) = sum_{i<=k} L(i), fast via halving recurrences
# base values from L(1..6): 1,2,3,5,6,8  => prefix: 1,3,6,11,17,25
# -------------------------
_S_BASE = [0, 1, 3, 6, 11, 17, 25]


@lru_cache(None)
def S(k: int) -> int:
    if k <= 6:
        return _S_BASE[k]
    if k & 1 == 0:
        n = k // 2
        return 3 * S(n) + S(n + 1) - 4
    else:
        n = k // 2
        return S(n) + 3 * S(n + 1) - 4


def find_length_and_offset(n: int) -> tuple[int, int]:
    """
    For n>=1, find minimal k with S(k)>=n, then offset = n - S(k-1) - 1.
    """
    l, r = 1, 1
    while S(r) < n:
        r <<= 1
    while l < r:
        mid = (l + r) >> 1
        if S(mid) >= n:
            r = mid
        else:
            l = mid + 1
    k = l
    offset = n - S(k - 1) - 1
    return k, offset


# -------------------------
# 4) small base table for k<=4: positions of lex-sorted distinct factors
#    for both starting bits 0/1.
# -------------------------
def build_base_table():
    N = 256  # plenty for k<=4
    t = [tm_bit(i) for i in range(N)]
    base = {}
    for k in range(1, 5):
        for bit in (0, 1):
            words = set()
            for i in range(N - k + 1):
                if t[i] == bit:
                    words.add(tuple(t[i : i + k]))
            lex = sorted(words)
            pos_list = []
            for w in lex:
                for i in range(N - k + 1):
                    if tuple(t[i : i + k]) == w:
                        pos_list.append(i)
                        break
            base[(k, bit)] = pos_list
    return base


BASE = build_base_table()


# -------------------------
# 5) Pos(k, offset, bit): unranking -> starting index in Thue–Morse
# -------------------------
@lru_cache(None)
def p(k: int, offset: int, bit: int) -> int:
    if k <= 4:
        return BASE[(k, bit)][offset]

    a = (k + 1) // 2          # ceil(k/2)
    b = k // 2 + 1            # floor(k/2)+1
    cnt_even = L(a)
    cnt_odd = L(b)

    if bit == 1:
        if offset < cnt_even:
            return 2 * p(a, offset, 1)
        else:
            return 2 * p(b, offset - cnt_even, 0) + 1
    else:
        if offset < cnt_odd:
            return 2 * p(b, offset, 1) + 1
        else:
            return 2 * p(a, offset - cnt_odd, 0)


# -------------------------
# 6) Precompute mu^h blocks' numeric values modulo MOD
#    pow2[h] = 2^(2^h) mod MOD
#    val[h][b] = value of mu^h(b) as binary integer mod MOD
# -------------------------
MAXH = 40  # enough (we only need up to ~30 for lengths <= 1e9)
pow2 = [0] * (MAXH + 1)
val = [[0, 0] for _ in range(MAXH + 1)]

pow2[0] = 2 % MOD          # 2^(1)
val[0][0] = 0
val[0][1] = 1

for h in range(1, MAXH + 1):
    pow2[h] = (pow2[h - 1] * pow2[h - 1]) % MOD
    # mu^{h}(0) = mu^{h-1}(0) || mu^{h-1}(1)
    val[h][0] = (val[h - 1][0] * pow2[h - 1] + val[h - 1][1]) % MOD
    # mu^{h}(1) = mu^{h-1}(1) || mu^{h-1}(0)
    val[h][1] = (val[h - 1][1] * pow2[h - 1] + val[h - 1][0]) % MOD


def substring_value(pos: int, length: int) -> int:
    """
    Compute integer value of T[pos:pos+length) in binary, modulo MOD,
    using dyadic aligned blocks.
    """
    end = pos + length
    cur = pos
    acc = 0
    while cur < end:
        rem = end - cur
        h_rem = rem.bit_length() - 1               # max h with 2^h <= rem
        lb = cur & -cur                            # largest power-of-two dividing cur
        h_div = lb.bit_length() - 1
        h = h_rem if h_rem < h_div else h_div      # aligned block exponent
        b = tm_bit(cur >> h)                       # which mu^h(block) to use
        acc = (acc * pow2[h] + val[h][b]) % MOD
        cur += 1 << h
    return acc


def A_mod(n: int) -> int:
    if n == 0:
        return 0
    k, off = find_length_and_offset(n)
    pos = p(k, off, 1)
    return substring_value(pos, k)


if __name__ == "__main__":
    ans = 0
    n = 1
    for _ in range(18):
        n *= 10
        ans = (ans + A_mod(n)) % MOD

    # print last 9 digits (zero-padded if needed)
    print(f"{ans:09d}")