浮点数

浮点运算的使用者可以被描绘为一个智商（IQ）正态分布（或者说钟形曲线）图的模样，因为他们对于浮点运算的理解和使用进程往往会经历三个阶段，形成一种典型的分布态势：

初级程序员通常会无处不在的使用浮点数，好像它是一种神奇的无限精度的数据类型。
然后他们发现0.1 + 0.2并不等于0.3，或者遇到其他的奇怪问题，他们会感到困惑和恐慌，开始误以为每次运算都会加入一些随机的误差，因此在很多年里，他们避免使用任何实数（浮点数）类型。
最后，他们终于鼓起勇气去读IEEE-754浮点数标准，理解浮点数是如何工作的，并开始适应如何正确地使用浮点数。

不幸的是，仍有许多人停留在对于浮点数运算的误解阶段。他们潜在的认知误区导致他们对浮点数运算存有种种错误认知——以为它基本上是不精确和不稳定的，而且比整数运算慢。

但这些都是错误的认知。首先，浮点运算其实往往比整数运算更快，这是因为有专门的指令用于加速浮点数的运算。其次，实数（如浮点数）的表示方式是有严谨标准的，并且在进行运算时，也遵循简单且确定的规则，比如四舍五入等舍入规则。这使得我们可以可靠地处理计算误差。

实际上，浮点数的使用是非常可靠的。甚至一些高级编程语言，特别是JavaScript，就根本没有整数类型。在JavaScript中，只存在一种数值类型，这种类型在内部以64位双精度浮点数（64-bit double）的形式储存。由于浮点数运算的工作方式，所有在−2^53到2^53之间的整数和涉及它们的计算结果，都可以精确地储存，所以从程序员的角度来看，几乎没有要区分单独的整数类型的需求。

一个值得注意的例外是，当需要对数字进行位操作时，通常浮点数单元（floating-point units，负责浮点数运算的协处理器）是不支持的。因此，这就需要将浮点数转换为整数。这个过程在JavaScript启动的浏览器中非常常见，以至于ARM增加了一条名为“FJCVTZS”的特殊指令，其名称可以解读为“将JavaScript中的浮点数向零舍入，转化为有符号整数”。这个指令的功能就是将浮点数以和JavaScript完全相同的方式转化为整数。这是软硬件之间相互反馈作用的一个有趣例子：软件应用（这里指大量的JavaScript代码对数值进行位操作的需求）驱动硬件设计（这里指ARM处理器新增的FJCVTZS指令）来更好地满足运算需求。

除非你是一个只使用实数类型来模拟整数算术的JavaScript开发者，否则你可能需要一个关于浮点数算术更深入的指导，所以我们将从一个更宽泛的主题开始。

实数表示

如果你需要处理实数（非整数），你有几种不同的选择，适用性各不相同。在直接跳到浮点数（这是本文的主要内容）之前，我们想讨论一下现有的其他选择以及它们背后的动机 — 毕竟，避免使用浮点算术的人确实有他们的理由。

符号表达式

首先，最麻烦的方法储存的不是结果值本身，而是存储产生这些值的代数表达式。

这里有一个简单的例子。在某些应用中，如计算几何学，除了加、减和乘三种运算之外，你还需要进行无需取整的除法运算，产生一个有理数。这个有理数可以用两个整数的比率（分数）来精确表示出来：

struct r {
    int x, y;
};

r operator+(r a, r b) { return {a.x * b.y + a.y * b.x, a.y * b.y}; }
r operator*(r a, r b) { return {a.x * b.x, a.y * b.y}; }
// r operator/(r a, r b) { return {a.x * b.x, a.y * b.y}; } // 原作者的实现
r operator/(r a, r b) { return {a.x * b.y, a.y * b.x}; }
bool operator<(r a, r b) { return a.x * b.y < b.x * a.y; }
// ...and so on, you get the idea

这个比率可以被化简到不可约分的形式，这样它就会变为一种独一无二的表示形式。

struct r {
    int x, y;
    r(int x, int y) : x(x), y(y) {
        if (y < 0)
            x = -x, y = -y;
        int g = gcd(x, y);
        x /= g;
        y /= g;
    }
};

这就是诸如WolframAlpha和SageMath这样的计算机代数系统的工作方式：它们仅操作符号表达式，避免将任何东西计算为实数值。

用这种方法，你可以得到绝对精度，它在有限的范围内运作良好，如仅支持有理数。但是这需要很大的计算成本，因为通常情况下，你需要以某种方式存储产生当前结果的所有历史操作，并在每次进行新的操作时考虑这些历史操作。当历史操作过多时，这会变得难以管理。

定点数

另一种方法是坚持使用整数，但将它们视为乘以一个固定常数。这个方法基本上相当于将使用的单位换成了一个更大的单位，这样可以将实数运算问题转变为整数运算问题。

因为有些值无法精确表示，所以这种方法使计算得到的结果不精确：你需要将结果四舍五入到最近的可表示的值。

这种方法常用于财务软件中，那里真的需要一种直接的方式来管理舍入错误，以便最终的数字相加。例如，纳斯达克在其股票列表中使用的基本单位是1/10000美元，这样可以确保在所有交易中刚好可以有4位小数的精度。

struct money {
    uint v; // 1/10000th of a dollar
};

std::string to_string(money) {
    return std::format("${0}.{1:04d}", v / 10000, v % 10000);
}

money operator*(money x, money y) { return {x.v * y.v / 10000}; }

