Project Euler 839

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

Project Euler 839

题目

Beans in Bowls

The sequence \(S_n\) is defined by \(S_0 = 290797\) and \(S_n = S_{n - 1}^2 \bmod 50515093\) for \(n > 0\).

There are \(N\) bowls indexed \(0,1,\dots ,N-1\). Initially there are \(S_n\) beans in bowl \(n\).

At each step, the smallest index \(n\) is found such that bowl \(n\) has strictly more beans than bowl \(n+1\). Then one bean is moved from bowl \(n\) to bowl \(n+1\).

Let \(B(N)\) be the number of steps needed to sort the bowls into non-descending order.
For example, \(B(5) = 0\), \(B(6) = 14263289\) and \(B(100)=3284417556\).

Find \(B(10^7)\).

解决方案

我们定义势函数\(\displaystyle{\Phi(a)=\sum_{i=0}^{N-1} i\cdot a_i.}\)一次操作只改变 \(a[n],a[n+1]\)：\(\Delta\Phi = n\cdot(-1) + (n+1)\cdot(+1)=1.\)所以每一步势函数恰好 \(+1\)。因此如果最终稳定状态为 \(b\)（非降序），那么\(B(N)=\Phi(b)-\Phi(a).\)

终止性也立刻得到：总豆数 \(T=\sum a_i\) 不变，而\(0\le \Phi(a)\le (N-1)\cdot T\)有上界，且每步 \(+1\)，因此必有限步停止。问题就转化为两件事：计算算初始 \(\Phi(a)\)。 求最终序列 \(b\) 并算 \(\Phi(b)\)。

我们引入前缀和\(\displaystyle{P_k=\sum_{i=0}^{k-1} a_i, k=0,1,\dots,N,}\)则\(a_i=P_{i+1}-P_i.\)非降条件 \(a_i\le a_{i+1}\) 等价于\(P_{i+1}-P_i \le P_{i+2}-P_{i+1}\iff 2P_{i+1}\le P_i+P_{i+2}.\)因此 稳定态 \(\Leftrightarrow\) 前缀和序列 \(P\) 满足离散凸性：\(2P_k\le P_{k-1}+P_{k+1} (1\le k\le N-1).\)

再看一次操作对 \(P\) 的影响：把 \(1\) 个豆从 \(n\) 移到 \(n+1\) 等价于前 \(n+1\) 个碗里的总豆数减少 \(1\)；其他前缀总豆数不变。也就是\(P_{n+1}\leftarrow P_{n+1}-1\)，其余\(P_k\)$。因此整个过程是在仅通过降低某些中间点 \(P_k\) 来把离散曲线压成凸的，且端点 \(P_0=0\)、\(P_N=T\) 始终固定。

这就解释了为什么最终状态具有很强的结构：它对应于在不超过初始前缀和曲线的前提下，取最大的凸序列（直观上就是把曲线向下压到刚好凸为止）。一旦得到最终凸的 \(P^{\ast}\)，答案序列就是\(b_i=P^{\ast}_{i+1}-P^{\ast}_i,\)并天然非降。

接下来要做的是在线性时间内直接构造最终的 \(b\)。

考虑把某一段连续位置（长度 \(l\)，总和 \(s\)）在最终状态里安排成非降且总和不变。若只要求尽量平，那么最平的整数形态必然是\(\lfloor s/l\rfloor\) 出现 \(l - (s \bmod l)\) 次，\(\lceil s/l\rceil=\lfloor s/l\rfloor+1\) 出现 \(s \bmod l\) 次，并且为了非降，所有较小值在前、较大值在后。

也就是说该段最终形态一定是\(q,q,\dots,q,q+1,\dots,q+1\)，其中 \(q=\left\lfloor \dfrac{s}{l}\right\rfloor\)。

