Blog. You could have invented Fenwick trees [blog/fenwick]
- CAIMEOX
- 2025-05-08
Blog. You could have invented Fenwick trees [blog/fenwick]
- CAIMEOX
- 2025-05-08
本文译自 Brent Yorgey 的 Functional Pearl 论文 You could have invented Fenwick trees,并做了一些解释补充、错误修正,以及将原文的 Haskell 代码都翻译到了 MoonBit 代码。
芬威克树(Fenwick trees),亦称二叉索引树(binary indexed trees),是一种精巧的数据结构,旨在解决这样一个问题:如何在维护一个数值序列的同时,支持在亚线性时间内完成更新操作与区间查询。其实现简洁高效,然亦颇为费解,多由一些针对索引的、不甚直观的位运算构成。本文将从线段树(segment trees)入手——这是一种更为直接、易于验证的纯函数式解法——并运用等式推理,阐释芬威克树的实现乃是一种优化变体,此过程将借助一个 MoonBit 嵌入式领域特定语言(EDSL)来处理无限二进制补码数。
假设我们拥有一个包含 $n$ 个整数的序列 $a_1, a_2, \ldots, a_n$,并希望能够任意交错执行以下两种操作,如图 1 所示: 注意,更新操作表述为将现有值增加 $v$;我们亦可通过增加 $v - u$ (其中 $u$ 为旧值) 来将给定索引处的值设为一个新值 $v$。 倘若我们仅将整数存储于一个可变数组中,那么更新操作可在常数时间内完成,但区间查询则需要与区间大小成线性关系的时间,因其必须遍历整个区间 $[i, j]$ 以累加其值。 为改善区间查询的运行时效,我们或可尝试缓存(至少部分)区间和。然则,此举必须审慎为之,因为在更新某个索引处的值时,相关的缓存和亦须保持同步更新。例如,一种直截了当的方法是使用一个数组 $P$,其中 $P_i$ 存储前缀和 $a_1 + \cdots + a_i$;$P$ 可通过一次扫描在常数时间内预计算得出。如此一来,区间查询便十分快捷:我们可通过计算 $P_j - P_{i-1}$ 在常数时间内获得 $a_i + \cdots + a_j$(为方便起见,我们设 $P_0 = 0$,如此该式即便在 $i=1$ 时亦能成立)。不幸的是,此时更新操作却需要线性时间,因为改变 $a_i$ 将导致对所有 $j \ge i$ 的 $P_j$ 进行更新。 是否存在一种数据结构,能够使这两种操作均在亚线性时间内完成?(在继续阅读下一段之前,您不妨稍作停顿,思考一番!)这并非仅仅是一个学术问题:该问题最初是在算术编码(Rissanen & Langdon, 1979; Bird & Gibbons, 2002)的背景下被提出的。算术编码是一类将消息转换为比特序列以供存储或传输的技术。为最大限度地减少所需的比特数,通常需要为出现频率较高的字符分配较短的比特序列,反之亦然;这就引出了维护一个动态字符频率表的需求。每当处理一个新字符时,我们就更新该表;而为了将单位区间细分为与各字符频率成比例的连续段,我们需要查询该表以获取累积频率(Ryabko, 1989; Fenwick, 1994)。 那么,我们能否使两种操作都在亚线性时间内完成呢?答案自然是肯定的。一种简单的技巧是将序列划分为 $\sqrt{n}$ 个桶,每个桶的大小为 $\sqrt{n}$,并创建一个大小为 $\sqrt{n}$ 的附加数组来缓存每个桶的和。更新操作仍然保持在 $O(1)$ 时间内,因为我们只需更新给定索引处的值以及相应桶的和。区间查询现在则可在 $O(\sqrt{n})$ 时间内完成:为求和 $a_i + \cdots + a_j$,我们手动累加从 $a_i$ 到其所在桶末尾的值,以及从 $a_j$ 所在桶开头到 $a_j$ 的值;对于所有介于两者之间的桶,我们则可以直接查找它们的和。 我们可以通过引入额外的缓存层级,以略微增加更新操作的开销为代价,进一步加快区间查询的速度。例如,我们可以将序列划分为 $\sqrt[3]{n}$ 个“大桶”,然后将每个大桶进一步细分为 $\sqrt[3]{n}$ 个“小桶”,每个小桶包含 $\sqrt[3]{n}$ 个值。每个桶的和都被缓存;现在每次更新需要修改三个值,而区间查询则在 $O(\sqrt[3]{n})$ 时间内完成。 以此类推,我们最终会得到一种基于二分治策略的区间和缓存方法,使得更新和区间查询操作的时间复杂度均为2 $O(\log n)$ 。具体而言,我们可以构建一棵完全二叉树3,其叶节点存储序列本身,而每个内部节点则存储其子节点的和。(这对于许多函数式程序员而言应是一个耳熟能详的概念;例如,手指树(finger trees)(Hinze & Paterson, 2006; Apfelmus, 2009)便采用了类似的缓存机制。)由此产生的数据结构通常被称为线段树,4 大抵是因为每个内部节点最终都缓存了底层序列的一个(连续)段的和。 图 2 展示了一棵基于长度为 $n=16$ 的示例数组构建的线段树(为简化起见,我们假设 $n$ 是 2 的幂,尽管将其推广到非 2 的幂的情况亦不难)。树的每个叶节点对应数组中的一个元素;每个内部节点都带有一个灰色条块,显示其所代表的底层数组段。 我们来看如何利用线段树来实现前面所述的两种操作,并使其均能在对数时间内完成。 图 4 阐释了计算区间 $[4 \ldots 11]$ 之和的过程。蓝色节点是我们递归遍历的节点;绿色节点是其范围完全包含在查询范围内、无需进一步递归而直接返回其值的节点;灰色节点则与查询范围不相交,返回零。本例中的最终结果是绿色节点值的总和,$1 + 1 + 5 + (-2) = 5$(不难验证,这确实是区间 $[4 \ldots 11]$ 内的值的和)。 在这个小例子中,看似我们访问了相当大比例的节点,但一般而言,我们访问的节点数量不会超过大约 $4 \log n$ 个。图 5 更清晰地展示了这一点。整棵树中最多只有一个蓝色节点可以拥有两个蓝色子节点,因此,树的每一层最多包含两个蓝色节点和两个非蓝色节点。我们实际上执行了两次二分查找,一次用于寻找查询区间的左端点,另一次用于寻找右端点。 线段树是解决该问题的一种颇为优良的方案:正如我们将在第 2 节中看到的那样,它们与函数式语言契合良好;同时它们也易于进行强大的泛化,例如支持惰性传播的区间更新以及通过共享不可变结构实现持久化更新历史(Ivanov, 2011b)。 芬威克树,或称二叉索引树(Fenwick, 1994; Ivanov, 2011a),是该问题的另一种解决方案。它们在通用性上有所欠缺,但在内存占用方面却极其精简——除了存储树中值的数组之外,几乎不需要任何额外空间——并且实现速度极快。换言之,它们非常适用于诸如底层编码/解码例程之类的应用场景,这些场景不需要线段树所能提供的任何高级特性,但却追求极致的性能。 本节展示了一个典型的 Java 语言芬威克树实现。可见,其实现异常简洁,大体由若干小型循环构成,每个循环迭代仅执行少量的算术与位运算。然而,这段代码究竟在做什么、其工作原理为何,却不甚明了!稍加审视, 我们的目标,并非为此已解决的问题编写优雅的函数式代码。相反,我们的目标在于运用一种函数式的领域特定语言来处理位串,并结合等式推理,从第一性原理出发,推导并阐释这段令人费解的命令式代码——展示函数式思维与等式推理之力,即便用于理解以其他非函数式语言编写之代码亦能奏效。在对线段树(第 2 节)建立更深入的直觉之后,我们将看到芬威克树如何能被视为线段树的一种变体(第 3 节)。随后,我们将绕道探讨二进制补码表示法,开发一个适宜的位操作 DSL,并解释 值得注意的是,本文及后续部分均采用基于 1 的索引方式,即序列中的首个元素索引为 1。此选择之缘由,后文将予以阐明。 本文的 $\log$ 表示以 2 为底的对数,原文使用 $\lg$。 原文是 balanced binary tree,翻译为平衡二叉树或有歧义,此处使用完全二叉树以避免歧义。 此处术语的使用存在一些混淆。截至本文撰写之时,维基百科关于线段树(Wikipedia Contributors, 2024)的条目讨论的是一种用于计算几何的区间数据结构。然而,谷歌搜索“segment tree”的大部分结果都来自算法竞赛领域,其中它指的是本文所讨论的数据结构(例如,参见 Halim et al., 2020, Section 2.8 或 Ivanov, 2011b)。这两种数据结构基本上是不相关的。 1. Introduction [blog/fenwick/introduction]
Code 1.1. Implementing Fenwick trees with bit tricks in Java [blog/fenwick/java_fenwick]
class FenwickTree {
private long[] a;
public FenwickTree(int n) { a = new long[n+1]; }
public long prefix(int i) {
long s = 0;
for (; i > 0; i -= LSB(i)) s += a[i]; return s;
}
public void update(int i, long delta) {
for (; i < a.length; i += LSB(i)) a[i] += delta;
}
public long range(int i, int j) {
return prefix(j) - prefix(i-1);
}
public long get(int i) { return range(i,i); }
public void set(int i, long v) { update(i, v - get(i)); }
private int LSB(int i) { return i & (-i); }
}
range
、get
和 set
函数尚属直观,但其余函数则令人困惑。我们可以观察到,prefix
和 update
函数均调用了另一个名为 LSB
的函数,该函数不知何故对一个整数及其负数执行了按位逻辑与运算。事实上,LSB(x)
计算的是 $x$ 的最低有效位(least significant bit),即它返回最小的 $2^k$,使得 $x$ 的第 $k$ 位为 1。然而,LSB
的实现原理,以及最低有效位在此处用于计算更新和前缀和的方式与缘由,均非显而易见。LSB
函数的实现(第 4 节)。借助于该 DSL,我们将推导出在芬威克树与标准二叉树之间来回转换的函数(第 5 节)。最终,我们将能够推导出在芬威克树内部移动的函数:先将其转换为二叉树索引,执行显而易见的操作以实现二叉树内的预期移动,然后再转换回来。通过等式推理将转换过程融合消除,最终将揭示隐藏其后的LSB
函数,一如预期(第 6 节)。
下面代码展示了一个在 MoonBit 中实现的简单线段树,其中使用了一些处理索引范围的工具函数,如图 8 所示。我们将线段树存储为一个递归的代数数据类型。
并使用与前一节中递归描述直接对应的代码实现了 区间工具函数: 尽管此实现简单且相对易于理解,但与仅将值序列存储于数组中相比,它引入了相当大的开销。我们可以通过将线段树的所有节点存储在一个数组中,采用如图 9 所示的标准从左至右、广度优先的索引方案,来更巧妙地利用空间(例如,这种方案或其类似方案常用于实现二叉堆)。根节点标号为 1;每当我们向下移动一层,我们就在现有二进制表示后追加一位:向左子节点移动时追加 0,向右子节点移动时追加 1。 因此,每个节点以二进制表示的索引记录了从根节点到达该节点的路径上左右选择的序列。从一个节点移动到其子节点,只需执行一次左位移并视情况加 1 即可;从一个节点移动到其父节点,则执行一次右位移。这便定义了一个从正自然数到无限二叉树节点的双射。若我们将线段树数组标记为 $s_1 \ldots s_{2n-1}$,那么 $s_1$ 存储所有 $a_i$ 的和,$s_2$ 存储 $a_i$ 前半部分的和,$s_3$ 存储后半部分的和,依此类推。$a_1 \ldots a_n$ 本身则存储为 $s_n \ldots s_{2n-1}$。 关键在于,既然沿树递归下降对应于对索引的简单操作,我们所讨论的所有算法都可以直接转换为处理(可变)数组的代码:例如,我们不再存储指向当前子树的引用,而是存储一个整数索引;每当需要向左或向右下降时,我们只需将当前索引乘以 2,或者乘以 2 再加 1。以数组存储树节点还带来了额外的可能性:我们不必总是从根节点开始向下递归,而是可以从某个感兴趣的特定索引出发,反向沿树向上移动。 那么,我们如何从线段树过渡到芬威克树呢?我们从一个看似无关紧要的观察开始:并非所有存储在线段树中的值都是必需的。当然,从某种意义上说,所有非叶节点都是“不必要的”,因为它们代表的是可以轻易从原始序列重新计算出来的缓存区间和。这正是线段树的核心思想:缓存这些“冗余”的和,以空间换时间,使我们能够快速执行任意的更新和区间查询,代价是将所需的存储空间加倍。 但这并非我所指!实际上,存在另一组我们可以舍弃的值,并且这样做仍然能够保持更新和区间查询的对数运行时效。舍弃哪些值呢?很简单:只需舍弃每个右子节点中的数据即可。图 10 展示了我们一直在使用的示例树,但其中每个右子节点的数据已被删除。注意,“每个右子节点”既包括叶节点也包括内部节点:我们舍弃与其父节点关系为右子节点的所有节点相关联的数据。我们将数据被丢弃的节点称为非活跃(inactive)节点,其余节点(即左子节点和根节点)称为活跃(active)节点。我们亦称,以这种方式使其所有右子节点都变为非活跃状态的树为疏树(thinned trees)。 更新一棵疏树颇为容易:只需像以前一样更新相同的节点,忽略对非活跃节点的任何更新即可。但我们如何应答区间查询呢?不难看出,剩余的信息足以重构被丢弃的信息(您或许愿意尝试说服自己这一点:能否在不参考任何先前图示的情况下,推断出图 10 中灰色节点应有的值?)。然而,仅此观察本身并不能为我们提供一个良好的计算区间和的算法。 关键在于考虑前缀和。正如我们在引言以及 2.1 中 给定一棵疏树,原始数组的任意前缀和(并由此引申出任意区间和)均可在对数时间内计算得出,且仅需使用活跃节点的值。 证明 出人意料的是,在处理前缀查询这种特殊情况时,第 1 节中描述并在 2.1 中实现的原始区间查询算法竟能原封不动地奏效!也就是说,当遇到当前节点的范围完全包含在查询范围内的基准情形时——此时我们会返回当前节点的值——这种情况只会在活跃节点上发生。 首先,根节点本身是活跃的,因此查询整个范围是有效的。接下来,考虑我们在某个节点上并递归到其两个子节点的情形。左子节点总是活跃的,因此我们只需考虑递归到右子节点的情况。右子节点的范围不可能完全包含在查询范围内:由于查询范围总是形如 $[1, j]$ 的前缀,若右子节点的范围完全包含在 $[1, j]$ 内,那么左子节点的范围也必定如此——这意味着父节点的范围(即其子节点范围的并集)也必定完全包含在查询范围内。但在那种情况下,我们会直接返回父节点的值,而不会递归到右子节点。因此,当我们确实递归到一个右子节点时,我们最终可能会返回 0,或者可能进一步递归到它的两个孙子节点,但无论如何,我们绝不会尝试去查看右子节点本身的值。 图 11 阐释了在线段树上执行前缀查询的过程。注意,被访问的右子节点要么是蓝色要么是灰色;唯一的绿色节点都是左子节点。 2. Segment tree [blog/fenwick/segtree]
Code 2.1. Simple segment tree implementation in MoonBit [blog/fenwick/segtree_code]
update
和 rq
;随后,get
和 set
亦可基于它们来实现。将此代码推广至适用于存储来自任意交换幺半群(若不需要 set
操作)或任意阿贝尔群(即带逆元的交换幺半群,若需要 set
操作)的值的线段树亦非难事——但为简单起见,我们在此保持原状,因为这种推广对我们的主线故事并无增益。enum SegTree {
Empty
Branch(Int, Range, SegTree, SegTree)
} derive(Eq, Show)
pub fn SegTree::update(self : SegTree, i : Index, v : Int) -> SegTree {
match self {
Empty => Empty
Branch(a, rng, l, r) if rng.contains(i) =>
Branch(a + v, rng, l.update(i, v), r.update(i, v))
_ => self
}
}
fn rq(self : SegTree, q : Range) -> Int {
match self {
Empty => 0
Branch(a, rng, l, r) =>
if disjoint(rng, q) {
0
} else if subset(rng, q) {
a
} else {
l.rq(q) + r.rq(q)
}
}
}
pub fn SegTree::get(self : SegTree, i : Index) -> Int {
self.rq((i, i))
}
pub fn SegTree::set(self : SegTree, i : Index, v : Int) -> SegTree {
self.update(i, v - self.get(i))
}
typealias Index = Int
priv type Range (Index, Index) derive(Eq, Show)
fn subset(r1 : Range, r2 : Range) -> Bool {
let (a1, b1) = r1._
let (a2, b2) = r2._
a1 >= a2 && b1 <= b2
}
fn contains(self : Range, i : Index) -> Bool {
subset((i, i), self)
}
fn disjoint(r1 : Range, r2 : Range) -> Bool {
let (a1, b1) = r1._
let (a2, b2) = r2._
b1 < a2 || b2 < a1
}
range
函数的实现所见,如果我们能够计算任意 $k$ 的前缀和 $P_k = a_1 + \cdots + a_k$,那么我们就能通过 $P_j - P_{i-1}$ 来计算区间和 $a_i + \cdots + a_j$。Theorem 2.2. [blog/fenwick/thm/thm1]
我们应如何在内存中实际存储一棵删减后的线段树?如果我们仔细观察图 10,一种策略便显而易见:只需将每个活跃节点“向下滑动”并向右移动,直至其落入底层数组的一个空位中,如图 12 所示。这在活跃节点与范围 $1 \ldots n$ 内的索引之间建立了一一对应的关系。理解这种索引方案的另一种方式是使用树的后序遍历,跳过非活跃节点,并为遍历过程中遇到的活跃节点赋予连续的索引。我们也可以通过以“右倾”风格绘制树来可视化这一结果(图 13),将每个活跃节点与其存储位置所在的数组槽垂直对齐。 这种将删减后线段树中的活跃节点存储在数组中的方法,正是所谓的芬威克树。有时我也会称其为芬威克数组,意在特别强调其底层的数组数据结构。 诚然,这是一种巧妙的空间利用方式,但关键问题在于如何实现更新和区间查询操作。我们为线段树实现的这些操作,无论是直接在递归数据结构上操作,还是在树存储于数组中时通过对索引进行简单操作,都是通过递归地沿树向下遍历来完成的。然而,当将删减后树的活跃节点存储在芬威克数组中时,哪些数组索引上的操作能够对应于在树中移动,这并非先验可知。为攻克此难题,我们首先需要绕道进入一个用于处理二进制补码数值的领域特定语言。 3. Fenwick trees [blog/fenwick/fenwick_tree]
通常用于实现芬威克树的位技巧依赖于二进制数的二进制补码表示法,该表示法允许以统一的方式表示正数和负数;例如,一个全由 1 组成的位串表示 -1。因此,我们现在转向开发一个嵌入 MoonBit 的领域特定语言,用于操作二进制补码表示。 首先,我们定义一个位 (Bit) 的类型,附带用于求反、逻辑与和逻辑或的函数(注:还有 接下来,我们必须定义位串,即位的序列。与其固定一个特定的位宽,不如使用无限位串会更加优雅1。使用 Lazy List 来表示潜在的无限位串似乎颇具诱惑力,但这会导致一系列问题。例如,无限列表的相等性是不可判定的,而且通常没有办法将一个无限的位列表转换回一个 实际上,我们想要的位串是那些最终恒定的串,即那些最终趋于无限长的全零尾部(代表非负整数)或全一尾部(代表负整数)的串。每个这样的串都有一个有限的表示,因此在 MoonBit 中直接编码最终恒定的位串,不仅能去除垃圾,还能辅助编写一些优雅且停机的算法。 尽管我们已经消除了垃圾值,但仍存在一个问题:同一个值可能存在多个不同的表示。例如, 使用 让我们从一些用于在 trait 让我们用 QuickCheck 来测试一下,验证我们的转换函数: 现在我们可以开始在 一个位序列的最低有效位(Least Significant Bit, LSB)可以定义如下: 注意,我们为 按位逻辑与可以直截了当地定义。注意,我们只需要两种情况;如果输入的有限部分长度不同,与 按位取反同样显然: 上述函数遵循熟悉的模式。我们可以很容易地将它们推广到作用于任意元素类型的最终恒定流,然后用泛型的 我们用通常的进位传播算法来实现加法,并附带一些针对 不难发现这个加法定义是能够停机并产生正确结果的;但我们也可以通过用 QuickCheck 来尝试一下,从而获得相当的信心: 最后,下面这个求反的定义对于任何学过二进制补码算术的人来说可能都很熟悉;我将其留作一个练习,供感兴趣的读者证明对于所有 现在我们拥有了揭开芬威克树实现中第一个谜团的工具。 $$ \forall x : \text{Bits}. \quad \text{lsb}(x) = \text{and}(x, \text{neg}(x)) $$ 证明 通过对 $x$ 进行归纳(译者注:为了便利证明书写,我们将 最后,为了表达我们将在下一节中开发的索引转换函数,我们的 DSL 中还需要一些其他东西。首先是一些用于设置和清除单个位以及测试特定位是否已设置的函数: 我们还需要的是左移和右移,以及一个通用的 部分读者或许能认出无限二进制补码位串即为二进数,也即 p 进数在特定情况 $p=2$ 时的形式,但我们的论述并不依赖于理解此关联。 4. Two’s complement binary [blog/fenwick/two_complement]
to_enum
和 from_enum
的实现,用于实现和 Int
类型的互转,在此省略):enum Bit {
O
I
} derive(Eq)
pub fn Bit::not(self : Bit) -> Bit {
match self {
O => I
I => O
}
}
impl BitAnd for Bit with land(self, other) {
match (self, other) {
(O, _) => O
(I, x) => x
}
}
impl BitOr for Bit with lor(self, other) {
match (self, other) {
(I, _) => I
(O, x) => x
}
}
Int
—— 我们如何知道何时停止?(译者注:Lazy List 在 MoonBit 中也容易出现栈溢出的问题,比原文使用的 Haskell List
还要更坏一些) 事实上,这些实际问题源于一个更根本的问题:无限位列表对于二进制补码位串来说是一个糟糕的表示,因为它包含“垃圾”,即那些不对应于我们预期语义域中值的无限位列表。例如,cycle([I, O])
是一个永远在 I
和 O
之间交替的无限列表,但它并不代表一个有效的二进制补码整数编码。更糟糕的是非周期列表,比如在每个素数索引处为 I
而其他地方均为 O
的列表。enum Bits {
Rep(Bit)
Snoc(Bits, Bit)
} derive(Eq)
Rep(b)
代表一个由无限个位 b
组成的序列,而 Snoc(bs, b)
则代表位串 bs
后跟一个末位 b
。我们使用 Snoc
而非 Cons
,是为了与我们通常书写位串的方式相匹配,即最低有效位在最后。另请注意在 Snoc
的 Bits
字段上使用了严格性注解;这是为了排除仅使用 Snoc
构成的无限位列表,例如 bs = Snoc(Snoc(bs, O), I)
。换言之,构造类型为 Bits
的非底层值的唯一方法是拥有一个有限的 Snoc
序列,最终由一个 Rep
终止。Snoc(Rep(O)), O)
和 Rep(O)
都代表包含全零的无限位串。在 Haskell 中可以通过精心构造一个双向模式(Pickering et al., 2016)来解决这个问题,但是 MoonBit 并没有 ViewPatterns
的概念,所以我们必须手动调整模式匹配,下面的 pat_match
函数可以从一个 Bits
中匹配出一个 Bits
和一个 Bit
,make
则是其构造对偶,可以从一个 Bits
和一个 Bit
中构造出一个 Bits
。pub fn to_snoc(self : Bits) -> Bits {
match self {
Rep(a) => Snoc(Rep(a), a)
self => self
}
}
pub fn Bits::pat_match(self : Bits) -> (Bits, Bit) {
guard self.to_snoc() is Snoc(bs, b)
(bs, b)
}
pub fn Bits::make(self : Bits, b : Bit) -> Bits {
match (self, b) {
(Rep(b), b1) if b == b1 => Rep(b)
(bs, b) => Snoc(bs, b)
}
}
.pat_match()
进行匹配时,会潜在地将一个 Rep
展开一步成为 Snoc
,这样我们就可以假装 Bits
值总是由 make
构造的。反之,使用 make
构造一个 Bits
值时,如果我们恰好将一个相同的位 b
追加到一个现有的 Rep b
上,它将什么也不做。这确保了只要我们坚持使用 make
和 pat_match
而从不直接使用 Snoc
,Bits
值将始终被规范化,使得末端的 Rep(b)
立即跟随着一个不同的位。
读者可能会认为 pat_match
中的 guard
并不安全,实际上它确实足以处理 Bits
类型的所有可能输入。然而,为了获得能停机的算法,我们通常会包含一个或多个针对 Rep
的特殊情况。Bits
和 Int
之间转换以及显示 Bits
(仅用于测试目的)的函数开始。pub fn Bits::to_bits(n : Int) -> Bits {
match n {
0 => Rep(O)
-1 => Rep(I)
n => Bits::make(Bits::to_bits(div(n, 2)), Bit::to_enum(n % 2))
}
}
pub fn Bits::from_bits(self : Bits) -> Int {
match self {
Rep(O) => 0
Rep(I) => -1
bs => {
let (bs, b) = bs.pat_match()
2 * Bits::from_bits(bs) + Bit::from_enum(b)
}
}
}
Show
的实现:impl Show for Bits with to_string(self) {
fn go {
Rep(b) => b.to_string().repeat(3) + "..."
s => {
let (bs, b) = s.pat_match()
b.to_string() + go(bs)
}
}
go(self).rev()
}
impl Show for Bits with output(self, logger) {
logger.write_string(self.to_string())
}
test "from . to == id" {
let b = Snoc(Rep(O), O) |> Snoc(I) |> Snoc(O) |> Snoc(I)
inspect!(b, content="...0000101")
let b1 = Snoc(Rep(I), O) |> Snoc(I)
inspect!(b1, content="...11101")
let b26 = Bits::to_bits(26)
inspect!(b26, content="...00011010")
let b_30 = Bits::to_bits(-30)
inspect!(b_30, content="...11100010")
inspect!(b_30.from_bits(), content="-30")
@qc.quick_check_fn!(fn(x) { Bits::from_bits(Bits::to_bits(x)) == x })
// +++ [100/0/100] Ok, passed!
}
Bits
上实现一些基本操作了。首先,递增和递减可以递归地实现如下:pub fn inc(self : Bits) -> Bits {
match self {
Rep(I) => Rep(O)
s =>
match s.pat_match() {
(bs, O) => Bits::make(bs, I)
(bs, I) => Bits::make(bs.inc(), O)
}
}
}
pub fn dec(self : Bits) -> Bits {
match self {
Rep(O) => Rep(I)
s =>
match s.pat_match() {
(bs, I) => Bits::make(bs, O)
(bs, O) => Bits::make(bs.dec(), I)
}
}
}
pub fn lsb(self : Bits) -> Bits {
match self {
Rep(O) => Rep(O)
s =>
match s.pat_match() {
(bs, O) => Bits::make(bs.lsb(), O)
(_, I) => Bits::make(Rep(O), I)
}
}
}
Rep(O)
添加了一个特殊情况,以确保 lsb
是全函数。严格来说,Rep(O)
没有最低有效位,因此定义 lsb(Rep(O)) = Rep(O)
似乎是合理的。test "lsb" {
inspect!(Bits::to_bits(26), content="...00011010")
inspect!(Bits::to_bits(26).lsb(), content="...00010")
inspect!(Bits::to_bits(24), content="...00011000")
inspect!(Bits::to_bits(24).lsb(), content="...0001000")
}
pat_match
的匹配会自动扩展较短的一个以匹配较长的那个。impl BitAnd for Bits with land(self, other) {
match (self, other) {
(Rep(x), Rep(y)) => Rep(x & y)
(xs, ys) =>
match (xs.pat_match(), ys.pat_match()) {
((bs, b), (bs1, b1)) => Bits::make(bs.lsb() & bs1.lsb(), b & b1)
}
}
}
pub fn Bits::inv(self : Bits) -> Bits {
match self {
Rep(x) => Rep(x.not())
s =>
match s.pat_match() {
(bs, b) => Bits::make(bs.inv(), b.not())
}
}
}
zipWith
来实现 land
,用 map
来实现 inv
。然而,就目前的目的而言,我们不需要这种额外的通用性。Rep
的特殊情况。
(译者注:原文中 op_add
的实现是非穷尽的,下面是修正后的版本):impl Add for Bits with op_add(self, other) {
match (self, other) {
(xs, Rep(O)) | (Rep(O), xs) => xs
(Rep(I), Rep(I)) => Bits::make(Rep(I), O)
_ =>
match (self.pat_match(), other.pat_match()) {
((xs, I), (ys, I)) => Bits::make((xs + ys).inc(), O)
((xs, x), (ys, y)) => Bits::make(xs + ys, x | y)
}
}
}
test "add" {
@qc.quick_check_fn!(fn {
(x, y) => x + y == Bits::from_bits(Bits::to_bits(x) + Bits::to_bits(y))
})
// +++ [100/0/100] Ok, passed!
}
x : Bits
,都有 $x + \text{neg } x = \text{Rep(O)}$。impl Neg for Bits with op_neg(self) {
self.inv().inc()
}
Theorem 4.1. [blog/fenwick/thm/thm2]
pat_match(bs, b)
和 make(bs, b)
记为 $bs : b$,运算名的 op_
前缀均省略。
pub fn set_to(self : Bits, idx : Int, b1 : Bit) -> Bits {
match (idx, self.pat_match()) {
(0, (bs, _)) => Bits::make(bs, b1)
(k, (bs, b)) => Bits::make(bs.set_to(k - 1, b1), b)
}
}
pub fn Bits::set(self : Bits, idx : Int) -> Bits {
self.set_to(idx, I)
}
pub fn Bits::clear(self : Bits, idx : Int) -> Bits {
self.set_to(idx, O)
}
pub fn test_helper(self : Bits, i : Int) -> Bool {
loop i, self.pat_match() {
0, (_bs, b) => b == I
k, (bs, _b) => continue k - 1, bs.pat_match()
}
}
pub fn odd(self : Bits) -> Bool {
self.test_helper(0)
}
pub fn even(self : Bits) -> Bool {
not(self.odd())
}
while
组合子,它迭代一个给定的函数,返回第一个使得谓词为假的迭代结果。pub fn shr(self : Bits) -> Bits {
self.pat_match().0
}
pub fn shl(self : Bits) -> Bits {
Bits::make(self, O)
}
pub fn while_[A](p : (A) -> Bool, f : (A) -> A, x : A) -> A {
loop x {
x => if p(x) { continue f(x) } else { x }
}
}
在推导我们的索引转换函数之前,我们必须处理一个略显棘手的事实。在传统的二叉树索引方案中,如图 9 所示,根节点的索引为 1,每个左子节点的索引是其父节点的两倍,每个右子节点的索引是其父节点的两倍加一。回想一下,在一棵删减后的线段树中,根节点和每个左子节点都是活跃的,而所有右子节点都是非活跃的。这使得根节点成为一个尴尬的特例——所有活跃节点都有偶数索引,除了索引为 1 的根节点。这使得检查我们是否处于一个活跃节点变得更加困难——仅仅查看最低有效位是不够的。 解决这个问题的一个简单方法是直接给根节点赋予索引 2,然后继续使用相同的方案标记其余节点——每个左子节点是其父节点的两倍,每个右子节点是其父节点的两倍加一。这导致了如图 14 所示的索引方式,就好像我们只是取了以 1 为根的树的左子树,并忽略了右子树。当然,这意味着大约一半的可能索引被省略了——但这不成问题,因为我们只会将这些索引作为中间步骤使用,最终会被融合掉。 图 15 展示了一棵二叉树,其中节点以两种不同的方式编号:每个节点的左侧显示其二叉树索引(根节点索引为 2)。每个节点的右侧显示其在芬威克数组中的索引,如果它有的话(非活跃节点右半部分简单地灰色显示)。下方的表格显示了从芬威克数组索引(顶行)到二叉树索引(底行)的映射。作为一个更大的例子,图 16 在更深一层的二叉树上展示了同样的事情。 我们的目标是找到一种方法来计算给定芬威克索引对应的二叉树索引,反之亦然。仔细观察图 16 中的表格,有个模式非常突出。首先,底行的所有数字都是偶数,这恰恰是因为二叉树的编号方式使得所有活跃节点都有偶数索引。其次,我们可以看到偶数 $32, 34 \ldots 46$ 按顺序出现在所有奇数位置上。这些正是树的叶节点,实际上,芬威克数组中每隔一个节点就是来自原始树的叶节点。与它们交替出现的,在偶数位置上的,是数字 $16, 8, 18, 4, \ldots$,它们对应于所有非叶节点;但这些恰好是图 15 中表格底行的二叉树索引序列——因为高度为 4 的树中的内部节点本身构成了一个高度为 3 的树,且节点以相同的顺序出现。 这些观察引出了 5.1 中所示的递推关系,用于计算存储在长度为 $2^n$ 的芬威克数组中的节点的二叉树索引序列 $b_n$:$b_0$ 就是单元素序列 $[2]$,否则 $b_n$ 就是偶数 $2^{n+1}, 2^{n+1} + 2, \ldots, 2^{n+1} + 2^n - 2$ 与 $b_{n-1}$ 交错排列的结果(译者注:代码实现中使用位运算代替乘方运算,在数学推理中,也将 我们可以检验这确实重现了 $n=4$ 时观察到的序列: 令 $s ! k$ 表示列表 $s$ 中的第 $k$ 项(从 1 开始计数),在代码中记为 有了这些准备,我们可以将芬威克到二叉树的索引转换函数定义为
$$ \text{f2b}(n, k) = b_n ! k. $$ 当然,由于 $b_n$ 的长度为 $2^n$,这个函数仅在范围 $[1, 2^n]$ 内有定义。 我们现在可以简化 $$
\begin{align*}
& \text{f2b}(n, (2 \cdot j)) \\
= & b_n ! (2 \cdot j) & \{ \text{f2b} \} \\
= & (\text{interleave}(\text{map}((2 \cdot), [2^n \ldots 2^n + 2^{n-1} - 1]), b_{n-1})) ! (2 \cdot j) & \{ \text{b} \} \\
= & b_{n-1} ! j & \{ \text{interleave-!} \text{ 引理 } \} \\
= & \text{f2b}((n - 1), j). & \{ \text{f2b} \}
\end{align*}
$$ 其中 $(2\cdot)$ 是函数 $$
\begin{align*}
& \text{f2b}(n, (2 \cdot j - 1)) \\
= & b_n ! (2 \cdot j - 1) & \{ \text{f2b} \} \\
= & (\text{interleave}(\text{map}((2 \cdot), [2^n \ldots 2^n + 2^{n-1} - 1]), b_{n-1})) ! (2 \cdot j - 1) & \{ \text{b} \} \\
= & \text{map}((2 \cdot), [2^n \ldots 2^n + 2^{n-1} - 1]) ! j & \{ \text{interleave-!} \text{ 引理 } \} \\
= & 2 \cdot (2^n + j - 1) & \{ \text{map, 代数} \} \\
= & 2^{n+1} + 2j - 2 & \{ \text{代数} \}
\end{align*}
$$ 因此,我们有 $$
\text{f2b}(n, k) =
\begin{cases}
\text{f2b}(n - 1,k / 2) & k \text{ 为偶数} \\
2^{n+1} + k - 1 & k \text{ 为奇数}
\end{cases}
$$ 注意,当 $n = 0$ 时,我们必须有 $k = 1$,因此 $\text{f2b}(0, 1) = 2^{0+1} + 1 - 1 = 2$,(译者注:原文写的 Note that when n = 0, we must have k = 1, and hence, f2b 0 1 = $2^{0+1} + 1 - 1 = 1$, as required, so this definition is valid for all n ≥ 0. 计算有误。根据推导,$f2b(0, 1) = 2 \neq 1$,等于 2 才能对应上我们新的根节点假设)所以这个定义对所有 $n \ge 0$ 都有效。现在将 $k$ 唯一地分解为 $2^a \cdot b$,其中 $b$ 是奇数。那么通过归纳我们可以看到 $$ \text{f2b} (n, 2^a \cdot b) = \text{f2b}(n - a, b) = 2^{(n-a)+1} + b - 1. $$ 换句话说,计算 容易使用 QuickCheck 验证这在范围 $[1, 2^4]$ 上产生与原始的 现在我们转向推导 $$ \text{b2f} (n, 2^c + d) = 2^{n-c+1} \cdot (d + 1). $$ 换句话说,给定输入 $2^c + d$,我们减去最高位 $2^c$,加 1,然后左移 $n - c + 1$ 次。同样,存在一个更简单的方法:我们可以先加 1(注意因为 $d \lt 2^{c-1}$,加 1 不会干扰 $2^c$ 处的位),然后左移足够多次,使最左边的位移到位置 $n + 1$,最后移除它。即: 验证: 5. Index conversion [blog/fenwick/index_conv]
Code 5.1. Recurrence for sequence of binary tree indices in a Fenwick array [blog/fenwick/interleaving]
pub fn interleave[T](sel : List[T], other : List[T]) -> List[T] {
match (sel, other) {
(Nil, _) => Nil
(Cons(x, xs), ys) => Cons(x, interleave(ys, xs))
}
}
fn b(i : Int) -> List[Int] {
match i {
0 => Cons(2, Nil)
n =>
Int::until(1 << n, (1 << n) + (1 << (n - 1))).map(fn { x => x * 2 })
|> @immut/list.from_iter
|> interleave(b(n - 1))
}
}
interleave
记为中缀运算符 $\curlyvee$)。test "b(4)" {
inspect!(
b(4),
content="@list.of([32, 16, 34, 8, 36, 18, 38, 4, 40, 20, 42, 10, 44, 22, 46, 2])",
)
}
s.nth(k)
(译者注:这和 MoonBit 标准库的 nth
存在差异,它是从 0 开始的索引)。下面的代码指出了两个关于索引和交错相互作用的简单引理,即 $(xs \curlyvee ys) ! (2 \cdot j) = ys ! j$ 和 $(xs \curlyvee ys) ! (2 \cdot j - 1) = xs ! j$(只要 $xs$ 和 $ys$ 长度相等)。// suppose xs.length() == ys.length()
interleave(xs, ys).nth(2 * j) == ys.nth(j)
interleave(xs, ys).nth(2 * j + 1) == xs.nth(j)
f2b
的定义如下。首先,对于偶数输入,我们有fn { x => 2 * x }
的简写,$[2^n \ldots 2^n + 2^{n-1} - 1]$ 是闭区间 $[2^n, 2^n + 2^{n-1} - 1]$ 的列表。对于奇数输入有:f2b
包括只要输入是偶数就重复除以 2(即右位移),然后最终减 1 并加上一个 2 的幂。然而,要知道最后要加哪个 2 的幂,取决于我们移动了多少次。思考这个问题的一个更好的方法是在开始时加上 $2^{n+1}$,然后让它随着其他所有位一起移动。(译者注:原文在这里说的并不是很清楚,在这里补充一下我的解释:我们想要求解 $X = 2^{n-a+1} + b$,可以先乘上一个因子 $2^a$ 得到 $2^a X = 2^{n+1} + 2^ab$,这样表达式中便有一个 $k = 2^a \cdot b$ 存在,它是函数的参数。这样一来等式的右边即是 $k$ 加上一个 $2^{n+1}$,而 $X = \dfrac{2^{n+1} + 2^a \cdot b}{2^a}$,除法在这里可以使用位移操作实现)因此,我们使用我们的 Bits
DSL 得到 f2b
的最终定义。单独定义 shift
函数将使我们的一些证明更加紧凑。pub fn shift(n : Int, se : Bits) -> Bits {
while_(even, shr, se.set(n))
}
pub fn f2b(n : Int, se : Bits) -> Bits {
shift(n + 1, se) |> dec
}
f2b(4, k)
相同的结果。b2f(n, _)
,它将从二叉树索引转换回芬威克索引。b2f(n, _)
应该是 f2b(n, _)
的左逆,也就是说,对于任何 $k \in [1, 2^n]$,我们应该有 $\text{b2f}(n, \text{f2b}(n, k)) = k$。
如果 $k$ 是 f2b
的一个输入,我们有 $k = 2^a \cdot b \le 2^n$,因此 $b-1 \lt 2^{n-a}$。故给定输出 $\text{f2b}(n, k) = m = 2^{n-a+1} + b - 1$,
$m$ 的最高位是 $2^{n-a+1}$,其余位代表 $b - 1$。所以,一般地,给定某个作为 f2b(n, _)
输出的 $m$,我们可以唯一地将其写为 $m = 2^c + d$,其中 $d \lt 2^{c-1}$;那么pub fn unshift(n : Int, se : Bits) -> Bits {
while_(fn { x => not(x.test_helper(n)) }, shl, se) |> Bits::clear(_, n)
}
pub fn b2f(n : Int, se : Bits) -> Bits {
unshift(n + 1, se.inc())
}
test "id" {
let f = fn {
size =>
Int::until(1, 1 << size)
.map(fn {
x =>
(f2b(size, x |> Bits::to_bits) |> b2f(size, _) |> Bits::from_bits) ==
x
})
.all(fn { x => x })
}
assert_true!(f(4))
assert_true!(f(5))
}
我们现在终于可以推导出在芬威克数组索引上移动所需的运算了,方法是从二叉索引树上的操作开始,并通过与芬威克索引的转换进行变换。首先,为了融合掉由此产生的转换,我们需要几个引理。 对于所有为奇数(即以 I 结尾)的 $bs : \text{Bits}$, 证明 两者均由定义直接可得。 以下两条定理对所有 Bits 值均成立: 证明 简单的 最后,我们需要一个关于将零位移入和移出值右侧的引理。 对于所有 $0 \lt x \lt 2^{n+2}$, $$ \text{while}(\text{not}(\text{test}(\cdot, n + 1)), \text{shl}, \text{while}(\text{even}, \text{shr}, x)) = \text{while}(\text{not}(\text{test}(\cdot, n + 1)), \text{shl}, x) $$ (注:$\text{test}$ 在代码中是 证明 直观上,这表明如果我们先移出所有零位,然后左移直到第 $n+1$ 位被设置,这可以通过完全忘记右移来获得相同的结果;移出零位然后再将它们移回来应该是恒等操作。 形式上,证明通过对 $x$ 进行归纳。如果 $x = xs : \text{I}$ 为奇数,等式立即成立,因为 $\text{while}(\text{even}, \text{shr}, x) = x$。否则,如果 $x = xs : \text{O}$,在左侧,O 会被 $\text{shr}$ 立即丢弃,而在右侧,$xs : \text{O} = \text{shl}(xs)$,并且由于 $xs \lt 2^{n+1}$,额外的 $\text{shl}$ 可以被吸收到 $\text{while}$ 循环中。剩下的就是归纳假设。 有了这些引理在手,让我们看看如何在芬威克数组中移动以实现 让我们看看如何推导出这种行为。
要在二叉索引方案下找到一个节点的最近活跃父节点,我们首先向上移动到直接父节点(通过将索引除以二,即执行一次右位移);然后只要当前节点是右子节点(即索引为奇数),就继续向上移动到下一个直接父节点。这产生了如下定义: 这就是为什么我们使用了根节点索引为 2 的略显奇怪的索引方案——否则这个定义对于任何活跃父节点是根节点的节点都不起作用! 现在,要推导出芬威克索引上的相应操作,我们通过与芬威克索引的转换进行共轭,计算如下。为了使计算更易读,每一步中被重写的部分都用下划线标出。 $$
\begin{align*}
& \underline{\text{b2f}(n, \text{activeParentBinary}(\text{f2b}(n, \cdot)))} \\ & \{ \text{展开定义} \} \\
= & \text{unshift}(n + 1, \underline{\text{inc}(\text{while}(\text{odd}, \text{shr}}, \text{shr}(\text{dec}(\text{shift}(n + 1, \cdot)))))) \\ & \{ \text{引理 (while-inc-dec)} \} \\
= & \text{unshift}(n + 1, \text{while}(\text{even}, \text{shr}, \text{inc}(\underline{\text{shr}(\text{dec}}(\text{shift}(n + 1, \cdot)))))) \\ & \{ \text{引理 (shr-inc-dec); shift(n+1, x) 总是奇数} \} \\
= & \text{unshift}(n + 1, \text{while}(\text{even}, \text{shr}, \underline{\text{inc}(\text{shr}}(\text{shift}(n + 1, \text{dec}(\cdot)))))) \\ & \{ \text{引理 (shr-inc-dec)} \} \\
= & \text{unshift}(n + 1, \underline{\text{while}(\text{even}, \text{shr}, \text{shr}}(\text{inc}(\text{shift}(n + 1, \cdot))))) \\ & \{ \text{while(even, shr, shr(x)) = while(even, shr, x) 在偶数输入上} \} \\
= & \underline{\text{unshift}(n + 1}, \text{while}(\text{even}, \text{shr}, \text{inc}(\text{shift}(n + 1, \cdot)))) \\ & \{ \text{unshift} \} \\
= & \text{clear}(n + 1, \underline{\text{while}(\text{not}(\text{test}(n + 1, \cdot)), \text{shl}, \text{while}(\text{even}, \text{shr}}, \text{inc}(\underline{\text{shift}(n + 1, \cdot)})))) \\ & \{ \text{引理 (shl-shr); shift} \} \\
= & \text{clear}(n + 1, \text{while}(\text{not}(\text{test}(n + 1, \cdot)), \text{shl}, \text{inc}(\text{while}(\text{even}, \text{shr}, \text{set}(n + 1, \cdot))))) \\
\end{align*}
$$ 在最后一步中,由于输入 $x$ 满足 $x \le 2^n$,我们有 $\text{inc}(\text{shift}(n + 1, x)) \lt 2^{n+2}$,因此引理 shl-shr 适用。 从右到左阅读,我们刚刚计算出的流水线执行以下步骤: 直观地说,这看起来很像加上 LSB!一般来说,要找到 LSB,必须移过连续的 0 位,直到找到第一个 1;问题是如何跟踪移过了多少个 0 位。 为了使这一点更形式化,我们首先定义一个辅助函数 对于所有 $x : \text{Bits}$, $\text{add}(x, \text{lsb}(x)) = \text{atLSB}(\text{inc}, x)$ 且 $\text{subtract}(x, \text{lsb}(x)) = \text{atLSB}(\text{dec}, x)$。 (注:代码中 $\text{atLSB}$ 的定义为 证明 对 $x$ 进行简单的归纳即可。 我们可以将“带哨兵的移位”方案与 令 $n \ge 1$ 且令 $f : \text{Bits} \to \text{Bits}$ 是一个函数,使得 那么对于所有 $0 \lt x \lt 2^n$,
$$ \text{unshift}(n + 1, f(\text{shift}(n + 1, x))) = \text{atLSB}(f, x). $$ 证明相当繁琐但并非特别具有启发性,因此我们省略它(一个包含完整证明的扩展版本可以在作者的网站上找到)。然而,我们确实注意到 加上 LSB 是沿芬威克索引树向上移动到最近活跃父节点的正确方法,即
$$\begin{align*} & \text{activeParentFenwick}(x) \\ &= \text{b2f}(n, \text{activeParentBinary}(\text{f2b}(n, x))) \\ &= \text{add}(x, \text{lsb}(x)) \end{align*} $$
在范围 $[1, 2^n)$ 上处处成立。(我们排除了 $2^n$,因为它对应于芬威克索引方案下的树根。) 证明 $$
\begin{aligned}
& \text{b2f}(n, \text{activeParentBinary}(\text{f2b}'(n, x))) \\
= & \text{unshift}(n + 1, \text{inc}(\text{shift}(n + 1, x))) & \{ \text{先前的计算} \} \\
= & \text{atLSB}(\text{dec}, x) & \{ \text{引理 (sentinel)} \} \\
= & \text{add}(x, \text{lsb}(x)) & \{ \text{引理 (add-lsb)} \}
\end{aligned}
$$ 我们可以进行类似的过程来推导前缀查询的实现(据称它涉及减去 LSB)。同样,如果我们想计算 $[1, j]$ 的和,我们可以从芬威克数组中的索引 $j$ 开始,它存储了结束于 $j$ 的唯一段的和。如果索引 $j$ 处的节点存储了段 $[i, j]$,我们接下来需要找到存储结束于 $i-1$ 的段的唯一节点。我们可以重复这样做,一路累加段的和。 从图 20 中寻找灵感,我们可以看到我们想要做的是找到我们最近的非活跃父节点的左兄弟,也就是说,我们向上走直到找到第一个作为右子节点的祖先,然后移动到它的左兄弟。在二叉索引方案下,这可以简单地实现为: 减去 LSB 是将芬威克树向上移动到覆盖当前段的之前段的活动节点的正确方法、
$$\begin{align*} & \text{prevSegmentFenwick}(x) \\ &= \text{b2f}(n, \text{prevSegmentBinary}(\text{f2b}(n, x))) \\ &= \text{subtract}(x, \text{lsb}(x)) \end{align*} $$
在范围 $[1, 2^n)$ 上处处成立。 证明 $$
\begin{aligned}
& \text{b2f}(n, \text{prevSegmentBinary}(\text{f2b}(n, x))) \\
& \{ \text{展开定义} \} \\
= & \text{unshift}(n + 1, \underline{\text{inc}(\text{dec}}(\text{while}(\text{even}, \text{shr}, \text{dec}(\text{shift}(n + 1, x)))))) \\
& \{ \text{引理 (while-inc-dec)} \} \\
= & \underline{\text{unshift}(n + 1}, \text{while}(\text{even}, \text{shr}, \text{dec}(\text{shift}(n + 1, x)))) \\
& \{ \text{unshift} \} \\
= & \text{clear}(n + 1, \underline{\text{while}(\text{not}(\text{test}(n + 1, \cdot)), \text{shl} , \text{while}(\text{even}, \text{shr}}, \text{dec}(\text{shift}(n + 1, x)))))) \\
& \{ \text{引理 (shl-shr)} \} \\
= & \underline{\text{clear}(n + 1, \text{while}(\text{not}(\text{test}(n + 1, \cdot)), \text{shl}}, \text{dec}(\text{shift}(n + 1, x)))) \\
& \{ \text{unshift} \} \\
= & \text{unshift}(n + 1, \text{dec}(\text{shift}(n + 1, x))) \\
& \{ \text{引理 (sentinel)} \} \\
= & \text{atLSB}(\text{dec}, x) \\
& \{ \text{引理 (add-lsb)} \} \\
= & \text{subtract}(x, \text{lsb}(x)) \\
\end{aligned}
$$ 6. Deriving Fenwick operations [blog/fenwick/deriving]
Theorem 6.1. shr-inc-dec [blog/fenwick/thm/thm3]
Theorem 6.2. while-inc-dec [blog/fenwick/thm/thm4]
Bits
归纳证明。例如,对于 inc
的情况,两侧的函数都会丢弃连续的 1 位,然后将第一个 0 位翻转为 1。Theorem 6.3. shl-shr [blog/fenwick/thm/thm5]
test_helper
,$\cdot$ 在此处是匿名函数简写: fn { x => not(test_helper(x, n + 1)) }
)update
和 query
;我们先从 update
开始。在实现 update
操作时,我们需要从一个叶节点开始,沿着路径向上到达根节点,更新沿途所有活跃的节点。实际上,对于任何给定的叶节点,其最近的活跃父节点恰好是存储在过去对应于该叶节点的槽中的节点(参见图 13)。因此,要更新索引 $i$,我们只需从芬威克数组中的索引 $i$ 开始,然后重复找到最近的活跃父节点,边走边更新。回想一下,用于 update
的命令式代码就是这样做的,通过在每一步加上当前索引的 LSB 来找到最近的活跃父节点:pub fn update(self : FenwickTree, i : Int, delta : Int) -> Unit {
let mut i = i
while i < self._.length() {
self._[i] += delta
i += FenwickTree::lsb(i)
}
}
pub fn active_parent_binary(self : Bits) -> Bits {
while_(odd, shr, self.shr())
}
lsb
函数本身通过递归栈来跟踪;找到第一个 1 位后,递归栈展开并将所有递归经过的 0 位重新追加回去。上述流水线代表了一种替代方法:将位 $n+1$ 设置为“哨兵”来跟踪我们移动了多少;右移直到第一个 1 确实在个位上,此时我们加 1;然后通过左移将所有 0 位移回,直到哨兵位回到 $n+1$ 的位置。这个过程的一个例子如图 19 所示。当然,这只适用于值足够小,以至于哨兵位在整个操作过程中不会受到干扰的情况。at_lsb
,它执行一个“在 LSB 处”的操作,即它移出 0 位直到找到一个 1,应用给定的函数,然后恢复 0 位。fn at_lsb(self : Bits, f : (Bits) -> Bits) -> Bits {
match self {
Rep(O) => Rep(O)
s =>
match s.pat_match() {
(bs, O) => bs.at_lsb(f) |> Bits::make(O)
(bs, I) => Bits::make(bs, I) |> f
}
}
}
Theorem 6.4. add-lsb [blog/fenwick/thm/thm6]
at_lsb
)atLSB
形式化地关联起来,通过以下(相当复杂的)引理:Theorem 6.5. sentinel [blog/fenwick/thm/thm7]
inc
和 dec
都符合 $f$ 的标准:只要 $n \ge 1$,对某个 $0 \lt x \lt 2^n$ 进行递增或递减不会影响第 $(n+1)$ 位,并且对一个小于 $2^n + 2^{n-1}$ 的数进行递增或递减的结果将是一个小于 $2^{n+1}$ 的数。我们现在可以将所有部分组合起来,证明在每一步加上 LSB 是实现 update
的正确方法。Theorem 6.6. [blog/fenwick/thm/thm8]
fn prev_segment_binary(self : Bits) -> Bits {
while_(even, shr, self).dec()
}
Theorem 6.7. [blog/fenwick/thm/thm9]
据我所知,从历史上看,芬威克树并非如本文所呈现的那般,作为线段树的一种优化而被开发出来。本文所叙仅仅是一种虚构的思想别史(但希望能有所启发),旨在彰显函数式思维、领域特定语言以及等式推理在探索不同结构与算法之间关系时的力量。作为未来工作,探索前文提及的线段树的某些泛化推广,考察是否能从中推导出支持额外操作的类芬威克结构,将不失为一桩趣事。 7. Conclusion [blog/fenwick/conclusion]