Project Euler 728

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

Project Euler 728

题目

Circle of Coins

Consider \(n\) coins arranged in a circle where each coin shows heads or tails. A move consists of turning over \(k\) consecutive coins: tail-head or head-tail. Using a sequence of these moves the objective is to get all the coins showing heads.

Consider the example, shown below, where \(n=8\) and \(k=3\) and the initial state is one coin showing tails (black). The example shows a solution for this state.

For given values of \(n\) and \(k\) not all states are solvable. Let \(F(n,k)\) be the number of states that are solvable. You are given that \(F(3,2) = 4\), \(F(8,3) = 256\) and \(F(9,3) = 128\).

Further define:

\[\displaystyle S(N) = \sum_{n=1}^N\sum_{k=1}^n F(n,k)\]

You are also given that \(S(3) = 22\), \(S(10) = 10444\) and \(S(10^3) \equiv 853837042 \pmod{1\,000\,000\,007}\)

Find \(S(10^7)\). Give your answer modulo \(1\,000\,000\,007\)

解决方案

我们把硬币状态编码为二进制向量：第 \(i\) 枚硬币 \(c_i\in\{0,1\}\)（\(i=0,1,\dots,n-1\)），规定 \(0=\text{H},1=\text{T}\)。那么整体状态为 \(c=(c_0,c_1,\dots,c_{n-1})\in\mathbb F_2^n\)。

一次操作等价于把某个长度为 \(k\) 的连续区间全部异或上 \(1\)。也就是说，允许加上模 \(2\) 的窗口向量：\(v_s = e_s+e_{s+1}+\cdots+e_{(s+k-1)\% n}\)。其中 \(s\in\{0,1,\dots,n-1\}\)，\(e_s\) 是 \(\mathbb F_2^n\) 里的标准基向量：也就是只有第 \(s\) 个分量为 \(1\)，其余全是 \(0\)。

也就是说，所有可达状态都可以表示为若干个 \(v_s\) 的异或组合。于是：从全零（全 H）出发能够到达的状态数，恰好就是可解状态数 \(F(n,k)\)。即\(F(n,k)\) 等于所有这些 \(v_s\) 通过异或叠加所能生成的不同状态总数。

不过，若直接去统计这些组合究竟能产生多少种不同结果，计算与分析都会变得不够直观。为此，下面引入一种更强的结构化不变量，用它来刻画一次操作对状态的影响，从而更容易得到 \(F(n,k)\) 的闭式表达。

定义相邻差分（模 \(2\)）：\(d_i = c_i\oplus c_{i+1} (0\le i\le n-1,\ c_n=c_0)\)，得到一个长度 \(n\) 的向量 \(d=(d_0,\dots,d_{n-1})\)。

翻转一段连续 \(k\) 枚硬币时：区间内部相邻两个硬币 同时翻转，因此它们的异或不变；只有区间两端边界的相邻关系改变。

设翻转区间从位置 \(s\) 开始，即翻转集合为\(\{s,s+1,\dots,s+k-1\}\pmod n\)，则一次操作只会翻转差分序列中恰好两个位置：\(d_{(s-1)\%k}\)和\(d_{(s+k-1)\%k}\)二者在环上相距 \(k\)。

换句话说，在差分空间中，每一步是：选择某个 \(i\)，同时翻转 \(d_i\) 与 \(d_{(i+k)\%n}\)。

记\(g=\gcd(n,k)\)，考虑在集合 \({0,1,\dots,n-1}\) 上做加 \(k\)的跳跃：\(i\mapsto i+k\pmod n\)。这个置换会把点分成 \(g\) 个互不相交的循环，每个循环长度为\(m=\dfrac{n}{g}\)而差分操作翻转 \(d_i\) 与 \(d_{(i+k)\%n}\)只会发生在同一个循环内，因此不同循环互相独立。

在一个长度为 \(m\) 的循环上，每次操作翻转两个点。于是该循环内的差分位总异或（即所有位之和 \(\bmod 2\)）保持不变，因为每次翻转两位，奇偶不变。因此该循环内可达差分序列构成一个维数为 \(m-1\) 的子空间，大小 \(2^{m-1}\)。

合并 \(g\) 个循环，因此总可达差分序列数为：\((2^{m-1})^g = 2^{g(m-1)} = 2^{n-g}\)

于是得到结论：从全零差分序列出发，可达的差分序列恰好有 \(2^{n-g}\) 种。

不过，从差分回到原序列：还差一个“整环翻转”的判定。给定差分序列 \(d\)，原硬币序列 \(c\) 是否唯一？

由递推\(c_{i+1}=c_i\oplus d_i (0\le i\le n-1,\ c_n=c_0)\)可知，一旦确定 \(c_0\)，其余全部确定；而 \(c_0\) 有两种选择，所以：每个差分序列对应恰好 2 个硬币状态：\(c\) 与 \(c\oplus \mathbf 1\)，其中 \(\mathbf 1=(1,1,\dots,1)\)。

