Project Euler 767

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

Project Euler 767

题目

Window into a Matrix II

A window into a matrix is a contiguous sub matrix.

Consider a \(16\times n\) matrix where every entry is either \(0\) or \(1\). Let \(B(k,n)\) be the total number of these matrices such that the sum of the entries in every \(2\times k\) window is \(k\).

You are given that \(B(2,4) = 65550\) and \(B(3,9) \equiv 87273560 \pmod{1\,000\,000\,007}\).

Find \(B(10^5,10^{16})\). Give your answer modulo \(1\,000\,000\,007\).

解决方案

令\(R=16\)。固定某个相邻行对 \((i,i+1)\)（其中 \(0\le i\le 14\)），以及窗口起点列 \(c\)（其中 \(0\le c\le n-k\)），窗口和为\(\displaystyle{\sum_{t=0}^{k-1}\bigl(A_{i,c+t}+A_{i+1,c+t}\bigr)=k.}\)

将窗口右移一列（从 \(c\) 变为 \(c+1\)），两边相减，消去公共部分得到：\(A_{i,c}+A_{i+1,c}=A_{i,c+k}+A_{i+1,c+k}.\)因此对每个固定的行对 \((i,i+1)\)，序列\(s_{c}=A_{i,c}+A_{i+1,c}\in\{0,1,2\}\)满足\(s_{c}=s_{c+k},\)即在列方向具有周期 \(k\)。

令两列分别为长度 \(R\) 的 \(0/1\) 向量 \(v=(v_0,\dots,v_{R - 1})\) 与 \(w=(w_0,\dots,w_{R - 1})\)，来自列 \(c\) 与列 \(c+k\)。上面的结论对所有 \(i=0,\dots,R - 2\) 给出\(v_i+v_{i+1}=w_i+w_{i+1}.\)

令差分 \(d_i=v_i-w_i\in\{-1,0,1\}\)，则\((v_i+v_{i+1})-(w_i+w_{i+1})=d_i+d_{i+1}=0,\)所以\(d_{i+1}=-d_i (0\le i\le R - 2).\)

因此 \(\{d_i\}\) 必须交替取相反数。只有两种可能：

  1. \(d_i\equiv 0\)，即 \(v=w\)（两列完全相同）；

  2. \(d_i\) 在 \(\pm1\) 间交替，这将强迫 \(v,w\) 都是严格交替的列：\(0101\cdots\)或\(1010\cdots,\)并且 \(w\) 恰好是 \(v\) 的逐位翻转。

于是得到关键分类：

  • 若某列 不是严格交替列，那么同余类中所有列都必须与它 完全相同；
  • 若某列是严格交替列，则同余类中每一列都只能是两种交替列之一（\(0101\cdots\) 或 \(1010\cdots\)）。

假设 \(k\mid n\)，令\(m=\dfrac{n}{k}.\)每个同余类（余数 \(r\in{0,1,\dots,k-1}\)）对应 \(m\) 列：\(r,r+k,r+2k,\dots,r+(m-1)k.\)

此前可知：每个同余类只有两种形态：

  • 等列类：该同余类的 \(m\) 列全都相同（列向量任意）；
  • 交替类：该同余类的每一列都必须是两种交替列之一（\(0101\cdots\) 或 \(1010\cdots\)）。

下面按“等列类”的数量计数。设有 \(l\) 个同余类是等列类，则共有\(\dbinom{k}{l}\)种选择方式。

剩下 \(k-l\) 个同余类为交替类。对一个交替类来说，\(m\) 列各自可独立选两种交替列之一，共 \(2^m\) 种，但其中有 \(2\) 种是“全都选同一种交替列”，这会使该同余类实际上也成为等列类（会被计入 \(l\) 的情况），因此必须剔除这两种：\(2^m-2\)。所以这部分贡献为\((2^m-2)^{k-l}.\) 。

现在只看那 \(l\) 个等列类形成的“前 \(k\) 列代表块”中的对应 \(l\) 列，组成一个 \(R\times l\) 的 \(0/1\) 矩阵 \(C\)。

对任意相邻行对 \((i,i+1)\)，一个 \(2\times k\) 窗口的总和必须是 \(k\)。其中交替类的每一列在任意相邻行对上贡献恒为 \(1\)，因此 \(k-l\) 个交替类总贡献为 \(k-l\)，那么等列类必须补足剩余的 \(l\)：\(\displaystyle{\sum_{j=0}^{l-1}\bigl(C_{i,j}+C_{i+1,j}\bigr)=l, i=0,1,\dots,R - 2.}\)

令第 \(i\) 行在等列类中的行和为\(\displaystyle{r_i=\sum_{j=0}^{l-1}C_{i,j},}\)则约束变为\(r_i+r_{i+1}=l.\)这意味着行和交替互补：\(r_{i+1}=l-r_i.\)

