使用 O (1) 时间复杂度计算比 x 大的最小的 2 的整数次幂
提到快速计算比给定整数 x
大(或相等)的最小的 2 ** n
的方法,一般来讲,我们都会这么算:
1 | function minPowerOf2(x) { |
或者有的数学功底不扎实的人可能就会这么写:
1 | function minPowerOf2(x) { |
看看这个 for
,看看这个倍增的变量,时间复杂度显然是 的。但是它实在是太慢了!有没有 的办法呢?
答案是:有。
1 | const epsR = 2 ** 53; // eps^-1 |
跑跑看:
1 | > minPowerOf2(1000) |
有点反直觉,但是确实工作了。这是什么原理呢?
背景知识介绍#
首先我们一眼可以看出 epsR
是 Number.MAX_SAFE_INTEGER + 1
。这个数字是 IEEE 754 标准下双精度浮点数(64 bits)所能精确表示的最大整数:
[63] | [62:52] | [51:0] |
---|---|---|
S(符号位) | E(阶码位) | M(小数位) |
0 表示正,1 表示负 | 1~2046 | 任意 |
其中,指数 的所有比特位为全 0 或全 1 时有特殊定义:
0b00000000000
=0x000
用于表示 signed zero (当 时)和 subnormals (当 时);0b11111111111
=0x7ff
用于表示 ∞ (当 时)和 NaNs (当 时)。
不过,以上两种情况都和本算法关系不是特别大,可以不用考虑。
二进制和十进制换算#
对于二进制小数数 ,我们可以将其按照如下规则转换为十进制数 :
- 举例,不能精确表出的:,是一个二进制无限循环小数。
- 不能表出的:根据 信息论, 能表示的数字最多只有 个,无论是整数还是小数。
那么我们可以知道,52 个小数位最多能表示出 个数,而 到 (包含首尾)总共就包括 个整数。这也是为什么常量 MAX_SAFE_INTEGER
的值是 2 ** 53 - 1
:比 MAX_SAFE_INTEGER
小的正整数可以被表出,但是不能保证每一个比它大的正整数都能被表出。
这里有一个问题:为什么能表示出的 个数正好与 到 上的每一个整数一一对应呢?而不是对应 x.5
或是别的什么小数?其实是因为当双精度浮点数 落在 到 之间时, 的阶码位 是固定的,为 。而小数位 的精度(后文统称最小粒度)是 ,所以 的最小粒度正好是 。因此在这个范围内每个能被 1 整除的数字都可以与某一个双精度浮点数 一一对应。同理,随着指数从 52 逐渐变小, 的最小粒度也会逐渐变小为 (对应指数为 51)、(对应 50)…… 最终,在 0 和 之间达到 ,也即 JS 中的 Number.EPSILON
。综上所述,比 MAX_SAFE_INTEGER
小的正整数是绝对可以被精确表出的。
原理分析#
TL;DR 其实就相当于找双精度浮点数的最高位
为了方便理解这些代码在干什么,我找到了 GitHub Gist 里的现成的可以显示双精度浮点数的二进制表示的轮子,并根据 gist 评论区指正内容稍作修改:
11e971a7
来看一个例子:x = 3
。
首先,使用上面提到的代码将有用的信息都显示在屏幕上:
Expression | Value | Binary (IEEE754 Double) |
---|---|---|
x | 3 | 0 10000000000 100..(45x0)..0000 |
epsR | 9007199254740992 | 0 10000110100 000..(45x0)..0000 |
epsR * x | 27021597764222976 | 0 10000110101 100..(45x0)..0000 |
epsR * x + x | 27021597764222980 | 0 10000110101 100..(45x0)..0001 |
deviation | 4 | 0 10000000001 000..(45x0)..0000 |
显然,相比于 epsR
,epsR * x
在 (阶码位)和 (小数位)上都有变化: - 增大了 Math.floor(Math.log2(x))
- 处理了 增加完之后剩下的那部分(epsR * x - 2 ** Math.floor(Math.log2(x))
)
这两步对于常规的双精度浮点数都是完全 OK 的。然而 epsR
实在是太大,导致在处理 epsR * x + x
时,到第二步出现了问题:
对于 这种情况,到第四步之后,想要表示 27021597764222980
, 的全部和 的第一位已经表示了 epsR * x
,所以 的剩下位应该表示 x
。然而 的剩下位不可能精确表示出 x
了:此时,阶码位 已经全部确定,使得十进制表达式中第二项固定为 ,要想表示出 3 则需要 的最小粒度变成 。然而因为 总共只有 52 位,所以只能苦哈哈地舍弃一部分精度,转换成比 epsR * x + x
大的、能精确表示的最接近 epsR * x + x
的数字。
—— 停一下!为什么到 epsR * x + x
时才出现误差?epsR * x
也超过 Number.MAX_SAFE_INTEGER
了啊?
—— 这是因为当 x 在 也即 范围内时,Math.floor(Math.log2(x))
的取值范围为 ,所以 epsR * x - 2 ** Math.floor(Math.log2(x))
都可以被小数位 正确处理。
具体来讲,前文提到了 且根据十进制表示可知 ,所以对于任何比 大的 的正整数倍数 ,必然可以构造 ,其中 —— 这样, 的最小粒度 乘以新增加的 ,可以得出 的最小粒度必然小于等于 1(且是 2 的整数次幂),所以对于每一个这样的 都存在一个双精度浮点数与之一一对应。
但是到了 epsR * x + x
时,这个性质不再被满足,就有可能出现误差。为了更好地展示这一点,我们可以观察一下 epsR * x + x
的 (小数位)的最后几位是在何时被抛弃的。
理想情况下,不丢失精度的表示方式:
1 | ...00 1 |
实际的表达方式:
1 | ...01 |
那么,只剩下最后一块拼图,就是为什么舍弃精度后 epsR * x + x
会转换成比 epsR * x + x
大的最小的双精度浮点数,而不是比 epsR * x + x
小的最大的双精度浮点数?(注:这里讨论的都是 x
不是 2 的正整数次幂的情况)
其实答案非常简单。不难看出此时阶码位 是由 epsR * x
确定的,因为 + x
部分太小了(),不足以影响到 。此时,存在 ,故此时该双精度浮点数的最小粒度是 也即「比 x 大的最小的 2 的正整数次幂」。所以 x
必然大于最小粒度的一半。
上述文字论述,可以转化为更直观的图表。在数轴上,存在如下几个数:
注意,
epsR * x + H/2
不是双精度浮点数,图示仅仅是为了方便参考大小关系。
那么自然而然,舍入时就会舍入到「比 epsR * x + x
大的最小的双精度浮点数」上了,而不会舍入到更远的「比 epsR * x + x
小的最大的双精度浮点数」。
所以我们可以从这个例子出发,进行一些普适的论证:刚刚说到的这个新的数字就会等于 epsR * x
加上一个比 x
大的最小的能够满足 的最小粒度的数字了(我们就叫它数字 吧);随着 中的 不断变大,能满足 的最小粒度的 也会不断变大。又因为最终目标是 ,自然要想取到 的最小值,等价于取到 的最小值,也就是 。
所以,得到结论: epsR * x + x
大的最接近它的双精度浮点数」里分离出来,自然就再减去一个 epsR * x
就可以了。因为「比 epsR * x + x
大的最接近它的双精度浮点数」和 epsR * x
都可以被双精度浮点数精确表出,所以这一步不会损失精度。
最后就得到了 ,也就是比 x
大的最小的 2 的正整数次幂。
补充说明:当 x
本身就是 2 的正整数次幂时,上述所有计算都不会有精度损失,所以最后返回的还是 x
。因此本文标题为了严谨,给出的说法是「比给定整数 x
大(或相等)的最小的 2 的正整数次幂」。