提到快速计算比给定整数 x 大(或相等)的最小的 2 ** n 的方法,一般来讲,我们都会这么算:

1
2
3
function minPowerOf2(x) {
return Math.pow(2, Math.ceil(Math.log(x) / Math.log(2)));
}

或者有的数学功底不扎实的人可能就会这么写:

1
2
3
4
5
6
7
8
function minPowerOf2(x) {
let power = 1;
for(;;) {
if(power >= x)
return power;
power *= 2;
}
}

看看这个 for,看看这个倍增的变量,时间复杂度显然是 Θ(logn) 的。但是它实在是太慢了!有没有 Θ(1) 的办法呢?

答案是:有。

1
2
3
4
5
6
7
const epsR = 2 ** 53;  // eps^-1
function minPowerOf2(x) {
const base = epsR * x;
const deviation = Math.abs(base + x - base);
if(deviation === 0) return Math.abs(x); // 也可以自己实现 abs() 函数
return deviation;
}

跑跑看:

1
2
> minPowerOf2(1000)
1024

有点反直觉,但是确实工作了。这是什么原理呢?

背景知识介绍#

首先我们一眼可以看出 epsRNumber.MAX_SAFE_INTEGER + 1。这个数字是 IEEE 754 标准下双精度浮点数(64 bits)所能精确表示的最大整数:

[63] [62:52] [51:0]
S(符号位) E(阶码位) M(小数位)
0 表示正,1 表示负 1~2046 任意

其中,指数 E 的所有比特位为全 0 或全 1 时有特殊定义:

  • 0b00000000000 = 0x000 用于表示 signed zero (当 F=0 时)和 subnormals (当 F0 时);
  • 0b11111111111 = 0x7ff 用于表示 (当 F=0 时)和 NaNs (当 F0 时)。

不过,以上两种情况都和本算法关系不是特别大,可以不用考虑。

二进制和十进制换算#

对于二进制小数数 101.01b,我们可以将其按照如下规则转换为十进制数 5.25

5.25=1×22+0×21+1×20+0×21+1×22
所以我们可以将双精度浮点数的十进制形式表示为:
(1)S×2E1023×(1+n=152(Mn×2n))
从双精度浮点数的二进制定义及十进制表示里,我们能发现当阶码位 E 越小时,所能精确表示的数字在数轴上就越密;随着 E 逐渐增大,能精确表示的数字在数轴上分布越来越稀疏。这是因为 (1+n=152(Mn×2n)) 的取值范围恒为 [0,2) 。事实上,在双精度浮点数所能表示的最小数字和最大数字之间,超过 99%小数点后不超过 16 位的小数无法被双精度浮点数表出或精确表出

  • 举例,不能精确表出的:0.1D=0.00˙011˙B,是一个二进制无限循环小数。
  • 不能表出的:根据 信息论Nbits 能表示的数字最多只有 2N 个,无论是整数还是小数。

那么我们可以知道,52 个小数位最多能表示出 252 个数,而 2522531 (包含首尾)总共就包括 252 个整数。这也是为什么常量 MAX_SAFE_INTEGER 的值是 2 ** 53 - 1:比 MAX_SAFE_INTEGER 小的正整数可以被表出,但是不能保证每一个比它大的正整数都能被表出。

这里有一个问题:为什么能表示出的 252 个数正好与 2522531 上的每一个整数一一对应呢?而不是对应 x.5 或是别的什么小数?其实是因为当双精度浮点数 P 落在 2522531 之间时,P 的阶码位 E 是固定的,为 52+1023=1076。而小数位 M 的精度(后文统称最小粒度)是 252,所以 P 的最小粒度正好是 2E1023×252=1。因此在这个范围内每个能被 1 整除的数字都可以与某一个双精度浮点数 P 一一对应。同理,随着指数从 52 逐渐变小,P 的最小粒度也会逐渐变小为 12(对应指数为 51)、14(对应 50)…… 最终,在 0 和 20=1 之间达到 1252,也即 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

