Project Euler 832

题目

Mex Sequence

In this problem \(\oplus\) is used to represent the bitwise exclusive or of two numbers.

Starting with blank paper repeatedly do the following:

1. Write down the smallest positive integer \(a\) which is currently not on the paper;

2. Find the smallest positive integer \(b\) such that neither \(b\) nor \((a \oplus b)\) is currently on the paper. Then write down both \(b\) and \((a \oplus b)\).

After the first round \(\{1,2,3\}\) will be written on the paper. In the second round \(a=4\) and because \((4 \oplus 5),(4 \oplus 6)\) and \((4 \oplus 7)\) are all already written \(b\) must be \(8\).

After \(n\) rounds there will be \(3n\) numbers on the paper. Their sum is denoted by \(M(n)\).

For example, \(M(10) = 642\) and \(M(1000) = 5432148\).

Find \(M(10^{18})\). Give your answer modulo \(1\,000\,000\,007\).

解决方案

使用以下程序暴力打印前面几行的规律，并得到结果：

st = set()
for i in range(100):
    x = 1
    while x in st:
        x += 1
    y = x + 1
    while y in st or x ^ y in st:
        y += 1
    w = x + y + (x ^ y)
    print(i, x, y, x ^ y, w)
    st |= {x, y, x ^ y}

1 1 2 3 6

2 4 8 12 24
3 5 10 15 30
4 6 11 13 30
5 7 9 14 30

6 16 32 48 96
7 17 34 51 102
8 18 35 49 102
9 19 33 50 102
10 20 40 60 120
11 21 42 63 126
12 22 43 61 126
13 23 41 62 126
14 24 44 52 120
15 25 46 55 126
16 26 47 53 126
17 27 45 54 126
18 28 36 56 120
19 29 38 59 126
20 30 39 57 126
21 31 37 58 126

22 64 128 192 384
23 65 130 195 390
24 66 131 193 390
25 67 129 194 390
26 68 136 204 408
27 69 138 207 414
28 70 139 205 414
29 71 137 206 414
30 72 140 196 408
31 73 142 199 414
32 74 143 197 414
33 75 141 198 414
34 76 132 200 408
35 77 134 203 414
36 78 135 201 414
37 79 133 202 414
-----------------
38 80 160 240 480
39 81 162 243 486
40 82 163 241 486
41 83 161 242 486
42 84 168 252 504
43 85 170 255 510
44 86 171 253 510
45 87 169 254 510
46 88 172 244 504
47 89 174 247 510
48 90 175 245 510
49 91 173 246 510
50 92 164 248 504
51 93 166 251 510
52 94 167 249 510
53 95 165 250 510
-----------------
54 96 176 208 480
55 97 178 211 486
56 98 179 209 486
57 99 177 210 486
58 100 184 220 504
59 101 186 223 510
60 102 187 221 510
61 103 185 222 510
62 104 188 212 504
63 105 190 215 510
64 106 191 213 510
65 107 189 214 510
66 108 180 216 504
67 109 182 219 510
68 110 183 217 510
69 111 181 218 510
-----------------
70 112 144 224 480
71 113 146 227 486
72 114 147 225 486
73 115 145 226 486
74 116 152 236 504
75 117 154 239 510
76 118 155 237 510
77 119 153 238 510
78 120 156 228 504
79 121 158 231 510
80 122 159 229 510
81 123 157 230 510
82 124 148 232 504
83 125 150 235 510
84 126 151 233 510
85 127 149 234 510

86 256 512 768 1536
87 257 514 771 1542
88 258 515 769 1542
89 259 513 770 1542

可以观察到如下规律：

所有的迭代都可以划分成阶段，其中第\(i\)个阶段包含了\(4^i\)轮的迭代。（\(i\)从\(0\)开始）。

其中，第\(i\)个阶段有如下特点：

