浮点数和大数(转自王晓华)

--write by zhuwx 2019-06-24 20:58:16 +0800 CST

点击量:17

整数的范围

“玩”算法,对各种类型的“数”能表达的数据范围要心里有数。

  • 8 位有符号整数能表示的范围是 −128 到 127,无符号时能表示的范围是 0 ~ 255。
  • 16 位的有符号整数能表示的范围是 −32768 到 32767,无符号时能表示的范围是 0 ~ 65535。
  • 32 位的有符号整数能表示的范围是 −2,147,483,648 到 2,147,483,647,无符号时能表示的范围是 0 ~ 4,294,967,295。
  • 64 位的有符号整数能表示的范围是 −9,223,372,036,854,775,808 到 9,223,372,036,854,775,807,无符号时能表示的范围是 0 ~ 18,446,744,073,709,551,615。

什么叫心里有数?32 位整数最大能表示大概 42 亿多一点,如果算法要表示某地的人口数,设计数据模型时用 32 位整数也就够了,但是要表示全球人数就不行了。一光年大约 9,460,730,472,580,800 米,从地球到仙女座星系的距离是 254 万光年,已经超出了 64 位整数能表达的范围。我们在设计数据结构的时候,要对问题的规模做到心中有“数”,不要用错数据类型。

大整数库

超过 64 位的超级大数怎么办?如果你要解决的问题需要这么大范围的数,可以考虑使用大整数。Python、Lisp 等语言都内建了大数计算机制,Java 也有 BigInteger 和 BigDecimal 分别用于大整数和大实数的计算。然而,C/C++ 没有内建的机制或类型支持大数计算,但是 C/C++ 有很多高性能的第三方大数库可供选择。

GMP 开源大数库是一个任意精度的大整数运算库,它包括了任意精度的整数、浮点数的各种基本运算操作。它是一个 C 语言的库,但是官方提供了 C++ 的包装类,根据其官方的资料宣称,GMP 是目前地球上速度最快的大数库。

Miracl 库是 Shamus Software Ltd 开发的一个大数库,它的使用许可针对教育科学研究或者非商业目的应用是免费的。它是一个 C 语言的库,同时提供了几个较为简单的 C++ 包装类。在功能上它不但提供了高精度的大整数和分数的各种数学运算操作而且还提供了很多密码学算法中的功能模块,如 SHA、AES、DSA 等中的一些底层操作,据说局部算法使用了汇编代码实现,所以速度应该还是不错的。

除了专用的大数库,很多密码学软件中也会包含大数计算相关的内容,比如开源的密码学计算库 Crypto++ 中就包含了整套的大数计算库,OpenSSL 中也有对大整数计算的支持。

自己设计大整数

OJ 上我遇到过这类题目,还好只是要实现一个大数加法。设计大整数常用两种数据存储方式,一种是字符串方式,如 46744073709551615,每一个字符就是十进制数字的一个位;字符串方式的优点是简单、直观,缺点当然是计算时的效率不高,特别是乘法和除法的计算,另外,字符串方式的另一个缺点是占用存储空间比较大,四个字节的字符串最大只能表示到 9999。设计大整数的另一种存储方式是采用 “2n2n 进制数”的形式表示大数,对于不支持 64 位整数的系统,可采用 216216 进制,用一个 16 位整数表示大数的一个位,这样每个位的有效数字是 0 ~ 0xFFFF,可以对比十进制整数一个位的有效数字 0 ~ 9 来理解,超过了 0xFFFF 就需要进位。在数据结构存储时,可以采用 16 位无符号整数的数组来存储大数的每个位,假如某个用 216216 进制表示的大数数组末两个元素是:

0xFFFF     0xFFFF

可以理解为是 216216 进制数的“99”,如果给它再加 1,就会产生进位,得到 216216 进大数“100”:

              FFFF        FFFF
+                            1
--------------------------------
=       1        0           0

为什么选用 216216 进制,而不是直接用 232232 进制呢?是因为计算会溢出,用 32 位整数临时存储两个 216216进制数位的计算不会溢出(乘法最终也是转化成加法实现的)。如果系统支持 64 位整数,就可以直接用 232232 进制表示大数。

大整数的加法是最简单的算法,下面我们用 232232 进制大数为例,演示一下大数加法算法实现(只考虑同号相加,异号相加可以调整为减法实现),我们的大数数据模型定义为:

#define MAX_BI_LEN   256

typedef struct 
{
    unsigned int length;  //大数在2^32进制下的数字位数    
    unsigned int value[MAX_BI_LEN];  //大数每一位的值
}BigInt;

这个数据结构定义没啥好说的,看注释就好了。加法算法实现首先将其中一个加数赋值给 result,并将 result 的位数先调整为两个加数中位数最长的那个值,然后开始按位加,并用 carry 记录进位。两个 32 位整数相加,存入 64 位的中间数 sum,sum 的低 32 位是两个加数位的和,高 32 位中如果有溢出就是进位。carry 初始时是 0,保证第一次加法计算时不会多加。最后循环完成各位相加后,如果 carry 不是 0,说明当前最高为相加又产生进位了,需要将结果也调整增加一位。