显然,相比于 epsRepsR * xE(阶码位)和 M(小数位)上都有变化: - E 增大了 Math.floor(Math.log2(x)) - M 处理了 E 增加完之后剩下的那部分(epsR * x - 2 ** Math.floor(Math.log2(x))

这两步对于常规的双精度浮点数都是完全 OK 的。然而 epsR 实在是太大,导致在处理 epsR * x + x 时,到第二步出现了问题:

对于 x=3 这种情况,到第四步之后,想要表示 27021597764222980E 的全部和 M 的第一位已经表示了 epsR * x,所以 M 的剩下位应该表示 x。然而 M 的剩下位不可能精确表示出 x 了:此时,阶码位 E 已经全部确定,使得十进制表达式中第二项固定为 210771023=254,要想表示出 3 则需要 M 的最小粒度变成 254。然而因为 M 总共只有 52 位,所以只能苦哈哈地舍弃一部分精度,转换成比 epsR * x + x 大的、能精确表示的最接近 epsR * x + x 的数字。

—— 停一下!为什么到 epsR * x + x 时才出现误差?epsR * x 也超过 Number.MAX_SAFE_INTEGER 了啊?

—— 这是因为当 x 在 [0,epsR1] 也即 [0,2531] 范围内时,Math.floor(Math.log2(x)) 的取值范围为 [0,52],所以 epsR * x - 2 ** Math.floor(Math.log2(x)) 都可以被小数位 M 正确处理。

具体来讲,前文提到了 M[0,2) 且根据十进制表示可知 M0mod252,所以对于任何比 epsR 大的 epsR 的正整数倍数 P,必然可以构造 P=(1)0×253×A×M,其中 A=Math.floor(Math.log2(x)) —— 这样,M 的最小粒度 252 乘以新增加的 2A(A[0,52]),可以得出 P 的最小粒度必然小于等于 1(且是 2 的整数次幂),所以对于每一个这样的 P 都存在一个双精度浮点数与之一一对应

但是到了 epsR * x + x 时,这个性质不再被满足,就有可能出现误差。为了更好地展示这一点,我们可以观察一下 epsR * x + xM(小数位)的最后几位是在何时被抛弃的。

理想情况下,不丢失精度的表示方式:

1
2
3
...00      1
// ↑ ↑
// 最末位 虚构的比最末位还要末的位

实际的表达方式:

1
2
3
...01
// ↑
// 最末位,四舍五入后发生进位

那么,只剩下最后一块拼图,就是为什么舍弃精度后 epsR * x + x 会转换成比 epsR * x + x 大的最小的双精度浮点数,而不是比 epsR * x + x 小的最大的双精度浮点数?(注:这里讨论的都是 x 不是 2 的正整数次幂的情况)

其实答案非常简单。不难看出此时阶码位 E 是由 epsR * x 确定的,因为 + x 部分太小了(x[0,2531]),不足以影响到 E。此时,存在 E1023=53+log2x,故此时该双精度浮点数的最小粒度是 21+log2x=2×2log2x 也即「比 x 大的最小的 2 的正整数次幂」。所以 x 必然大于最小粒度的一半。

上述文字论述,可以转化为更直观的图表。在数轴上,存在如下几个数:

epsR乘以x附近的双精度浮点数.jpg

注意,epsR * x + H/2 不是双精度浮点数,图示仅仅是为了方便参考大小关系。

那么自然而然,舍入时就会舍入到「比 epsR * x + x 大的最小的双精度浮点数」上了,而不会舍入到更远的「比 epsR * x + x 小的最大的双精度浮点数」。

所以我们可以从这个例子出发,进行一些普适的论证:刚刚说到的这个新的数字就会等于 epsR * x 加上一个比 x 大的最小的能够满足 M 的最小粒度的数字了(我们就叫它数字 B 吧);随着 253×A 中的 A 不断变大,能满足 M 的最小粒度的 B 也会不断变大。又因为最终目标是 minB=253×A×M,自然要想取到 B 的最小值,等价于取到 M 的最小值,也就是 252

所以,得到结论:

B=253×A×52=2A+1=2Math.ceil(Math.log2(x))
那么我们要想把 B 从「比 epsR * x + x 大的最接近它的双精度浮点数」里分离出来,自然就再减去一个 epsR * x 就可以了。因为「比 epsR * x + x 大的最接近它的双精度浮点数」和 epsR * x 都可以被双精度浮点数精确表出,所以这一步不会损失精度。

最后就得到了 B,也就是比 x 大的最小的 2 的正整数次幂。

补充说明:当 x 本身就是 2 的正整数次幂时,上述所有计算都不会有精度损失,所以最后返回的还是 x。因此本文标题为了严谨,给出的说法是「比给定整数 x 大(或相等)的最小的 2 的正整数次幂」。

来源:https://blog.jiejiss.com/