计算机中的浮点数

在计算机中对于数的表达，都是使用二进制的方式。对于浮点数也是，不同于整数的表示，浮点数的表示复杂不少。目前而言，大多数的电脑都是使用IEEE 754标准所定义的表示方法。本文的目的是，讲解IEEE 754的基本格式，并且来解释一些浮点数的问题。

定点法来表示二进制数

先不管定点法的意思是什么。先来回想一下十进制表示小数的方法，这个非常简单，比如说123.456。类似地，可以得到一个通用化的表达。对于任意的十进制数:

$$
d_{m}d_{m-1}d_{m-2}\dots d_{1}d_{0}.d_{-1}d_{-2}\dots d_{-n}
$$

可以表示为：

$$
d = \sum_{i=-n} ^{m} 10^ {i} * d_{i}, 其中d_i \in [0,9]
$$

这看着有点抽象，比如说十进制数12.34，就可以表示为: $1\times10^{1}+2\times10^{0}+3\times10^{-1}+4\times10^{-2}=12.34$. 我们对一个数采用类似的办法来表示一个二进制浮点数。如下：

$$
b = \sum_{i=-n} ^{m} 2^ {i} \times b_{i},其中b_i=0 或者 1
$$

比如说，二进制数101.11，就可以表示为：$1\times 2^{2} + 0\times2^{1}+1\times2^{0}+1\times2^{-1}+1\times2^{-2}$

对于这种方法来表示一个二进制数，我们可以观察到，小数部分都是以2的倍数作为分母的，比如说0.5就是二分之一，0.25就是四分之一。我们可以明确的表达这种2的倍数作为分母的数字，对于如0.3，0.1这种数，因此只能近似的来表达。不过我暂时还不知道何种算法来求一个浮点数的近似表达的二进制串。

此外，这种方法还具有一个显著的缺点，无法表示特别大的数，适用范围有限。比如说，如果想要表示$5*2^{100}$这样的数，在101(二进制形式的5)后面还需要跟着一大堆的零。显然是不合适的。

上面这种形式来表示二进制数叫做定点法，因为将小数点固定在了某一个位。我们希望能够引入一种方法，可以表示更大范围的数字。这种方法就是我们熟悉的科学计数法。

IEEE浮点表示法

科学计数法我们很熟悉，如123.456就可以表示为$1.23456*10^{2}$。科学计数法可以表示更大的范围。首先先给出IEEE 754对于浮点数的基本表示，下图是32位浮点数的基本形式（float）：

下图是64位浮点数的形式(double)：

IEEE标准用$V=(-1) ^ {s} \times M ^{E}$这样的形式来表示一个数，64和32位浮点数的基本结构一样，下文以32位浮点数来讲解。

符号(sign): 第一个bit来表示符号，0为正数，1为负数，上面公式中的s。
尾数（fraction）：从0到22 bit作为尾数，总共23 bit，上面公式中的M。
阶码（exponent）: 从23到30 bit作为解码，总共8 bit，上面公式中的E。

这确实非常的不够直观，没有具体的例子来讲述，为了内容的完整性，将概念讲完后再来举例子阐述。

IEEE 标准中，fraction部分有不同的情况，取决于阶码（exponent）是否0。对于浮点数，IEEE 标准定义了如下几种格式：

Normalized 大多数情况下的单精度浮点数都是这种形式的。这种形式之下，阶码（exponent）被解释为$E = e- Bias$.下面我们对这些符号作解释。
- e：表示8 bit可以表示的整数，在Normalized，e 的取值范围就是（0，255），PS：上图中输了Normalized 情况下阶码字段不能为255或者1，所以取值范围是（0，255）。注意，e并不是阶码中的值，E才是。
- Bias：公式$Bias=2^{k-1}-1$，k是阶码的位数，所以得到在单精度浮点数中，Bias = 127，相类似的，在双精度中Bias = 1023。
综上所述，我们来计算阶码的取值范围，带入到$E = e- Bias$,得到了阶码的取值范围为[-126,127]。
对于fraction这部分，M = 1+f (上图中灰色的部分)。f这部分就是前面已经说到过的，小数的二进制表示，如0.5=100...0，0.25= 010...00。没理解没关系，待会再来说具体的例子。
De normalized当阶码都为0的时候，，此时$E = 1- Bias$，注意和上面的区别。于是，De normalized情况中，阶码总是为126，通常用来表示非常小的值。因为$2^{-126}$真的是一个非常小的数。对于fraction这部分，M = f,注意和前面的区别。
当f部分全部都是0的时候，此时就是表示浮点数的0了。因为符号位的关系，这里会存在 +0和 -0，但是问题不大。
特殊情况在阶码都是全为1的时候，视为特殊值。如果此时的f部分为0，那么就是就无穷了。符号号组合起来就可以分别表示正无穷和负无穷。如果f不为0，那么就是表示NaN，NaN的意思是 = Not a Number。

