Project Euler 701

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

Project Euler 701

题目

Random connected area

Consider a rectangle made up of \(W \times H\) square cells each with area \(1\).

Each cell is independently coloured black with probability 0.5 otherwise white. Black cells sharing an edge are assumed to be connected.

Consider the maximum area of connected cells.

Define \(E(W,H)\) to be the expected value of this maximum area. For example, \(E(2,2)=1.875\), as illustrated below.

You are also given \(E(4, 4) = 5.76487732\), rounded to \(8\) decimal places.

Find \(E(7, 7)\), rounded to \(8\) decimal places.

解决方案

由于每个格子独立且黑白概率都为 \(\dfrac12\)，所以所有 \(2^{WH}\) 种染色方案等概率。于是期望可写为\(\displaystyle E(W,H)=\frac{1}{2^{WH}}\sum_{\omega} M(\omega),\)其中 \(\omega\) 遍历全部染色方案，\(M(\omega)\) 表示该方案的最大黑色连通块面积。

因此问题的本质变成：对所有染色方案，精确统计其最大黑色连通块面积并求和。困难在于共有 \(2^{HW}\) 种方案，无法直接枚举。

可见，这题是一个典型的局部生成、全局连通的问题。如果我们按某种顺序逐格处理（例如从上到下、从左到右），那么每一步只新增一个格子的黑白信息。问题在于：虽然当前只处理一个格子，但它可能改变一个很大的连通块，甚至把两个大块合并成更大的块。

如果想在扫描过程中保持信息完整，就要回答一个关键问题：在当前，过去区域中哪些信息会影响未来，哪些信息已经不再重要。

对于连通块问题，可见未来只会通过当前处理边界与过去发生联系。因此，真正需要保留的是边界上的连通信息，而不是整个已处理区域的形状。这就是连通性剖面DP的核心出发点。

因此，我们固定按如下顺序扫描格子：按行扫描：先第 \(0\) 行从左到右，再第 \(1\) 行从左到右，……，直到第 \(H-1\) 行。以宽度 \(W\) 为例，在任意时刻，扫描线与尚未处理区域的交界可以用一个长度为 \(W\) 的数组表示。记为\(\mathrm{lab}[0],\mathrm{lab}[1],\dots,\mathrm{lab}[W-1].\)

它表示：在每一列 \(j\) 的前沿位置上，若这里对应某个仍可能与未来连接的黑色连通块，则记下该连通块的编号；若该列前沿处没有活动黑块，则记为 \(0\)。

于是\(\mathrm{lab}[j]=0\)表示该列没有活动连通块延伸到前沿；\(\mathrm{lab}[j]=k>0\)表示该列前沿属于编号为 \(k\) 的活动连通块。注意这里的活动连通块指的是：已经在已处理区域中出现，但仍然有机会在未来继续扩展的黑色连通块。此外，还需要记录每个活动编号 \(k\) 当前已经包含了多少黑格，记作 \(\mathrm{s}[k]\)。

注意，我们最终要的是全局最大黑色连通块面积。但扫描尚未结束时，一个连通块可能还在增长，所以它的最终大小未定。那么什么时候一个连通块的大小可以认为已经最终确定？

关键在于，当某个活动连通块的编号不再出现在前沿数组 \(\mathrm{lab}\) 中时，它就不可能再与未来任何格子相连。原因在于未来区域与已处理区域的唯一接触界面就是前沿；如果一个连通块在前沿上已经没有任何出口，它就彻底封闭了。

于是该连通块的最终大小已经确定，恰为当前记录的 \(\mathrm{s}[k]\)，可以拿来更新当前已知最大连通块面积。这就说明我们不需要保存所有已经完成的连通块的形状，只需要在它们退出前沿的那一刻，把大小并入最大值统计即可。这一步保证了扫描法不会漏掉任何连通块，也不会错误估计最大值。

现在我们把状态定义朴素的三部分（后文会优化）：前沿连通结构（即 \(\mathrm{lab}\)）；活动连通块大小（即 \(\mathrm{s}\)）；当前为止的最大连通块面积（记作 \(\mathrm{best}\)）。

