整理一下编程语言中的整数和浮点数的相关内容,针对的情景是科学计算。

整数

整数模型

编程语言中的整数类型不同于数学意义上的整数,而只是它的一个有限子集,因为计算机为了计算效率,会使用固定的字节数来存储一个整数数据,例如 \(n\) 个字节,这意味这只有 \(2^{8n}\) 个不同状态,只能表示 $2^{8n} $ 个整数。

\(n\) 个字节所对应的 \(8n\) 比特的值依次记作 \(a_i \in \{0,1\}\),这里 \(i=0,\dots,8n-1\),那么通常有两类方案:

第一种方案是无符号整数,表示的值 \(V\)\[ V = 2^0 a_0 + 2^1 a_1 + \dots + 2^{8n-2} a_{8n-2} + 2^{8n-1} a_{8n-1} \] 表示的范围为 \[ [0,2^{8n}-1] \]

第二种方案是有符号整数,表示的值 \(V\)\[ V = 2^0 a_0 + 2^1 a_1 + \dots + 2^{8n-2} a_{8n-2} - 2^{8n-1} a_{8n-1} \] 表示的范围为 \[ [-2^{8n-1},2^{8n-1}-1] \] 注意这里为了调整范围,让最高位 \(a_{8n-1}\) 扮演了不同的角色,可以将其称为符号位:取0时 \(V\) 为正数,取1时 \(V\) 为负数。(这里直接绕过通常教材中关于原码-反码-补码的繁琐定义,从模运算理解会更加显然且本质)

不管哪一种方案,都存在有限的上下界,因此执行加减乘除运算时,都存在溢出的风险。 计算机通常会对运算的结果取模,将其重新映射到可表示范围中。 例如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <iostream>
#include <limits>

int main() {
unsigned int maxUnsignedInt = std::numeric_limits<unsigned int>::max();
std::cout << maxUnsignedInt << std::endl;
unsigned int overflowUnsigned = maxUnsignedInt + 1;
std::cout << overflowUnsigned << std::endl; // 0

int maxInt = std::numeric_limits<int>::max();
std::cout << maxInt << std::endl;
int overflowSigned = maxInt + 1;
std::cout << overflowSigned << std::endl; // 输出未定义(可能为负值)

return 0;
}

程序运行结果为

1
2
3
4
4294967295
0
2147483647
-2147483648

虽然实际都会进行模运算,但是在语法上对两种情况的规定可能是不一样的:

  • 对无符号整数的运算结果取模是C/C++语法标准规定的;
  • 对有符号整数的运算结果取模实际上是未定义行为。

编程语言支持

常见编程语言的整数支持如下:

  • C/C++:
    • int通常为32位整数;
    • 不建议使用long,因为它只是表示范围不小于int的整数类型,但是在不同平台的实现有差异:在Windows中是64位,但是在Linux中是32位;
    • long long通常为64位整数;
    • 可以使用int32_tint64_t等固定字节数的整数类型;
    • 不建议使用int8_tuint8_t,因为它们本质上是charunsigned char的别名,在输入输出时可能会被当作字符而非数值处理。
  • MATLAB:
    • MATLAB提供了固定字节的整数类型,例如int32, uint64等,但是不建议使用。
  • Python:
    • Python的整数实现不是基于固定字节的,而是一种动态的模型,可以支持任意大的整数,在内存允许的情况下不存在溢出问题,代价就是性能偏低;
    • Numpy提供了基于固定字节的整数类型,例如numpy.int32,可以带来更快的计算效率。
  • Fortran:
    • integer通常为32位整数;
    • 可以使用integer(kind=4)integer(kind=8) 指定不同的位数。

自动整除

在很多编程语言中,两个整数的除法默认会进行整除,只会得到一个整数值,例如在C语言中

1
2
3
4
int a = 1;
int b = 2;
double c = a / b; // c = 0.000000000000
double d = 2 / 3; // d = 0.000000000000

只有在相除的两个数的类型不全为整数时,才会进行通常意义下的除法,可以使用下面的做法来避免整除

1
2
3
4
int a = 1;
int b = 2;
double c = (a * 1.0) / b; // c = 0.500000000000
double d = 2.0 / 3; // d = 0.666666666666

