随机数的精度问题
上一篇文章谈到蒙特卡洛方法计算圆周率，其中产生高质量随机数的代码：
var arr = crypto.getRandomValues(new Uint32Array(1)); var x = arr[0] / 0xFFFFFFFF; x.toString(2); // 32位循环一次，因此熵仍然是32bit
按照随机数的规范，随机数的区间应该是[0, 1)，因此有必要校正为：
var arr = crypto.getRandomValues(new Uint32Array(1)); var x = arr[0] / (0xFFFFFFFF + 1); x.toString(2);
然而，上述方法产生的[0,1)之间的随机浮点数实质上只有32比特的熵，并不能视为一个足够精度的浮点数，也即这些小数不可能分布在一些间隙中。
这是因为将[0,1]区间等分为0xFFFFFFFF + 1份，而随机数只是落在等分点上。
那么如何让这些产生的随机数的间隙更加细腻呢？
1.双精度浮点数规范
根据双精度浮点数的规范，双精度浮点数使用64比特来储存，第1比特是符号位，第2至第12比特是2的指数(exponent)，在实际运算时需要减掉0x3ff，第13至第64比特是尾数(mantissa)。
除了0和无穷大这几个特例之外，浮点数的实际值是： (符号位) 2^(exponent – 0x3ff) * (1.mantissa)。
尾数(mantissa)部分共有52比特，如果用bit0~bit51表示其52比特，那么在转换为小数时的公式是 bit0 * 2^(-1) + bit1 * 2^(-2) + bit2 * 2^(-3) + … + bit51 * 2^(-52)。
下文的讨论不考虑符号位。下文中，exponent指存储的数值，而指数部分为实际值，也即exponent – 0x3ff。
浮点数如果要大于等于1，exponent显然要大于等于0x3ff。exponent的最大值是：(111 1111 1111)2进制，也即0x7FF，0x7FF在双精度浮点数中特别定义为无限大。
如果忽略尾数部分，只靠指数部分，那么浮点数能够表达的最大整数是2^(0x7FE – 0x3ff) = 2^1023。如果在浏览器的控制台中，输入Math.pow(2, 1023)和Math.pow(2, 1024)会明显看到这一分界点。
Math.pow(2, 1023); // 8.98846567431158e+307 Math.pow(2, 1024); // Infinity
虽然2^1023是如此之巨大，但实际上由于精度限制，2^1022到2^1023之间的整数无法全部精确表达。
实际上，2^1022下一个能够被表达的整数是2^1022 * (1 + 2^(-52))，也即2^1022 + 2^970。
2.连续整数的最大值
这里有一个很有趣的问题：双精度浮点数能够精确表达的连续整数的最大值是多少？
在研究这个问题之前，首先看整数100是如何在浮点数中编码的。
我们知道，100 > 64 = 2^6 且 100 < 128 = 2^7，那么exponent是：6 + 0x3ff = 0x6 + 0x3ff = 0x405 。
尾数部分 (100 – 2^6) / 2^6 = 36 / 2^6 = 32 / 2^6 + 4 / 2^6 = 2^(-1) + 2^(-4) = 0.5625。
用52bit来编码的话，就是：0x9 0000 0000 0000。
和符号位以及exponent合并起来，整个64bit是：0x4059 0000 0000 0000
这里有个在线的浮点数比特位查看器：http://www.binaryconvert.com/convert_double.html
因此虽然100落在64到128之间，但由于尾数部分的调整，所以仍然能够精确表达。
那么100的下一个浮点数是多少呢？
我们知道100是2^6 * (1 + 2^(-1) + 2^(-4))，那么下一个浮点数就是2^6 * (1 + 2^(-1) + 2^(-4) + 2^(-52))，它们之间的差值是2^6 * 2^(-52) = 2^(-46)。
所以随着指数部分的增加，尾数部分能够微调的最小步进越来越小。也就是说精度随着数的增大而下降。
所以当尾数部分能够精确微调的最小步进大于1时，就是双精度浮点数能够精确表达的连续整数的终点。
当指数部分为52时，尾数部分的最小步进为1 = 2^52 * 2^(-52)。
指数部分为52时的最大整数为：2^52 * (1 + 2^(-1) + 2^(-2) + 2^(-3) + … + 2^(-52))，下一个整数是2^53。这就是连续整数的终点。
所以在浏览器的控制台中，看到这种现象就不神奇了（在指数部分为53时，尾数能够调整的最小步进是2）
Math.pow(2, 53); // 9007199254740992 var i = 9007199254740993; i; // 9007199254740992 i = 9007199254740994; i; // 9007199254740994 i = 9007199254740995; i; // 9007199254740994 i = 9007199254740996; i; // 9007199254740996
3.随机数小数部分的精度
回到随机数问题。我们产生的随机数的目标区间是[0, 1)。然而在[0, 1)之间，双精度浮点数并没有那么简单。
由于双精度浮点数的小数部分的精度表达跟指数部分密切相关。当指数小于1时，指数部分对精度有加成作用。
其实很简单，在[0.5, 1)区间，能够精确表达的浮点数有2^52个，它们在区间内是等距的，而在[0.25, 0.5)区间，浮点数也有2^52个，它们在区间内也是等距的。对于[0.125, 0.25)等等等区间也是如此。因此随着浮点数变小，步进也越来越小。
如果我们使用52随机比特产生[1, 2)区间的随机数，很显然指数部分恒为0，尾数部分的步进是2^(-52)，双精度浮点数完全能够精确表达这范围的随机数。
如果使用52随机比特， 记为bit0~bit52，按公式：bit0 * 2^(-1) + bit1 * 2^(-2) + bit2 * 2^(-3) + … + bit51 * 2^(-52)得到的随机数，这些随机数最小值0， 最大值2^(-1) + 2^(-2) + 2^(-3) + … + 2^(-52)，也即二进制0.1111111111111111111111111111111111111111111111111111 (=0.9999999999999998)，符合在[0, 1)区间均匀分布。那么这些纯小数是否能完全被双精度浮点数表达？这关系到随机数的均匀性是否会因精度的损失而被变形。
答案是，[0, 1)区间以 2^(-52)为步进的纯小数能够完全被双精度浮点数表达。
比如非0的最小的小数，(二进制)0.0000000000000000000000000000000000000000000000000001，也即2^(-52)，它的双精度浮点数很容易就计算出为0x3CB0000000000000。(0x3CB – 0x3FF = -52)。这远远没有达到双精度浮点数的能力下限。双精度浮点数的exponent，最小值可以取到0x001，指数为0x001 – 0x3FF = -1022，也即最小的浮点数为2^(-1022)。
最大的小数，(二进制)0.1111111111111111111111111111111111111111111111111111，bit0为1，也即2^(-1)的系数为1，这系数可以放进指数部分里。那么 2^(-1) + bit1 * 2^(-2) + bit2 * 2^(-3) + … + bit51 * 2^(-52) = 2^(-1) * ( 1 + bit1 * 2^(-1) + bit2 * 2^(-2) + … + bit51 * 2^(-51))，稍微复习一下双精度浮点数的规范，很显然指数部分为-1， 尾数部分为1111111111111111111111111111111111111111111111111111 (51位)，补齐到52位，也即11111111111111111111111111111111111111111111111111110。 写成hex，那就是0x3FEFFFFFFFFFFFFE。这个过程叫约化。
因此无论是从哪一bit开始为非0值，它都能约化到指数部分，而指数部分的下限值为-52，不可能超过双精度浮点数的-1022下限。而剩余的比特的位数只能少于52位，约化后实际是左移位后作为尾数部分。尾数部分不足52位，在后面补齐0。
4. 52bit精度的蒙特卡洛法计算圆周率
经过上述讨论，按照浮点数的规范，随机数理论上间隙可以小至2^(-52)，精度将比2^(-32)得到极大提高。
为了产生52bit的精度，我们需要利用两个Uint32，第一个整数左移20位，第二个整数右移12位（舍弃12位），然后两数相加，得到52比特的随机整数（双精度浮点数形式）。因为第一个整数的左移20位运算已经超出32位整数范围，不能够再使用移位运算符。简单地，乘以2^20即可。52比特的随机整数，除以2^52，等到[0, 1)区间的纯小数，即为我们所需。
function getRandomFloat() {
var arr = crypto.getRandomValues(new Uint32Array(2));
var mantissa = arr[0] * 1048576 + (arr[1] >>> 12);
return mantissa * Math.pow(2,-52);
}
var m = 0, n = 10000;
for (var i = 0; i < n; ++i) {
var x = getRandomFloat();
var y = getRandomFloat();
if (x * x + y * y <= 1) {
m++;
}
}
var pi = 4 * m / n;
console.log(pi);
评论