那么这样确实可以做出一个完全正确的 DP：每次新加一个格子（黑或白）更新状态，最后把所有状态按方案数加权求和即可。但这个设计有一个致命性能问题：

相同的前沿结构与活动块大小（也就是未来转移方式完全相同），可能对应许多不同的 \(\mathrm{best}\)，如果把 \(\mathrm{best}\) 也放进状态键，就会把本来能合并的状态拆成很多份，状态数量暴涨。

也就是说，决定未来转移的信息与用于最终求值的信息被混在了一起，导致状态膨胀。因此我们需要进一步优化。

需要注意的是，未来如何转移，只取决于：前沿上各列属于哪些活动连通块（\(\mathrm{lab}\)），以及这些活动连通块当前多大（\(\mathrm{s}\)）。把这两部分统称为一个连通剖面。

现在，我们不再让每个剖面只对应一个方案数，而是对应一个数组\(\mathrm{cnt}_P[b]\)表示达到剖面\(P\)且当前最大连通块面积恰为\(b\)的方案数，\(0\le b\le WH.\)

这样做的效果是：哈希键只包含连通剖面（真正影响未来的部分），并且原本因 \(\mathrm{best}\) 不同而分裂的状态，被重新聚合到同一个键。最终，最大值信息作为附属分布放在值里，用长度仅 50 的数组维护即可。

另一个非常关键的优化是规范化编号。\(\text{lab}\) 数组里的连通块编号本质上只是“标签”，没有数学意义，因此像 \([2,2,0,5,5,0,0]\) 和 \([1,1,0,2,2,0,0]\) 这种状态虽然编号不同，但表示的前沿连通关系完全相同。为避免 DP 把这类同构状态重复计数，需要在每次转移后按从左到右的首次出现顺序，将非零编号依次重命名为 \(1,2,3,\dots\)，并同步重排对应的 \(\mathrm{s}\)。这样所有等价的连通结构都会被归并到同一个连通剖面，从而显著压缩状态数。

处理当前格子时，只需要考虑它与已处理区域的相邻关系。由于采用行优先扫描，当前格子只可能与两侧已处理格子相邻：左邻和上邻。设这两个方向对应的活动连通块编号分别为：左邻编号为\(L\)，上邻编号为\(U\)。其余方向（右、下）尚未处理，不会影响当前连通性。现在考虑如何进行转移。

若当前格子是白色，它不属于任何黑色连通块。这时前沿在该列的位置会被清空。若该列原本的上方编号为 \(U\)，则需要检查：清空之后，编号 \(U\) 是否还在前沿其他列中出现。

若不再出现，说明这个活动连通块已经彻底封闭，其最终大小为 \(\mathrm{s}[U]\)。此时：用 \(\mathrm{s}[U]\) 更新最大值（在分布层面体现为与 \(b\) 取 \(\max\)），然后该连通块从活动集合中移除。

这一步没有新增黑格，因此“当前格子形成/增大的连通块大小”可以看作 \(0\)。对最大值分布而言，只是可能因为封闭了某块而需要把它并入最大值统计（实现时可在剖面转移前或转移后处理，本质一致）。

当前格子为黑色时，它与左、上邻居中的黑色活动块发生连接。分类讨论如下。

  1. \(L=0,U=0.\)左、上都没有黑色活动块与之相连。于是当前格子形成一个新的黑色活动连通块，大小为 \(1\)。这一步会使当前可能的最大值至少达到 \(1\)。在分布上，相当于把每个 \(b\) 映射到 \(\max(b,1)\)。
  2. 恰有一个非零，例如 \(L\neq 0,U=0.\)当前格子接到某个已有活动连通块上，设编号为 \(x\)。那么该活动块大小增加 \(1\)，当前格子在前沿对应位置也标记为 \(x\)。如果更新后该活动块大小变成 \(m\)，则最大值分布的更新规则是：\(b \longmapsto \max(b,m).\)
  3. \(L=U\neq 0.\)左、上两个方向本来就属于同一个活动连通块。当前格子只是在该块内部继续扩展：活动块大小增加 \(1\)；并不会产生新的合并。设更新后大小为 \(m\)，最大值分布仍按\(b \longmapsto \max(b,m)\)更新。
  4. \(L\neq 0,U\neq 0,L\neq U.\)这是最关键的一类：当前格子把两个不同的活动连通块连接起来了。此时必须合并两个活动块。设它们大小分别为 \(\mathrm{s}[L]\) 和 \(\mathrm{s}[U]\)，那么合并后新块大小为\(m=\mathrm{s}[L]+\mathrm{s}[U]+1.\)并且前沿中所有属于这两个编号的位置，都必须统一成同一个编号（随后再规范化）。最大值分布仍按同样的规则更新：\(b \longmapsto \max(b,m).\)