在Python中也存在这种问题,例如Python2的/在某些情况下也会自动整除,但是在Python3中被改过来了,只有使用//才会进行整除。

浮点数

浮点数模型

科学计算中主要使用的数据类型都是满足IEEE标准的浮点数类型,符合 IEEE 规范的浮点数采用如下的二进制科学计数法表示,分别用固定的长度存储底数和指数,以及符号位三部分。

\[ V = (-1)^s \times (1 + M) \times 2^E \]

其中:

  1. \((-1)^s\) 代表符号位,\(s=0,1\)
  2. \(1+M\) 代表底数(存储时缺省了前面的 1),范围 \(0 \le M <1\)
  3. \(E\) 代表指数,可以有正负。

一些非法的值:

  1. \(E\) 全为 1,\(M\) 全为 0,代表正负无穷大 Inf(取决于符号位);
  2. \(E\) 全为 1,\(M\) 不全为 0,代表 NaN。

这部分参考微软的文档

单精度浮点数的长度为 4 个字节,在内存中具体分配为

其中

  1. 符号位使用 1 比特,记录 \(s=0,1\)
  2. 指数部分使用 8 比特,代表 \(E \in [-2^{7},2^{7}-1]\)
  3. 底数部分使用 23 比特,代表二进制的小数部分 \(M\),前面缺省一个 1,即底数 \(1+M\)

单精度浮点数可表示的最大值约为

\[ 2 \times 2^{127} = 2^{128} \approx 3.4 \times 10^{38} \]

双精度浮点数,长度为 8 个字节,在内存中的具体分配类似单精度浮点数,其中

  1. 符号位使用 1 比特,记录 \(s=0,1\)
  2. 指数部分使用 11 比特,代表 \(E \in [-2^{10},2^{10}-1]\)
  3. 底数部分使用 52 比特,代表二进制的小数部分 \(M\),底数 \(1+M\)

双精度浮点数可表示的最大值约为

\[ 2 \times 2^{1023} = 2^{1024} \approx 1.8 \times 10^{308} \]

四精度浮点数,长度为 16 个字节,在内存中的具体分配类似单精度浮点数,其中

  1. 符号位使用 1 比特,记录 \(s=0,1\)
  2. 指数部分使用 15 比特,代表 \(E \in [-2^{14},2^{14}-1]\)
  3. 底数部分使用 112 比特,代表二进制的小数部分 \(M\),底数 \(1+M\)

四精度浮点数可表示的最大值约为

\[ 2 \times 2^{16383} = 2^{16384} \approx 1.2 \times 10^{4932} \]

浮点数误差

浮点数与数学意义中的实数存在显著的差异,主要是浮点数的存储和运算规则导致的,例如浮点数存在上下界,浮点数运算不满足严格的加法和乘法的结合律和分配律等。

这种差异在大多数情况下只会产生机器精度的误差,但是在某些情况下会非常明显,产生非常不合理的结果,下面举几个例子。(例子参考MATLAB的官方文档)

舍入误差:浮点数的有限精度表示可能导致舍入误差。例如4/3不能精确表示为二进制分数,在存储和运算中都会被截断。

1
a = 1 - 3*(4/3 - 1) % 2.2204e-16

淹没:将两个量级差异非常大的数进行直接加减,可能会直接淹没量级过小的数据。

1
b = 1 + 1e-16 % 1

抵消:将两个非常相近的数进行直接相减时,得到的结果可能会损失很多有效位数,甚至得到错误的结果。

1
2
3
c1 = (100 + 1e-14) - (100 - 1e-14) % 2.8422e-14
c2 = (200 + 1e-14) - (200 - 1e-14) % 0
c3 = (2.0^53 + 1) - 2.0^53 % 0

这组结果表明

  • 第一个结果的相对误差已经接近百分之五十,有效位数明显降低。
  • 后两个结果更加极端,由于在对应量级附近的两个浮点数差异已经大于字面量的差异,因此实际处理时将其舍入为同一个浮点数进行相减。

在数值微分公式中,我们仍然会对一些非常相近的数进行直接相减,虽然浮点数模型表明这种做法会损失精度,但是这也是没有办法的做法。