下面是一个单精度浮点数的例子

一个32 bit的浮点数如下：

首先来看fraction部分，001...000对应的小数是：0.125，再看阶码部分。1000 0000 = 128，然后减去 bias。得到阶码为 128-127 = 1。符号位位0，表示的是正数。可得，这个数为$ 2 ^{1} * (1+0.125) = 2.25$。别忘了在normalized情况下，M = 1 +f 。

浮点数的加法

讲解浮点数的加法是为了接下来的讲解IEEE 754表示法在运算的时候所带来的的问题做准备的。这一块并不复杂。

回想十进制下科学计数法的加法。科学计数法加法的规则为：

首先将两个数都化为相同阶码
尾数与尾数相加，如果超过了10，那么阶码加一。

比如说$1.2 \times 10^{2} + 9.3\times10^{3} = 0.12 \times 10^{3} + 9.3 \times 10 ^{3} = 9.42 \times 10^{3}$ 。

PS：这个例子说明的是，我们在移动尾数的时候，都是移动小的那个。上面的1.2小数点左移以为成为了0.12，因为我们需要保证小数的形式为x.yyyyy这种形式，如果让大的往右移，就会出现xx.yyyy。破坏了格式。况且，我们在IEEE 754中的尾数格式就是需要为x.yyyy形式的。

PS：

这里有一个问题就是，对于normalized的浮点数，默认尾数为 1 +f。而在浮点加法的时候，对尾数右移（也就是小数点左移）。此时，在加法的时候如果我们还对该数解释为Normalized,那么尾数就还是被解释为 1 + f。比如说1.25，小数点左移之后变成了0.125，如果我们还是以1+f来解释，那么就成为1.125了，显然就不行了，我虽然没有找到正式的文档来说明这一点，但是就这里的而言，看起来在此时，1会被忽略。

IEEE 754的一些问题

这里，我所知道的主要就是两个经典的问题，一个是在大多数的编程语言中，0.1 + 0.2 != 0.3。前面说过了，二进制表示小数的决定了计算机只能准确的表达以2的倍数作为分母的小数。比如说0.25.0.5,1.25等等，对于其他的小数只能近似的表达。

在这里，说实话内容还不少。可以详细说说，先看下面代码，这里首先一个问题就是，在C语言中，甚至其他语言中都是类似的。如果我们对一个浮点数如果没有加f后缀，汇编编译器默认为这是double类型的数.因此下面代码中的几个小数都是double类型的。

#include "stdio.h"
int main(int agrc, char **argv) {
    printf("%.55lf\n",0.2);
    printf("%.55lf\n",0.1);
    printf("%.55lf\n",0.3);
    printf("%.55lf\n",0.2 + 0.1);
}

输出：上面的代码在Linx(WSL)和Windows下输出有不同。我以Linux作为环境，在WSL中输出如下：

0.2的64位表示比0.2稍大，0.1也是。所以两者的和也是稍微大于0.3。而0.3的64位表示稍微小于0.3。所以很自然的是相加后的和大于0.3。printf("%.55lf\n",0.1);在windows下的输出非常不一样（java的也相类似），所以就不多说了。

但是这里，这个对于单精度的情况又是有些不同，如下代码：

#include "stdio.h"
int main() {

    float a = 0.1f;
    float b = 0.2f;
    //float有效为保留到小数点后23位
    printf("%.23f\n",a+b);
    printf("%.23f\n",0.3f);
}
//输出
//0.30000001192092895507812
//0.30000001192092895507812

很奇怪，0.1f + 0.2f == 0.3f，和前面0.1 + 0.2 != 0.3又不一样了。是因为计算机总是以最近似的值来表示某个小数，在这里，说明 0. 1 + 0.2的和恰好和最接近0.3的近似的数相等。因为double的精度更高，所以近似于0.3的数并不是0.1+0.2的和。所以会不相等。不过不管怎么说，最主要的是知道，在计算机中，这些数并不能被准确的表达。

可以参考Floating point precision and equality in Java，提到了为什么0.1f + 0.2f = 0.3f。

解决办法

其实也很简单，我们只需要约掉这些低位的数字。比如说round(6)只保留6位即可。或者说，判断Abs(0.1 + 0.2) - 0.3 < foo。我们只需要foo在我们接受的误差内就行。Is floating point math broken?