从左到右处理原序列 \(a[i]\)。假设我们已经把前缀 \(a[0..i-1]\) 的最终形态用若干“等值块”表示出来（这些块按出现顺序排列，且块的取值非降）。现在加入新元素 \(x=a[i]\)：

  • 若 \(x\) 不小于当前分块表示中最后一个块的取值（也就是前缀最终形态的末尾值），显然可以直接在末尾新开一个取值为 \(x\) 的块（或把它并入末尾同值块），不会破坏非降。
  • 若 \(x\) 小于当前分块表示中最后一个块的取值，末尾出现下降，说明仅靠把 \(x\) 放在最后会违背“非降”。因此在最终状态里，必须把这一段尾部连同其左侧的一些块一起重新调整：在块的视角，就是把 \(x\) 与若干个相邻的末尾块合并成一个更长的连续区间，然后把该区间的总和按“尽量平均”的原则重新分配成若干个相等值（以及可能多出的 \(+1\)）的块，使得分配后的块值重新恢复为非降。

但合并多长才够？设当前准备重整的后缀段总和为 \(s\)、长度为 \(l\)，其平均值为 \(\mu=s/l\)。为了让它能接在左边块值 \(v\) 后面且不下降，必须保证该段的最小值 \(\lfloor\mu\rfloor\) 不小于 \(v\)。一个足够强、且便于判定的条件是\(v < \mu\)，这等价于\(v\cdot l < s.\)

如果 \(v\ge \mu\)，那么不管怎么在该段内部做“尽量平均”的整数分配，它都会产生某些值小于或等于 \(v\)，并且很容易违背单调非递减的性质；正确做法就是继续把左边块也并入后缀段，重新计算更大的 \(s,l\)，直到左边块值严格小于新平均值为止。

这就得到经典的单调栈的过程：

维护一个栈，元素为 \((v, c)\)（表示当前前缀最终态的若干等值块，值严格递增或递增合并）。处理新值 \(x\) 时，设 \(s=x,l=1\)，只要栈顶块值 \(v\) 满足\(v \ge \dfrac{s}{l}\iff v\cdot l \ge s,\)就把栈顶并入：\(s'=s+ v\cdot c,l'=l+c\)，并弹栈。结束后，把该合并段按尽量平均分裂为 \(q\) 与 \(q+1\) 两块压回栈。

栈最终给出最终序列 \(b\) 的“值-长度”分解。若某块值为 \(v\)，长度为 \(c\)，并且它覆盖的下标是 \(L,\dots,L+c-1\)，则\(\displaystyle{\sum_{i=L}^{L+c-1} i\cdot v= v\cdot \frac{(2L+c-1)c}{2}.}\)逐块累加即可算出 \(\Phi(b)\)。

初始 \(\Phi(a)=\sum i\cdot S_i\) 也能在线生成 \(S_i\) 时顺带累加，完全不需要存 \(a\)。最终\(B(N)=\Phi(b)-\Phi(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
52
53
54
55
56
57
58
59
60
61
N = 10 ** 7
MOD = 50515093

def solve(N: int) -> int:

    s = 290797

    vals = [s]
    cnts = [1]

    init_sum = 0

    for i in range(1, N):
        s = (s * s) % MOD
        init_sum += i * s

        last = vals[-1]
        if s > last:
            vals.append(s)
            cnts.append(1)
            continue
        if s == last:
            cnts[-1] += 1
            continue

        total = s
        cnt = 1
        while vals and vals[-1] * cnt >= total:
            total += vals[-1] * cnts[-1]
            cnt += cnts[-1]
            vals.pop()
            cnts.pop()

        q, r = divmod(total, cnt)
        n0 = cnt - r

        if n0:
            if vals and vals[-1] == q:
                cnts[-1] += n0
            else:
                vals.append(q)
                cnts.append(n0)

        if r:
            qp1 = q + 1
            if vals and vals[-1] == qp1:
                cnts[-1] += r
            else:
                vals.append(qp1)
                cnts.append(r)

    final_sum = 0
    L = 0
    for v, c in zip(vals, cnts):
        final_sum += v * (2 * L + c - 1) * c // 2
        L += c

    return final_sum - init_sum


print(solve(N))