[Data Structure]Binary Indexed Tree/Fenwick Tree

本文主要参考topcoder讲解BIT(Binary Index Tree)的文章，图片均来自链接的文章。

Binary Index Tree/Fenwick Tree是一种应用在区间查询的数据结构。其基本的思想是任何整数都可以分解成2的指数为基的线性组合；以此类推，任何一个区间也可以分解成不同的长度为2的指数的区间的并集。那么BIT是如何运作的呢？我们以数组arr(长度为N)为例，所求的为区间和，对应的BIT也是一个长度为N的数组tree，对于index = i，tree[i]存的是：

• sum(arr[i - 2^j + 1]: arr[i])，其中j为i的二进制表示从右到左第一个1的位置(从0开始)
• 换句话说，如果i的二进制表示为a1b，其中1为从右到左第一个1，b全为0。那么tree[i]存的是从arr[a0b + 1]到arr[a1b]的和
• 以12为例，其二进制为1100；去掉最低位的1为1000，十进制为8。那么tree[12] = arr[12] + arr[11] + arr[10] + arr[9]

如图所示，一个长度为16的数组对应的BIT数组每个索引所存的区间和如下图所示。注意BIT的数组下标从1开始，所以长度比原数组多1：

如果我们用多叉树的形式表示，如下图所示：

这个树具有以下特点：

• 每个节点所表示的区间依然如上述，即节点i表示[i - 2^j + 1, i]的区间, j为LSB(lowest significant bit)所在的位置
• 每个节点n的所有子节点表示的区间的是连续相邻的区间，并且这些区间的并集为和n表示的区间相等
• 每个节点n的子节点为n - 2^(j - 1), n - 2^(j - 2)..., n - 2^0, j为LSB所在的位置
• 每个节点n的父节点为n + 2^j, j为LSB所在的位置

对应更新和查询的操作为：

• 查询区间[0, i]的和: 每次查询tree[i]并且抹掉i的最后一个1直到1等于0。累加所有对应的tree[i]的值。
• 拿11(1011)举例，抹掉最后一个1之后为10(1010)。第一次我们查询tree[11] = arr[11] 
• 抹掉10(1010)最后一个1之后为8(1000)。第二次查询tree[10] = arr[10] + arr[9]
• 抹掉8(1000)最后一个1之后为0。最后一次查询tree[8] = tree[8] + tree[7] + ... + tree[1]
• 所以tree[11] + tree[10] + tree[8] = arr[11] + arr[10] + ... + arr[1]
• 查询区间[i, j]的和 = 区间和[0, i] - 区间和[0, j - 1]
• 单点更新arr[i]的值：我们要更新所有包含i的区间。所以我们每一次更新 i += 1 << j，j为i的最低位的1。
• 假设输入数组长度为16
• 拿i = 10(1010)举例，第一次我们更新tree[10] = arr[10] + arr[9]
• 更新i之后i = 10 + 2 = 12(1100)。第二次更新tree[12] = arr[12] + ... + arr[9]
• 更新i之后i = 12 + 4 = 16(10000)。第三次更新tree[16] = arr[16] + ... + arr[1]
• i之后会超过16，停止

实现的话有一个小技巧，我们可以通过i & -i来取i的最低为1的位。具体的原理可以参考这篇文章。实现起来代码非常的简洁。假设区间长度为N，单点更新和区间查询的复杂度为O(log N)。

假设我们要实现如下功能：

• 更新处于[0, N]区间的数, [0, N]区间初始为0
• 统计[i, j]的区间的区间和，其中i <= j <= N

那么实现代码如下：

2D Binary Index Tree

2D BIT和2D Segment Tree的思想是类似的。对于输入为m x n的数组，对应的2D BIT就是一个(m + 1) x (n + 1)的数组。每一行就是对应输入数组那一行的1D BIT。更新和查询依然是很类似的，我们要读/写对应行列的cell的数据。以下图为例，输入数组是一个16 * 8的2D矩阵，如果我们要update input矩阵的(4, 2)，也就是BIT中的(5, 3)：

• 要update的行为(5, 6, 8 16)
• 要update的列为(3, 4, 8)
• 要update的cell如下图蓝色所示

时间复杂度的话单点更新和查询就变为O(log m) * O(log n)。具体的实现可以参考这道求2D区间和的题目。

BIT与RMQ