这种方法看起来比较粗糙,对于货币这些十分精确的场景,在Java中可以使用Big Decimal，Python中Decimal。来自这里

这种误差并不是因为二进制来表示小数所独有的，十进制也有无法正确表达数字，比如说1/3 = 0.33333...33。1/3 * 3 != 0.333...3 。同样，因为截断误差,在取1/3的时候也是近似的。

大数吃小数问题

另外一个问题就是大数吃小数问题。前面说过科学计数法下，加法的运算规则，很重要一点是我们要保证两个数的阶码相同。一个经典的场景就是，一个很大的数加上一个很小的数，因为两者之间的阶码差距很大，于是小的数尾数在不断的右移动就直接变为0了。

#include "stdio.h"
int main(int argc,char *argv[],char *envp[]) {
    float a = 1e10f;
    float b = 3.14f;
    printf("%lf", a +b);
}
//输出
//10000000000.000000

因为3.14和1e10之间的阶码差的太大了，以至于3.14的尾数部分在右移的时候超过了23位，直接变成了0。因此在相加的时候，对1e10的尾数没有造成任何影响。更加正式来说，来自：

An important case occurs when the numbers differ widely in magnitude. If the exponents differ by more than 24, the smaller number will be shifted right entirely out of the mantissa field, producing a zero mantissa. The sum will then equal the larger number.

只要两者的倍数超过了$2^{24}$次，就会出现这种情况。这也就是导致了另外一个问题：浮点数不具有结合律.

#include "stdio.h"
int main(int argc,char *argv[],char *envp[]) {
    float a = 1e10f;
    float b = 3.14f;
    printf("%lf", (b+a)-a);
}
//输出
//0.000000

int main(int argc,char *argv[],char *envp[]) {
    float a = 1e10f;
    float b = 3.14f;
    printf("%lf", b + (a-a));
}
//输出
//3.140000

可以使用Kahan's Summation Formula来解决大数吃小数的问题。这里就不讨论了。以后再学吧

一段小程序

下面是一段C语言小程序，用于查看一个浮点数的二进制，以及得到fraction,exponent，s等信息。

#include "stdio.h"
#include "math.h"

typedef unsigned int uint_32;

union ufloat {
    int a;
    float b;
};

uint_32 show_binary(int a) {
    int len = sizeof(a) * 8;
    uint_32 result = 0;
    printf("binary format:");
    while (len > 0) {
        len--;
        uint_32 b = (1 << len) & a;
        result |= b;
        if (b)
            printf("%c",'1');
        else
            printf("%c",'0');
        if (len % 4 == 0)
            printf(" ");
    }
    printf("\n");
    return result;
}

void show_fields(uint_32 num) {
    uint_32 sign_mask = 0b10000000000000000000000000000000;
    uint_32 exp_mask =  0b01111111100000000000000000000000;
    uint_32 frac_mask = 0b00000000011111111111111111111111;

    uint_32 exp = (num & exp_mask) >> 23;
    uint_32 cursor = 1 << 22;
    double frac_result = 0;
    uint_32  frac = frac_mask & num;
    int index = -1;

    while (cursor > 0) {
        uint_32  foo = (cursor & frac) == 0 ? 0 : 1;
        frac_result += pow(2,index) * foo;
        index -= 1;
        cursor >>= 1;
    }

    if (!(num & sign_mask)) {
        printf("positive,");
    } else {
        printf("negative,");
    }

    if ( exp != 0 && exp != 255) {
        //regular format
        frac_result += 1;
        exp = exp -127;
    } else if ( exp == 0) {
        exp = -126;
    } else {
        if ( frac_result == 0) {
            printf("infinite number");
            return;
        } else {
            printf("Nan");
            return;
        }
    }
    printf("exponent:%d,fraction: %.23f\n",exp,frac_result);
}

int main(int argc,char *argv[],char *envp[]) {
    float a = 1e10f;
    float b = 3.14f;
    printf("%lf", b +a - a);
}

结语

本篇大概的介绍了一下计算机是如何来表示二进制数的，主要是IEEE 754标准的内容。现在大多数机器都是用这个。IEEE 754的作者William Kahan也因此获得了图灵奖（1989）。IEEE 754有一些小问题，比如说浮点数不具有结合律，0.1 + 0.2 != 0.3，以及大数吃小数的问题。在对于精度不是很重要的场合中，只要在可接受的误差内，用float也行。精度高就可以用BigDecimal。