但可解性并不一定对这两者同时成立，关键在于：是否能通过操作把某个状态变成其整体翻转 \(c\oplus\mathbf 1\)，这等价于：是否能生成向量 \(\mathbf 1\)。这取决于 \(\dfrac{k}{g}\) 的奇偶性。

我们并把硬币按步长 \(g\) 分成 \(g\) 组：第 \(i\) 组为位置 \(i, i+g, i+2g, \dots, i+(m-1)g\)，其中 \(m=n/g\)，且 \(i=0,1,\dots,g-1\)。

定义每组的异或：\(C_i = c_i\oplus c_{i+g}\oplus\cdots\oplus c_{i+(m-1)g}, i=0,1,\dots,g-1\)。这是一个长度 \(g\) 的“组特征向量” \(C=(C_0,\dots,C_{g-1})\)。

一次操作翻转连续 \(k\) 枚硬币，而 \(k=gl\)，其中\(l=\dfrac{k}{g}\)。翻转区间在每一组中覆盖的硬币数完全相同，恰好是 \(l\) 个。因此：每个 \(C_i\) 会被翻转 \(l\) 次（XOR 上 \(1\) 的次数为 \(l\)）。所以若 \(l\) 为偶数，则 \(C_i\) 不变；若 \(l\) 为奇数，则 \(C_i\) 取反。

于是得到两种情形。

情形 A：\(\dfrac{k}{g}\) 为偶数

此时每次操作都保持 \(C=(C_0,\dots,C_{g-1})\) 不变，因此 \(C\) 是严格不变量。总状态数为 \(2^n\)；且不同 \(C\) 的取值有 \(2^g\) 种；并且每个 \(C\) 类大小相同，为 \(2^{n-g}\)。

全 H 状态对应 \(C=(0,\dots,0)\)，因此可解状态数为该类大小：\(F(n,k)=2^{n-g}\)。

情形 B：\(\dfrac{k}{g}\) 为奇数

此时每次操作都会把所有 \(C_i\) 同时取反，因此 \(C\) 本身不再不变，但它在“整体取反”等价下保持不变：\(C \sim C\oplus(1,1,\dots,1)\)

因此“类数”减半：\(C\) 的取值共有 \(2^g\) 种，但每个类把 \(C\) 与其整体取反合并，因此类数为 \(2^{g-1}\)。于是可解状态数为：\(F(n,k)=\dfrac{2^n}{2^{g-1}}=2^{n-g+1}\).

综上所述，综合两种情况，\(F(n,k)\)可以写成：

\[ F(n,k)=2^{,n-\gcd(n,k)+\left(\frac{k}{\gcd(n,k)}\bmod 2\right)} \]

从定义：\(\displaystyle{S(N)=\sum_{n=1}^N\sum_{k=1}^n F(n,k)}\)，对每一对 \((n,k)\)，令\(g=\gcd(n,k), n=mg, k=lg, \gcd(m,l)=1\)，其中 \(1\le l\le m\)，且 \(1\le g\le \left\lfloor \dfrac{N}{m}\right\rfloor\)。代入闭式，有：\(F(mg,lg)=2^{mg-g}\times((l\bmod 2) + 1)\)。

因此

\[S(N)=\sum_{m=1}^{N}\sum_{g=1}^{\lfloor N/m\rfloor}2^{g(m-1)}\left(\sum_{\substack{1\le l\le m\\ \gcd(l,m)=1}} ((l\bmod 2) + 1)\right)\]

记 \(\displaystyle{W(m)=\sum_{\substack{1\le l\le m\\ \gcd(l,m)=1}} ((l\bmod 2) + 1)}\)。设与 \(m\) 互质的数中，奇数个数为 \(A(m)\)，偶数个数为 \(B(m)\)，则 \(A(m)+B(m)=\varphi(m)\)，并且\(W(m)=2A(m)+B(m)=\varphi(m)+A(m)\)。接下来我们分类讨论计算\(W(m)\)：

  • \(m\)为偶数。若 \(\gcd(l,m)=1\)，则 \(l\) 必为奇数（否则与 \(2\) 有公因子）。因此\(A(m)=\varphi(m), B(m)=0\Rightarrow W(m)=2\varphi(m)\)。
  • \(m>1\) 为奇数。 则互质的数在奇偶上严格对半： 因为映射 \(l\mapsto m-l\) 保持互质性，并且 \(m\) 奇数使得 \(l\) 与 \(m-l\) 奇偶相反，构成一一配对。故\(A(m)=B(m)=\dfrac{\varphi(m)}{2}\Rightarrow W(m)=\varphi(m)+\dfrac{\varphi(m)}{2}=\dfrac{3}{2}\varphi(m)\)
  • \(m=1\)。只有 \(l=1\)（奇），所以$W(1)=2。

