浮点数相关问题

伟大的 zjk 曾经说过：

“为什么不写整数计算几何”。

（此处应有一张图）

这确实很对，不会有过大的精度误差，实现起来也能方便很多。

但是总是有需要使用浮点数的时候，发现自己对 cpp 的浮点数以及一些隐式转换并不是很熟悉，遂在学习计算几何之前先做此记录。

一个三岁小孩都知道的事实是，我们现今所使用的计算机是以二进制方式储存数据的。

浮点数显然也需要使用二进制方式存储，目前比较通用的一种标准是 IEEE 754 浮点数：

对于一个浮点数 \(F\)，我们将其表示为 \(F = (-1)^S \times M \times 2^{E}\)。

其中 \(S\) 称作符号位，取值为 \(\{0, 1\}\)，用于表示正负。

\(M\) 表示有效数字，是一个取值为 \([1, 2)\) 的二进制小数，可以用于描述数的精度。

\(E\) 表示指数位，可以用于描述数的大小范围。

这么说很抽象，实际上你可以将其看作二进制的科学记数法，比如我们将 \(-5.0\) 表示为 \((-1.01 \times 2^{2})_{2}\)，那么它用以上的方式描述就可以写作 \((S,M,E) = (1,1.01,2)\)。

不过实际储存的时候还有一些其它的要求，这里我们以 float（假定为 32 位浮点数）为例。

IEEE 754 标准（实际上为 IEEE 754 binary32 格式）规定，对于浮点数的储存，需要满足以下要求：

  1. 在储存 \(M\) 时忽略掉小数点以及第一位（毕竟取值都固定为 \([1, 2)\) 了），占据 \(23\) 个 bit。
  2. \(E\) 的真实值为 \([-127, 128]\)，在存储时加上 \(127\)，换句话说 \(E\) 的实际取值为 \([0, 255]\)，占据 \(8\) 个 bit。

聪明的你一定能发现，这是有问题的，如果 \(F = 0.xxx\) 的形式呢？

所以还有规定：

  1. 如果 \(E\) 的 \(8\) 位全部为 \(0\)，那么 \(M\) 的第一位不为 \(1\)，而应当为 \(0\)。

同样的，还需要对 \(\infty,\text{NaN}\) 有定义。

  1. 如果 \(E\) 的 \(8\) 为全部为 \(1\)，如果 \(M\) 的 \(23\) 个 bit 全部为 \(0\)，那么根据 \(S\) 的不同表示正负无穷；否则表示 \(\text{NaN}\)。

IEEE 754 当中还规定了几种不同的格式(表格修改自 OI-wiki)：

格式 (IEEE 754 -) 位宽 可表示区间 可保证精度位数
binary32 (单精度浮点数) \(32\) \(\pm 3.4 \times 10^{38}\) \(6 \sim 9\)
binary64 (双精度浮点数) \(64\) \(\pm 1.8 \times 10^{308}\) \(15 \sim 17\)
binary64 extend \(\ge 80\) \(\pm 1.2 \times 10^{4932}\) \(\ge 18 \sim 21\)
binary128 \(128\) \(\pm 1.2 \times 10^{4932}\) \(33 \sim 36\)

可保证精度位数是什么等下再来说。

不同编译器对标准的实现一般都不一样，在这里只讨论 g++，使用版本 12.1.0，9.3.0 还没测。

  • float 对应 binary32 格式
  • double 对应 binary64 格式。
  • long double 对应 binary64 扩展格式（和 __float80 表现一致）

而如果我们想要使用 binary128 格式浮点数，对于 gcc 来说有 <quadmath.h> 中的 __float128 这一扩展，不过在此不做展开。

在 C++23 标准中，新增了 <stdfloat> 库，增加了几种定宽浮点数的支持，其中包括对 binary128 的支持（std::float128_t），不过它们和正常的 float, double 等似乎存在一定隔离。

可以强制转换，但隐式转化存在一些我不太清楚的风险，不过现在还用不到哈哈，具体可以参阅 cpp ref。

众所周知浮点数存在误差，但是是为什么呢？

实际上是因为十进制计数系统和二进制计数系统在转化时出现的不兼容导致的。

比如 0.2，它就没法被准确的表达为二进制浮点数。

用辗转相除法，我们会发现没法除到小数点不为 \(0\)。

不信？你可以试试看，也可以试一试这份代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// author : black_trees
// make yourself an idiot.

#include <iomanip>
#include <iostream>
#include <quadmath.h>

using namespace std;

int main() {
    float val1 = 0.2;
    cout << fixed << setprecision(128) << val1 << endl;

    double val2 = 0.2;
    cout << fixed << setprecision(128) << val2 << endl;

    long double val3 = 0.2;
    cout << fixed << setprecision(128) << val3 << endl;

    __float80 val4 = 0.2;
    cout << fixed << setprecision(128) << val4 << endl;

    // __float128 val5 = 0.2;
    // cout << fixed << setprecision(128) << val5 << endl;

    return 0;
}

可能的输出为(g++ 12.1.0, with -std=c++17 -O2)：

1
2
3
4
0.20000000298023223876953125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0.20000000000000001110223024625156540423631668090820312500000000000000000000000000000000000000000000000000000000000000000000000000
0.20000000000000001110223024625156540423631668090820312500000000000000000000000000000000000000000000000000000000000000000000000000
0.20000000000000001110223024625156540423631668090820312500000000000000000000000000000000000000000000000000000000000000000000000000

你可以显式地注意到误差的存在，至于后面那一段 \(0\)，因为类型并不能支持到那么多位，所以 std::setprecision 就将其处理为 \(0\) 了。

刚才所说的可保证精度位数，就是十进制下，小数点前后，能够被精确表示的数位数量，这个其实可以通过 \(\dfrac{1}{2^{\text{bit}(M)}}\) 来计算。

以上主要介绍的是关于浮点数的底层实现，现在来说一点平常代码里会用到的细节。

对于两个浮点数的比较，我们不会直接使用 <, >，而是应当设置一个值 \(eps(ilon)\)。

什么意思？我们实际上是希望通过 \(eps\) 来消除误差对于比较大小带来的影响，我们希望比较的时候只考虑精度正确的那些位。

虽然你直接使用的时候可能没什么问题，但那是特殊情况，一般来说 \(eps\) 的取值需要取决于使用的类型以及需求。

一些实现：

1
2
3
4
5
6
7
const double eps = 1e-8; // 10^{-8}
inline int sgn(double x) {
    return (x > eps) - (x < -eps);
}
bool cmp(double x, double y) {
    return sgn(x - y) < 0; // less than
}

tips: 有的时候经过计算，可能答案会极度接近 0 但又不是，比如 -0.0000000001，这将导致输出 -0.0，这种情况加一个 \(eps\) 就能防止出错了，不过感觉不咋常用啊。

ref:

最后更新: November 6, 2023