• \(a\in[4^i,2\cdot 4^{i}),b\in[2\cdot 4^{i},3\cdot 4^i),a\oplus b\in[3\cdot 4^i,4^{i+1})\)，这三个区间刚好连续且不重复地覆盖了\([4^i,4^{i+1})\)这个区间。
• 对于\(i>0\)，第\(i\)个阶段的所有迭代可以均匀分成\(4\)部分，并且这\(4\)部分都可以从第\(i-1\)阶段得出。对于这\(4\)部分\(a,b,a\oplus b\)的关系，可以列出如下表格：

\(\begin{array}{|c|c|c|} \hline a&b&a\oplus b\\ \hline [4^i+0\cdot 4^{i-1},4^i+1\cdot 4^{i-1})&[2\cdot 4^i+0\cdot 4^{i-1},2\cdot4^i+1\cdot 4^{i-1})&[3\cdot 4^i+0\cdot 4^{i-1},3\cdot4^i+1\cdot 4^{i-1})\\ \hline [4^i+1\cdot 4^{i-1},4^i+2\cdot 4^{i-1})&[2\cdot 4^i+2\cdot 4^{i-1},2\cdot4^i+3\cdot 4^{i-1})&[3\cdot 4^i+3\cdot 4^{i-1},3\cdot4^i+4\cdot 4^{i-1})\\ \hline [4^i+2\cdot 4^{i-1},4^i+3\cdot 4^{i-1})&[2\cdot 4^i+3\cdot 4^{i-1},2\cdot4^i+4\cdot 4^{i-1})&[3\cdot 4^i+1\cdot 4^{i-1},3\cdot4^i+2\cdot 4^{i-1})\\ \hline [4^i+3\cdot 4^{i-1},4^i+4\cdot 4^{i-1})&[2\cdot 4^i+1\cdot 4^{i-1},2\cdot4^i+2\cdot 4^{i-1})&[3\cdot 4^i+2\cdot 4^{i-1},3\cdot4^i+3\cdot 4^{i-1})\\ \hline \end{array}\)

并且，每一部分内部的变化量都和第\(i-1\)阶段相同，因此这里考虑使用分治的思想求解\(M(n)\)。

为了方便计算\(M(n)\)，我们将\(n\)进行分块，以\(4^j\)为一块进行整体求和。

令\(f(n,a,b,c,m)\)表示求三元组\((a,b,c)\)以及之后的\(n\)个项目之和，并且视作是以\((a,b,c)\)为开头进行的第\(m\)个阶段。这\(m\)个阶段中，所有数之和为\(\dfrac{(2a+4^m-1)\cdot 4^m}{2}+\dfrac{(2b+4^m-1)\cdot 4^m}{2}+\dfrac{(2c+4^m-1)\cdot 4^m}{2}\)。如果\(n<4^m\)，那么继续将剩下的\(n\)个迭代分成\(4\)部分进行递归求解。

代码

N = 10 ** 18
mod = 10 ** 9 + 7


def f(n, a, b, c, m):
    if n <= 0:
        return 0
    pw = 1 << (2 * m)
    if n > pw:
        return (f(pw, a, b, c, m) + f(n - pw, a << 2, b << 2, c << 2, m + 1)) % mod
    elif n == pw:
        return sum((x + x + pw - 1) * pw // 2 for x in [a, b, c]) %mod
    else:
        pw2 = pw >> 2
        ans = f(min(n, pw2), a, b, c, m - 1)
        ans += f(min(n - 1 * pw2, pw2), a + 1 * pw2, b + 2 * pw2, c + 3 * pw2, m - 1)
        ans += f(min(n - 2 * pw2, pw2), a + 2 * pw2, b + 3 * pw2, c + 1 * pw2, m - 1)
        ans += f(min(n - 3 * pw2, pw2), a + 3 * pw2, b + 1 * pw2, c + 2 * pw2, m - 1)
        return ans % mod


ans = f(N, 1, 2, 3, 0)
print(ans)