综上所述，有\(W\)的讨论结果：

\[ W(m)= \begin{cases} 2,&m=1\\ 2\varphi(m),&m\ge 2\land m\equiv 0\pmod2\\ \dfrac{3}{2}\varphi(m),&m\ge 3\land m\equiv 1\pmod2 \end{cases} \]

令\(T=\left\lfloor\dfrac{N}{m}\right\rfloor\)，则计算\(\displaystyle{\sum_{g=1}^{T}2^{g(m-1)}}\)时，

  • 当 \(m=1\) 时，指数为 \(0\)，所以该和为 \(T=N\)。
  • 当 \(m\ge 2\) 时，这是公比为 \(2^{m-1}\) 的几何级数：\(\displaystyle{\sum_{g=1}^{T}2^{g(m-1)}=\sum_{g=1}^{T}\left(2^{m-1}\right)^g\frac{2^{m-1}\left(2^{(m-1)T}-1\right)}{2^{m-1}-1}}\)。

因此最终可以将\(\displaystyle{S(N)=\sum_{m=1}^N W(m)\cdot \sum_{g=1}^{\lfloor N/m\rfloor}2^{g(m-1)}}\)写成显式形式：

  • \(m=1\) 项：\(W(1)=2\)，内和为 \(N\)，贡献 \(2N\)。
  • \(m\ge 2\) 项：\(\displaystyle{S(N)=2N+\sum_{m=2}^{N}W(m)\cdot\frac{2^{m-1}\left(2^{(m-1)\lfloor N/m\rfloor}-1\right)}{2^{m-1}-1}}\)。

我们可以进一步简化成：

\[ S(N)=2N+\sum_{m=2}^{N} \varphi(m)\cdot C(m)\cdot \frac{2^{m-1}\left(2^{(m-1)\lfloor N/m\rfloor}-1\right)}{2^{m-1}-1} \]

其中\(C(m)=\begin{cases}2,&m\equiv 0\pmod{2}\\\dfrac{3}{2},&m\equiv 1\pmod{2}\\\end{cases}\)。

代码

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
N = 10 ** 7
MOD = 1_000_000_007

def totient_sieve_linear(n: int):
    """phi[0..n] in O(n) time, using linear sieve."""
    phi = [0 for _ in range(n+1)]
    phi[1] = 1
    is_comp = bytearray(n + 1)
    primes = []

    for i in range(2, n + 1):
        if not is_comp[i]:
            primes.append(i)
            phi[i] = i - 1
        for p in primes:
            x = i * p
            if x > n:
                break
            is_comp[x] = 1
            if i % p == 0:
                phi[x] = phi[i] * p
                break
            else:
                phi[x] = phi[i] * (p - 1)
    return phi

def solve(N: int) -> int:
    inv2 = (MOD + 1) // 2

    # 1) phi
    phi = totient_sieve_linear(N)

    # 2) pow2[i] = 2^i mod MOD, i=0..N
    pow2 = [0 for _ in range(N+1)]
    pow2[0] = 1
    for i in range(1, N + 1):
        pow2[i] = (pow2[i - 1] * 2) % MOD

    # 3) batch inverse for denom(m) = 2^(m-1) - 1, m>=2
    #    store prefix products in pref, and later overwrite pref[m] as inv_denom[m]
    pref = [0 for _ in range(N+1)]
    pref[1] = 1
    zero = bytearray(N + 1)  # zero[m]=1 iff denom(m)==0

    acc = 1
    for m in range(2, N + 1):
        d = (pow2[m - 1] - 1) % MOD
        if d == 0:
            zero[m] = 1
            pref[m] = acc
        else:
            acc = (acc * d) % MOD
            pref[m] = acc

    inv_acc = pow(pref[N], MOD - 2, MOD)
    for m in range(N, 1, -1):
        if zero[m]:
            pref[m] = 0
        else:
            pref[m] = (inv_acc * pref[m - 1]) % MOD
            d = (pow2[m - 1] - 1) % MOD
            inv_acc = (inv_acc * d) % MOD

    # m = 1 special: W(1)=2, G(1)=N
    res = (2 * (N % MOD)) % MOD

    for m in range(2, N + 1):
        T = N // m  # floor(N/m)

        # W(m)
        ph = phi[m]
        if (m & 1) == 0:
            W = 2 * ph % MOD
        else:
            W = 3 * ph * inv2 % MOD

        # G(m)
        # ratio a = 2^(m-1)
        if zero[m]:
            G = T % MOD
        else:
            a = pow2[m - 1]
            exp = (m - 1) * T
            term = (pow2[exp] - 1) % MOD
            G = (a * term) % MOD
            G = (G * pref[m]) % MOD  # multiply inv(a-1)

        res += W * G
        res %= MOD

    return res

print(solve(N))