设当前某个剖面 \(P\) 的分布为 \(\mathrm{cnt}_P[b]\)。若本次转移形成（或更新）的黑色连通块大小为 \(m\)，则转移后的分布 \(\mathrm{cnt}'[x]\) 满足：如果 \(m=0\)（例如当前格子为白，且不触发额外最大值提升），则分布不变；如果 \(m>0\)，则原来所有 \(b<m\) 的方案都会变成新最大值 \(m\)；原来所有 \(b\ge m\) 的方案最大值保持不变。

因此有：\(\displaystyle\mathrm{cnt}'[m]\leftarrow \sum_{b=0}^{m-1}\mathrm{cnt}_P[b],\)以及对所有 \(b\ge m\)，有\(\mathrm{cnt}'[b]\leftarrow \mathrm{cnt}_P[b].\)这实际上是在做一个按 \(\max\) 运算推动的分布更新。

代码

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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
#include <bits/stdc++.h>
using namespace std;

const int W = 7;
const int H = 7;
const int BMAX = W * H;
typedef long long ll;
struct Profile
{
    uint32_t front;
    uint64_t sizes;
};

struct ProfileHash
{
    size_t operator()(Profile const &p) const noexcept
    {
        uint64_t x = p.front;
        uint64_t y = p.sizes;
        uint64_t h = x * 0x9e3779b97f4a7c15ULL;
        h ^= y + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        return (size_t)h;
    }
};

struct ProfileEq
{
    bool operator()(Profile const &a, Profile const &b) const noexcept
    {
        return a.front == b.front && a.sizes == b.sizes;
    }
};

static inline int get_lab(uint32_t front, int i)
{
    return (front >> (3 * i)) & 7;
}

static inline uint32_t set_lab(uint32_t front, int i, int v)
{
    uint32_t mask = 7u << (3 * i);
    return (front & ~mask) | ((uint32_t)v << (3 * i));
}

static inline int get_size(uint64_t sizes, int lab)
{
    return (sizes >> (6 * (lab - 1))) & 63;
}

static inline uint64_t set_size(uint64_t sizes, int lab, int v)
{
    uint64_t mask = 63ull << (6 * (lab - 1));
    return (sizes & ~mask) | ((uint64_t)v << (6 * (lab - 1)));
}

static inline bool label_present(uint32_t front, int lab)
{
    for (int i = 0; i < W; i++)
        if (get_lab(front, i) == lab)
            return true;
    return false;
}

static inline int max_label(uint32_t front)
{
    int k = 0;
    for (int i = 0; i < W; i++)
        k = max(k, get_lab(front, i));
    return k;
}

static inline Profile canonical(uint32_t front, uint64_t sizes)
{
    int mp[8];
    int ns[8];
    memset(mp, 0, sizeof(mp));
    memset(ns, 0, sizeof(ns));
    int k = 0;
    int lab[W];
    for (int i = 0; i < W; i++)
    {
        int x = get_lab(front, i);
        if (x == 0)
        {
            lab[i] = 0;
            continue;
        }
        if (mp[x] == 0)
        {
            mp[x] = ++k;
            ns[k] = get_size(sizes, x);
        }
        lab[i] = mp[x];
    }
    uint32_t nfront = 0;
    for (int i = 0; i < W; i++)
        nfront |= (uint32_t)lab[i] << (3 * i);
    uint64_t nsizes = 0;
    for (int i = 1; i <= k; i++)
        nsizes |= (uint64_t)ns[i] << (6 * (i - 1));
    return Profile{nfront, nsizes};
}

static inline Profile trans_white(Profile p, int c)
{
    uint32_t front = p.front;
    uint64_t sizes = p.sizes;
    int U = get_lab(front, c);
    front = set_lab(front, c, 0);
    if (U != 0 && !label_present(front, U))
        sizes = set_size(sizes, U, 0);
    return canonical(front, sizes);
}

static inline pair<Profile, int> trans_black(Profile p, int c)
{
    uint32_t front = p.front;
    uint64_t sizes = p.sizes;
    int L = (c > 0) ? get_lab(front, c - 1) : 0;
    int U = get_lab(front, c);
    int k = max_label(front);
    int new_sz = 0;

    if (L == 0 && U == 0)
    {
        int nid = k + 1;
        sizes = set_size(sizes, nid, 1);
        front = set_lab(front, c, nid);
        new_sz = 1;
    }
    else if (L != 0 && U == 0)
    {
        int s = get_size(sizes, L) + 1;
        sizes = set_size(sizes, L, s);
        front = set_lab(front, c, L);
        new_sz = s;
    }
    else if (L == 0 && U != 0)
    {
        int s = get_size(sizes, U) + 1;
        sizes = set_size(sizes, U, s);
        front = set_lab(front, c, U);
        new_sz = s;
    }
    else
    {
        if (L == U)
        {
            int s = get_size(sizes, L) + 1;
            sizes = set_size(sizes, L, s);
            front = set_lab(front, c, L);
            new_sz = s;
        }
        else
        {
            int root = min(L, U);
            int other = max(L, U);
            int s = get_size(sizes, root) + get_size(sizes, other) + 1;
            sizes = set_size(sizes, root, s);
            sizes = set_size(sizes, other, 0);
            for (int i = 0; i < W; i++)
                if (get_lab(front, i) == other)
                    front = set_lab(front, i, root);
            front = set_lab(front, c, root);
            new_sz = s;
        }
    }

    Profile q = canonical(front, sizes);
    return {q, new_sz};
}

static inline void add_dist(array<ll, BMAX + 1> &dst, array<ll, BMAX + 1> const &src, int m)
{
    if (m <= 0)
    {
        for (int b = 0; b <= BMAX; b++)
            dst[b] += src[b];
        return;
    }
    ll pref = 0;
    for (int b = 0; b < m; b++)
        pref += src[b];
    dst[m] += pref;
    for (int b = m; b <= BMAX; b++)
        dst[b] += src[b];
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    unordered_map<Profile, array<ll, BMAX + 1>, ProfileHash, ProfileEq> dp, ndp;
    array<ll, BMAX + 1> init{};
    init[0] = 1;
    dp[canonical(0u, 0ull)] = init;

    for (int r = 0; r < H; r++)
    {
        for (int c = 0; c < W; c++)
        {
            ndp.clear();
            ndp.reserve(dp.size() * 2 + 16);
            for (auto const &it : dp)
            {
                Profile st = it.first;
                auto const &dist = it.second;

                Profile s0 = trans_white(st, c);
                auto &d0 = ndp[s0];
                add_dist(d0, dist, 0);

                auto t1 = trans_black(st, c);
                Profile s1 = t1.first;
                int m = t1.second;
                auto &d1 = ndp[s1];
                add_dist(d1, dist, m);
            }
            dp.swap(ndp);
        }
    }

    ll num = 0;
    for (auto const &it : dp)
    {
        auto const &dist = it.second;
        for (int b = 0; b <= BMAX; b++)
            num += (ll)dist[b] * (ll)b;
    }
    double ans = (double)num;
    for (int i = 0; i < W * H; i++)
        ans /= 2.0L;
    cout.setf(std::ios::fixed);
    cout << setprecision(8) << ans << "\n";
    return 0;
}