Project Euler 588

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

Project Euler 588

题目

Quintinomial coefficients

The coefficients in the expansion of \((x+1)^k\) are called binomial coefficients.

Analoguously the coefficients in the expansion of \((x^4+x^3+x^2+x+1)^k\) are called quintinomial coefficients.(quintus= Latin for fifth).

Consider the expansion of \((x^4+x^3+x^2+x+1)^3\):

\(x^{12}+3x^{11}+6x^{10}+10x^9+15x^8+18x^7+19x^6+18x^5+15x^4+10x^3+6x^2+3x+1\)

As we can see \(7\) out of the \(13\) quintinomial coefficients for \(k=3\) are odd.

Let \(Q(k)\) be the number of odd coefficients in the expansion of \((x^4+x^3+x^2+x+1)^k\).

So \(Q(3)=7\).

You are given \(Q(10)=17\) and \(Q(100)=35\).

Find \(\sum_{k=1}^{18}Q(10^k)\).

解决方案

我们只关心每个系数的奇偶性，因此可以在模 \(2\) 的意义下工作：把所有系数都放到 \(\mathbb F_2=\{0,1\}\) 中。

在 \(\mathbb F_2[x]\) 里：

  • 系数为奇数 \(\Longleftrightarrow\) 模 \(2\) 后系数为 \(1\)
  • 系数为偶数 \(\Longleftrightarrow\) 模 \(2\) 后系数为 \(0\)

所以：\(Q(k)\) 就是多项式 \(P(x)^k\) 在 \(\mathbb F_2[x]\) 中的非零项数量。

最终，我们可以把多项式用一个集合\(S\)表示：\(\displaystyle{F(x)=\sum_{e\in S}x^e,}\)此时多项式\(F\)和\(G\)的加法就是集合的对称差：\(S_{F+G}=S_F\triangle S_G,\)乘法对应指数相加并统计奇偶次出现。

我们可以知道， \(\mathbb F_2[x]\) 满足 Frobenius 自同态性质：\(\displaystyle{\left(\sum a_i x^i\right)^2=\sum a_i x^{2i}.}\)也就是说：平方不会产生交叉项，只会把每一项的指数乘 \(2\)。 因此：\(P(x)^{2^t} = P(x^{2^t})=1+x^{2^t}+x^{2\cdot 2^t}+x^{3\cdot 2^t}+x^{4\cdot 2^t}.\)

这让我们可以用二进制快速幂计算 \(P(x)^n\)，同时只维护指数集合。

设 \(n\) 的二进制写作（从高到低）：\(n=\overline{B_1\texttt{00}B_2},\)其中 \(B_1,B_2\) 都是二进制串（允许含 \(\texttt{0}\) 与 \(\texttt{1}\)），并且中间至少出现一次分隔符 \(\texttt{00}\)。

等价于把 \(n\) 写成\(n = A\cdot 2^{t}+B,\)其中 \(t\ge 2\)（因为中间至少有两个 \(0\)），并且 \(B<2^{t-2}\)。

对应到多项式：\(P(x)^n = P(x)^{A\cdot 2^{t}}\cdot P(x)^B.\)可见有两点：

  1. 由平方翻倍性质，有\(P(x)^{A\cdot 2^{t}} = P(x^{2t})^A\)所以这个因子的所有指数都是 \(2^t\) 的倍数。
  2. 另一个因子 \(P(x)^B\) 的最高次数不超过 \(4B\)，而\(4B < 2^t.\)

于是当我们做乘法时，任意一项指数要么是\(2^t\)的倍数，要么是\(<2^t\)的数。这样的表示是唯一的，不会出现不同的配对得到同一个指数的碰撞，也就不会产生模 \(2\) 下的抵消。

因此得到结论：如果二进制中两段之间至少隔着 00，则\(Q(n)=Q(A)\cdot Q(B).\)

把这个过程反复应用：把 \(n\) 的二进制表示（忽略末尾的 \(0\)）按 00 切开：\(n =\overline{m_1\texttt{00}m_2\texttt{00}\cdots\texttt{00}m_r,}\)那么\(\displaystyle{Q(n)=\prod_{i=1}^r Q(m_i).}\)

这些块 \(m_i\) 的特点是：以 \(1\) 开头结尾，且不含子串 00（它们正对应 OEIS A247647 那类二进制串）。

对本题的 \(n=10^k\)，二进制里 \(0\) 非常多，所以会被切成很多短块；最大的块也很小，因此可以直接算块的 \(Q\) 并缓存。

对一个较小的整数 \(m\)，直接按二进制快速幂模拟即可：

  • 用集合 \(A\) 表示当前多项式的指数集合，初始 \(A={0}\)（多项式 \(1\)）
  • 扫描二进制位（从高到低）：
    1. 平方：指数翻倍\(A'\leftarrow\{2e:e\in A\}\)
    2. 如果当前位是 \(1\)，则再乘 \(P(x)=1+x+x^2+x^3+x^4\)，在集合上就是 XOR 叠加平移：\(A\leftarrow A\triangle (A+1)\triangle (A+2)\triangle (A+3)\triangle (A+4)\)
  • 最终 \(Q(m)=|A|\)

代码

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
from functools import lru_cache

# 计算一个“块”(二进制串)对应的 Q 值
# bits: 不含前导 0 的二进制字符串（通常也不含 "00"）
@lru_cache(None)
def q_block(bits: str) -> int:
    # A 代表当前多项式在 GF(2) 下的“非零项指数集合”
    A = {0}
    for c in bits:
        # 平方：指数翻倍
        A = {2 * e for e in A}

        # 若该位为 1：乘以 P(x)=1+x+x^2+x^3+x^4
        if c == '1':
            newA = set(A)
            for shift in (1, 2, 3, 4):
                for e in A:
                    t = e + shift
                    # XOR：出现奇数次保留，偶数次删除
                    if t in newA:
                        newA.remove(t)
                    else:
                        newA.add(t)
            A = newA

    return len(A)


def Q(n: int) -> int:
    """返回 Q(n)：(1+x+...+x^4)^n 的奇数系数个数。"""
    if n == 0:
        return 1  # 约定 P^0 = 1

    b = bin(n)[2:].rstrip('0')  # 忽略末尾 0（对应不断平方，不改变项数）
    if not b:
        return 1

    ans = 1
    for part in b.split('00'):
        part = part.strip('0')
        if part:
            ans *= q_block(part)
    return ans


total = 0
p = 10
for _ in range(18):
    total += Q(p)
    p *= 10
print(total)