void BigIntAdd(const BigInt& value1, const BigInt& value2, BigInt& result)
{
    BigInt result = value1;

    unsigned int carry = 0;
    /*先调整位数对齐*/
    if(result.length < value2.length)
        result.length = value2.length;
    for(unsigned int i = 0; i < result.length; i++)
    {
        unsigned __int64 sum = value2.value[i];
        sum = sum + result.value[i] + carry;
        result.value[i] = (unsigned long)sum;
        carry = (unsigned)(sum >> 32);
    }
    //处理最高位,如果当前最高位进位carry!=0,则需要增加大数的位数
    if(carry != 0)
    {
        result.value[result.length] = carry;
        result.length += 1;        
    }
}

是不是很简单? play_with_algo 中有一个 CBigInt 类,用的就是 232232 进制大数的形式,支持简单的加、减、乘、除运算,加、减、乘、除运算都是参考《编译原理》书中介绍的加法器、减法器、乘法器和除法器的原理实现的,效率不是很高,但是还可以一用。用于高性能计算就算了,通过其实现代码了解一下大数计算的原理还是可以的。

浮点数的精度

IEEE 定义了多种浮点数的二进制表示形式,以 IEEE 754 标准为例,单精度浮点数(float)用 32 位二进制位表示,其中 31 位是符号位,23 ~ 30 共 8 个比特位表示指数,0 ~ 22 共 23 个比特位表示有效数字。双精度浮点数(double)用 64 位二进制表示,其中最高为 63 位是符号位,52 ~ 62 共 11 个比特位表示指数位,0~ 51 共 52 个比特位表示有效数字。IEEE 还定义了用 80 个比特位表示的扩展双精度类型,不过大家用的机会不多,这里就不介绍了。

根据双精度浮点数的定义,其能表示的数的正数范围大概是 1.7E−308 〜 1.7E + 308(负数范围类似),既然这么大,那表示 254 万光年是否绰绰有余?当然,从数的范围来讲,确实是绰绰有余,但是考虑到精度,用双精度浮点数承载天文计算就还差点儿火候。双精度浮点数有 52 个比特位表示有效数字,换算成十进制的话也就是大约 16 个十进制位的有效数字,这个精度用作科学计算显然有点尴尬,科学计算,还得靠大数库。

为什么很多数都强调用 double,尽量不用 float?也是这个原因,float 类型的浮点数有效数字只有 23 个比特,换算成十进制差不多 7~8 个有效数位,精度实在是太低了。

浮点数的误差

浮点数的误差其实也是由精度引起的,以双精度浮点数为例,按照十进制方式其有效数字大约 16 位,换句话说,超过 16 位的小数用双精度数是区分不出来的。所以你要知道,这两个数计算机是区分不出来的,在计算机来看,它们两个是相等的:

    double aaaa = 0.1234567890123456789;
    double bbbb = 0.123456789012345678;

除了超出精度范围的数之外,在精度范围内的数也会存在误差,这个误差的主要原因是用二进制小数表示十进制小数时存在转换误差,看看下面的代码,你以为它们在内存中的值是什么?

    double first = 0.3;
    double second = 0.4;

    double sum = first + second;

来看看:

enter image description here

是不是和想象中的不一样?我们来看看 0.3 怎么转化为二进制,转化方法就是不停地乘 2,有进位到 1 了就去掉 1,并将对应的位设置为 1,如果没有进位,就将对应的位设置为 0。来看看这个过程:

0.3 X 2 = 0.6                              0

0.6 X 2 = 1.2                              1

0.2 X 2 = 0.4                              0

0.4 X 2 = 0.8                              0

0.8 X 2 = 1.6                              1        //0.6又出现了

可以看到 1001 会重复出现,无限循环,最后的结果是 0.0100110011001…,直到超出浮点数的有效位数。这就是我说的二进制小数不能精确表达十进制小数的原因。并不是所有的十进制小数都不能用二进制小数精确表达,那么什么样的值能精确表达呢,就是那些能用 1/2、1/4、1/8、1/16……准确累加的值,比如 0.5、0.75(1/2 + 1/4)或 0.625(1/2 + 1/8)这类的小数。

再来看看下面的代码,你以为结果如何:

    double r1 = 16.1 * 100 + 0.9 * 100;
    double r2 = 17.0 * 100;
    if (r1 == r2)
    {
        std::cout << '16.1 * 100 + 0.9 * 100 == 17.0 * 100' << std::endl;
    }

实际上我们期望的那个输出打印并没有出现,因为 r1 不等于 r2 !真实情况是 r1 = 1700.0000000000002、r2 = 1700.0000000000000。还是同样的道理,这就是我在课程中一直强调不要用浮点数判断相等的原因,大多数情况浮点数比较都没有太大问题,但是总有你想不到的那个情况出现,为了不给自己找麻烦,最好不要直接用浮点数直接判等,当然,直接判断不等于也尽量避免。那怎么判断两个浮点数相等呢?基础部分内容中有答案。

总结

只要涉及到用计算机表示数,都不要想当然,了解一点整数和浮点数与二进制存储方式之间的原理,有助于减少困惑和对着代码排除错误的时间,最后:

“世界上有 10 种人,一种是懂二进制的,一种是不懂二进制的。”