除了引入舍入误差，更大的问题是，当比例常数使用不当时，它们会变得毫无用处用。如果你处理的数字太大，那么表示它的整数值将会溢出；如果数字太小，它们将被四舍五入为零。有趣的是，前一种情况曾经在纳斯达克遇到过，当伯克希尔·哈撒韦的股票价格接近(2^32−1)/10000 = $429,496.7295时，它就无法再适用于一个无符号的32位整数来表示。

这个问题使得定点运算在需要同时处理小数和大数的应用中基本上都是不合适的，例如，在评估某些物理方程时：

$E = mc^2$

在这个特殊的例子中， $m$ 通常与一个质子的质量（ $1.67⋅10^{−27}$ kg）的数量级相同，而 $c$ 是光速（ $3⋅10^9$ m/s）。

浮点数

在大多数数值应用中，我们主要关注的是相对误差。我们希望计算结果与真实值的误差不超过0.01%，我们并不真正关心那个0.01%在绝对单位中等于多少。

浮点数通过存储最重要的一部分数字和数值的数量级来解决这个问题。更确切地说，它们用一个整数（称为有效数（significand）或尾数（mantissa））表示，并使用某个固定基数（base）的指数（exponent）进行缩放 - 最常见的是2或10。例如：

$1.2345 = \underbrace{12345}_{\text{mantissa}} × {\underbrace{10}_{\text{base}}}^{\overbrace{-4}^{\text{exponent}}}$

计算机工作在固定长度的二进制字上。因此，当设计一个用于硬件的浮点格式时，你会希望有一个固定长度的二进制格式，其中一些位被分配给尾数（用于更大的精度），一些位分配给指数（用于更大的范围）。

下面这个自己实现的浮点数操作可以帮助我们理解这个概念：

// DIY floating-point number
struct fp {
    int m; // mantissa
    int e; // exponent
};

这种方式可以表示形式为 $±m×2^e$ 的数字，其中 $m$ 和 $e$ 都是有界的，并且可能是负整数。这对应于负数或小数。这些数字的分布非常不均匀：在 $[0,1)$ 范围内的数字数量大致与 $[1,+∞)$ 范围内的数字数量一样。

请注意，一些数字的这种表示形式并不唯一。例如，数字1可以被表示为以下形式：

$1 × 2^0 = 2 × 2^{-1} = 256 × 2^{-8}$

以及其他28种不会使尾数溢出的方式。

这种多重表示可能对一些应用（如比较或哈希）构成问题。为了解决这个问题，我们可以使用一定的规则来规范化这些表示方式。在十进制中，标准形式是始终把小数点（逗号）放在第一个数字后面（例如6.022e23），对于二进制，我们可以做同样的规定：

$42 = 10101_2 = 1.0101_2 × 2^5$

遵循这个规则，则首位总是 1。存储它是冗余的，所以我们只是假装它在那里，而只存储其他位，这些位对应于 $[0,1)$ 范围内的某个有理数。现在，可表示的数字集大约是这样的：

$\{ ±(1+m)\cdot2^e | m = \frac{x}{2^{32}}, x \in [0, 2^{32}]\}$

由于现在 $m$ （尾数）是一个非负数，我们将把它变成一个无符号整数，并添加一个独立的布尔字段来表示数字的正负号。

struct fp {
    bool s;     // sign: "0" for "+", "1" for "-" 
    unsigned m; // mantissa
    int e;      // exponent
};

让我们尝试使用我们自己创建的浮点数实现一些算术操作，例如乘法。使用新的公式，结果应该是如下的形式：

$c = a \cdot b \\ = (s_a \cdot (1 + m_a) \cdot 2^{e_a} \cdot (s_b \cdot (1 + m_b) \cdot 2^{e_b}) \\ = s_a \cdot s_b \cdot (1 + m_a) \cdot (1 + m_b) \cdot 2^{e_a} \cdot 2^{e_b} \\ = \underbrace{s_a \cdot s_b}_{s_c} \cdot (1 + \underbrace{m_a + m_b + m_a \cdot m_b}_{m_c}) \cdot 2^{\overbrace{e_a + e_b}^{e_c}}$

通过分组计算。但是，还有两个需要注意的细节：

新的尾数（也叫有效数字或者小数部分）现在在 $[0,3)$ 的范围内。我们需要检查它是否大于1并且对表示进行规划化处理，这会用到下面的公式： $1+m=(1+1)+(m−1)=(1+\frac{m−1}{2})⋅2$
由于精度的限制，结果可能（很可能）无法精确表示。例如在乘法操作中，可以看到在计算两个大数字 $ma \cdot mb$ 时，其结果需要的位数是单个数需要的位数的两倍。为了应对这个问题，我们需要选择一个能够用我们的有效位数来表示的最接近的数。在具体做法上，常见的做法是采用四舍五入的方法进行处理。

为了正确处理尾数溢出问题，我们从尾数 $m$ 中预留一位来处理溢出。这就是为什么把 $m$ 限制在 $[0,2^{31})$ 范围内的原因。

fp operator*(fp a, fp b) {
    fp c;
    c.s = a.s ^ b.s;
    c.e = a.e + b.e;
    
    uint64_t x = a.m, y = b.m; // casting to wider types
    uint64_t m = (x << 31) + (y << 31) + x * y; // 62- or 63-bit intermediate result
    if (m & (1<<62)) { // checking for overflow
        m -= (1<<62); // m -= 1;
        m >>= 1;
        c.e++;
    }
    m += (1<<30); // "rounding" by adding 0.5 to a value that will be floored next
    c.m = m >> 31;
    
    return c;
}