我们可以看出，区间和的运算是一个可逆的过程。比如sum(0, 4) = 5, sum(5, 7) = 6, 那么sum(0, 7) = sum(0, 4) + sum(5, 7) = 11。同理我们可以反过来推出sum(5,7) = sum(0, 7) - sum(0, 4) = 11 - 5 = 6。但是并不是所有的区间查询都拥有这样的性质。一个很明显的例子，RMQ(Range Minimum Query)问题，我们可以推导出min(0, 7) = min(min(0, 4), min(5, 7))，但是这个过程不是可逆的，也就是说已知min(0, 7)和min(0, 4)的情况下，我么不一定能推导出min(5, 7)的值。我们之前实现区间和BIT的时候利用了区间和的这种性质，比如查询的时候我们查询两个从0开始的区间和然后相减。更新的时候我们直接在每个父节点的区间和上累加。但是对于RMQ问题，我们都无法进行这种操作，查询不必说，更新的话如果我们把arr中位于i位置的值oldVal更新为val，BIT的节点对应的值如何改变有以下几种情况：

1. 如果oldVal > val，那么我们从i对应的节点向上一直更新，即currNode.val = min(currNode.val, val)
2. 如果oldVal == val，我们不需要做任何事
3. 如果oldVal < val，情况会比较复杂：
1. 首先如果currNode.val != oldVal，代表oldVal不是当前区间的最小值，所以对于当前节点和父节点我们不需要做任何事
2. 如果currNode.val != oldVal，代表oldVal是当前区间的最小值，我们把oldVal更新为更大的val，那么当前区间的最小值是可能变的。那么我们需要去当前节点的所有子节点查询新的最小值，再更新当前区间的最小值。同理，对父节点我们也要做同样的事情，知道我们到达根节点或者发现不需要继续更新了。

更新的思路我们清楚了，但是如果按照Region Sum的方式存，我们还是没有办法解决查询的问题。为了解决这个问题，我们需要改变存储的方式，用两个BIT来存数据，具体的详情可以参照这篇paper。我们在这里简单介绍：

• 第一个BIT为我们之前介绍的BIT，节点i存的是[i - 2^j + 1, i]区间的最小值，j为LSB所在的位置
• 第二个BIT为反向的BIT，节点i存的是[i, i + 2^j - 1]区间的最小值，j为LSB所在的位置。节点n的父亲为i - (i & i), 儿子为n + 2^(j - 1), n + 2^(j - 2)..., n + 2^0, j为LSB所在的位置
• 此外我们还需要一个数组val来存原始的输入数组，val[i] = arr[i]。我们总共需要这三个数组存储。

两个BIT如下图所示(图片来自上面的论文，忽略每个节点存的值，我们只看结构)，我们实现的时候BIT2还是当做从1开始，所以第二个图会有些许不同:

对于区间[s, e]的查询，输入查询区间为[s - 1, e - 1]，因为我们存BIT从0开始:

1. 在BIT1中我们从s对应的节点开始向上爬，保证s不大于e。
2. 在BIT2中我们从e对应的节点开始向上爬，保证e不小于s。
3. 在两次攀爬的过程中我们最终停在的节点的index是一样的。证明如下: 我们可以把s和e写成二进制表达式，并且假设两个数第一个不一样的位从p + 1(从左到右)开始:
• s = a1a2a3...ap0bp+2...bn
• e = a1a2a3...ap1cp+2...cn
• 我们每次更新s += s & -s，我们最终会到达k = a1a2a3....ap10...0
• 我们每次更新e -= s & -s，我们最终也会到达k
4. 在1中我们经过的所有节点的index，我们去BIT2中取index对应节点对应区间的最小值
5. 在2中我们经过的所有节点的index，我们去BIT1中取index对应节点对应区间的最小值
6. 4和5对应的所有区间在加上val[k]的最小值，就是查询区间的最小值

以查询区间[4, 12]距离，对应BIT的index的查询区间为[5, 13]:

1. 我们首先在BIT1中从5开始爬，会经过5, 6最后停在8，在BIT2中对应的区间为[5, 5], [6, 7]
2. 同理在BIT2中从13开始爬，会经过13, 12最后停在8，在BIT1中对应的区间为[13, 13], [9, 12]
3. 那么val[8]加上1中[5, 5], [6, 7]和2中[13, 13], [9, 12]的区间正好对应了[5, 13]的区间，那么查询的结果就是[5, 13]的最小值

查询的时间复杂度为O(log n)，更新的话因为对于每个节点我们可能需要查一下所有的子节点，所以时间复杂度为O(log^2 n)，代码如下：