对所有奇数下标行（\(i=1,3,5,\dots,R - 1\)）做一次按位翻转（\(0\leftrightarrow 1\)），翻转后的这些行行和变为 \(l-r_i\)，于是每一行的行和都统一为同一个值 \(s\)。翻转是双射，所以计数不变。

因此等列类的合法选择数等于：对某个 \(s\in[0,l]\)，每一行独立从 \(l\) 列中选出 \(s\) 个位置放 \(1\)，共有 \(\dbinom{l}{s}^{R}\) 种；再对所有 \(s\) 求和：\(\displaystyle{\sum_{s=0}^{l}\binom{l}{s}^{R}.}\)

综合以上三部分，得到

\[ B(k,n)=\sum_{l=0}^{k}\binom{k}{l}(2^{n/k}-2)^{k-l}\sum_{s=0}^{l}\binom{l}{s}^{R}. \]

注意，这个公式的使用要满足\(k\mid n\)。

接下来讲解用卷积在\(O(k\log k)\)的时间计算\(\displaystyle{\sum_{s=0}^{l}\binom{l}{s}^{R}}\)。

记\(\displaystyle{F(l)=\sum_{s=0}^{l}\binom{l}{s}^{R}.}\)将组合数展开：\(\dbinom{l}{s}^{R}=\left(\dfrac{l!}{s!(l-s)!}\right)^{R}=\dfrac{l!^{R}}{s!^{R}(l-s)!^{R}}.\)于是\(\displaystyle{F(l)=l!^{R}\sum_{s=0}^{l}\frac{1}{s!^{R}(l-s)!^{R}}.}\)

令\(a_s=\dfrac{1}{s!^{R}},\)则内层和正是卷积系数：\(\displaystyle{c_l=\sum_{s=0}^{l}a_sa_{l-s}.}\)因此\(F(l)=l!^{R}c_l.\)只需构造多项式\(\displaystyle{A(x)=\sum_{s=0}^{k}a_s x^s,}\)则\(\displaystyle{A(x)^2=\sum_{l=0}^{2k}c_l x^l,}\) 所有 \(c_0,\dots,c_k\) 可以通过一次多项式平方在 \(O(k\log k)\) 得到。

最终需要在模\(P=1000000007\)下计算，但 \(P\) 不是 NTT 友好素数。做法是选取三个 NTT 友好素数（都形如 \(p=t\cdot 2^q+1\)）：

\[ p_1=998244353, p_2=985661441, p_3=943718401, \]

分别在这三个模数下做 NTT 得到卷积系数模 \(p_i\) 的结果，再用中国剩余定理把每个系数还原为一个小于 \(p_1p_2p_3\) 的整数。由于这里卷积系数的真实大小远小于 \(p_1p_2p_3\)，因此可视为精确还原，最后再对 \(P\) 取模即可。

代码

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
from flint import nmod_poly
from tools import CRT

K = 10 ** 5
N = 10 ** 16
MOD = 1000000007

def conv3_flint(a, k):
    p1 = 998244353
    p2 = 985661441
    p3 = 943718401

    A1 = nmod_poly(a, p1)
    A2 = nmod_poly(a, p2)
    A3 = nmod_poly(a, p3)

    C1 = A1 * A1
    C2 = A2 * A2
    C3 = A3 * A3

    c1 = [int(x) for x in list(C1)[: k + 1]]
    c2 = [int(x) for x in list(C2)[: k + 1]]
    c3 = [int(x) for x in list(C3)[: k + 1]]
    mods = [p1, p2, p3]
    res = [0] * (k + 1)
    for i in range(k + 1):
        res[i] = CRT(mods, [c1[i], c2[i], c3[i]]) % MOD
    return res

def B(k: int, n: int) -> int:
    fac = [1] * (k + 1)
    for i in range(1, k + 1):
        fac[i] = fac[i - 1] * i % MOD
    ifac = [1] * (k + 1)
    ifac[k] = pow(fac[k], MOD - 2, MOD)
    for i in range(k, 0, -1):
        ifac[i - 1] = ifac[i] * i % MOD

    fac16 = [0] * (k + 1)
    invfac16 = [0] * (k + 1)
    for i in range(k + 1):
        fac16[i] = pow(fac[i], 16,  MOD)
        invfac16[i] = pow(ifac[i],16, MOD)

    conv = conv3_flint(invfac16, k)

    F = [0] * (k + 1)
    for l in range(k + 1):
        F[l] = fac16[l] * conv[l] % MOD

    m = n // k
    x = (pow(2, m, MOD) - 2) % MOD

    powx = [1] * (k + 1)
    for i in range(1, k + 1):
        powx[i] = powx[i - 1] * x % MOD

    fk = fac[k]
    ans = 0
    for l in range(k + 1):
        choose = fk * ifac[l] % MOD * ifac[k - l] % MOD
        ans = (ans + choose * powx[k - l] % MOD * F[l]) % MOD
    return ans


print(B(K, N))