上面的这些事实意味着,在编写数值计算程序时,即使我们将相同的数学公式以略微不同的语法翻译到代码中,也有可能产生不同的结果! 在绝大多数情况下,这种数值差异可能只是机器精度量级的,并不会造成实质性的影响,但是在某些极端情况下,也可能得到非常不合理的结果。

编程语言支持

常见编程语言的浮点数支持如下:

  • C/C++:
    • 通常使用64位的浮点数double,浮点数字面量默认均为double类型;
    • 不要使用32位的浮点数float,它的精度对于科学计算而言是不够的;
    • 在语法上还支持更高精度的long double,但是语法标准只规定了它的精度不低于double,并没有规定到底是多少位。实际上,long double在不同编译环境下的实际位数是不一样的,可能不是IEEE标准的128位浮点数,在跨平台时容易出现问题。
  • MATLAB:
    • 默认情况下,数和矩阵元素等都是64位浮点数double,不需要特别处理。
  • Python:
    • 默认的浮点数类型虽然名称为float,但是实际是64位浮点数;
    • 某些第三方库可以提供更高精度的浮点数。
  • Fortran:
    • 通常使用64位浮点数real
    • 在某些环境中还支持更高的128位浮点数。

浮点数字面量

需要注意的是,编程语言对于数值字面量的解释规则与数学直觉不同,使用科学计数法(含e)的字面量总会被解释为双精度浮点数,无论它自身有没有小数部分,这是很多编程语言普遍采用的规则。

例如在C++中,下面的两个整数定义语句是不等价的,第一个语句存在doubleint的隐式类型转换

1
2
int n1 = 1e6;
int n2 = 1000000;

在通常情况下两者不会产生差异,但是在极端情况下仍然存在着两类隐患:浮点数误差和整数范围溢出。 其中浮点数误差的问题可能是比较隐蔽的,例如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <cstdint>
#include <iostream>

int main(int argc, char *argv[]) {
int64_t N1 = 9e18 + 1 - 9e18;
std::cout << N1 << std::endl; // 0

int64_t tmp = 9e18;
std::cout << tmp << std::endl; // 9000000000000000000
int64_t N2 = tmp + 1 - tmp;
std::cout << N2 << std::endl; // 1

return 0;
}

这里的大整数并不存在溢出问题,9e18仍然在int64_t的可表示范围内。 两者的区别在于:N1的定义右侧执行的是浮点数的加减法,N2的定义右侧执行的是整数的加减法。

在其它语言中也会出现同样的问题,例如在Python中

1
2
3
4
N1 = 9e18 + 1 - 9e18 # 0.0

tmp = int(9e18)
N2 = tmp + 1 - tmp # 1

在C++中还存在doublefloat字面量的区别,浮点数字面量在无后缀时视作double,在含有f后缀时视作float,例如下面的几个等式判断

1
2
3
4
5
6
7
float a1 = 2.2;
bool b1 = (a1 == 2.2); // false
bool c1 = (a1 == 2.2f); // true

float a2 = 2.25;
bool b2 = (a2 == 2.25); // true
bool c2 = (a2 == 2.25f); // true

这里第一个结果是因为2.2在转换为float时存在舍入误差,最后两个结果则是因为2.25恰好可以无损转换为float

浮点数常量

在数值实验中使用一些浮点数常量是很常见的需求,但是实际定义和使用浮点数常量时存在一些小坑,以pi为例进行说明。

很多语言例如MATLAB和Python都提供了pi可以直接使用,但是在C/C++中还需要我们自行实现。 下面这种随手给出的定义方式是不行的,它会引入远大于机器精度的截断误差

1
#define PI 3.1415926

我们可以参考官方的实现:在<math.h>头文件中,虽然pi也是通过宏直接定义的, 但是它使用了位数足够长的字面量,确保表达式的截断误差小于机器精度

1
#define M_PI 3.14159265358979323846 // <math.h>

上面的正确写法需要的位数实在太长了,而且通常情况下都不建议使用宏来定义常量, 我们可以利用反三角函数给出定义,它们的返回值也是足够精确的,例如

1
2
3
const double pi = acos(-1.0);
// or
const double pi = 4 * atan(1.0);

在C++20之后,我们还可以直接使用<numbers>提供的相关数学常量

1
constexpr double pi = std::numbers::pi;