编程语言中的整数和浮点数

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

整数

整数模型

编程语言中的整数类型不同于数学意义上的整数，而只是它的一个有限子集，因为计算机为了计算效率，会使用固定的字节数来存储一个整数数据，例如 $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$ 为负数。（这样的理解可以直接绕过通常教材中关于原码-反码-补码的繁琐定义，从模运算角度进行解释会更加显然且本质）

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

#include <iostream>
#include <limits>

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

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

    return 0;
}

程序运行结果为

4294967295
0
2147483647
-2147483648

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

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

编程语言支持

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

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

自动整除

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

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

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

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中被改过来了，只有使用//才会进行整除。

整除和取余

整除和取余是两个看起来非常简单明确的基本数学运算，但是在不同编程语言的实现中，其实存在着很多的差异，需要注意一下。

数论中的整除和取余

在数论中的整除和取余通常定义为：对于整数 $a,b \in \mathbb{Z}$，其中 $b \neq 0$，存在唯一的商 $q\in \mathbb{Z}$ 和余数 $r \in \mathbb{Z}$，使得
$$
a = b ,q + r
$$

对于余数要求 $0 \le r < b$，即余数是一个非负的且不超过除数的整数。

虽然初等数论主要关注的都是整数，通常不会涉及对实数的整除和取余运算，但是我们仍然可以直接将上述定义推广到实数中：对于实数 $a,b \in \mathbb{R}$，其中 $b \neq 0$，存在唯一的商 $q\in \mathbb{Z}$ 和余数 $r \in \mathbb{R}$，使得
$$
a = b ,q + r
$$

对于余数要求 $0 \le r < b$，即余数是一个非负的且不超过除数的实数。

数论中要求在任何情况下，商是一个整数，余数是一个非负的且绝对值不超过商的数。

计算机中的整除和取余

虽然在数论中已经提供了对于整数和实数在所有情况下的整除以及取余的定义，但是在计算机实践中并没有完全采纳这套定义，而是遵循了与之略微不同的规则，
在出现负数和出现浮点数的情况下，与数论中的定义会产生显著差异。
并且非常离谱的是，一共存在好几种候选的具体规则，不同的编程语言可能选择了不同的规则！

在计算机语言中，对于数 $a,b$（可能是整数或浮点数，这里不做区分），称满足如下要求的整数 $q$ 和数 $r$ 为商和余数
$$
a = b ,q + r, q \in \mathbb{Z}, |r| < |b|
$$

注意：这里的要求比数论中定义更弱，它没有要求 $r$ 非负。
事实上，这些要求不足以称为定义，因为在非整除的情况下，我们事实上可以找到两个满足要求的余数 $r_1 < 0 < r_2$，它们对应的商也不同，在数轴上这代表着在一个一般点两侧最靠近的两个整除点。到底选哪一个作为余数？我们必须明确说明这一点才能得到一个完整的定义。

考虑到对 $a,b \in \mathbb{N}$ 习惯用法（也就是数论中的定义）的兼容性，我们希望满足：在 $a,b>0$时，选择的余数 $r>0$。
但是在其它情况下呢？通常有如下几种策略：（这里按照余数的选取策略进行了分类，实际上也可以按照商的选取策略进行分类）

要求 $r$ 始终非负；
要求 $r$ 的符号始终与被除数 $a$ 一致；
要求 $r$ 的符号始终与除数 $b$ 一致；
要求 $r$ 选择 $r_1,r_2$ 中更靠近0的那个，也就是绝对值最小的那个。（这个并不能保证正整数整除的余数非负，通常只用于浮点数）

各种策略的示例如下表

取余策略	$4 \div 3$	$(-4) \div 3$	$4 \div (-3)$	$(-4) \div (-3)$
余数符号非负（数论）	商 1, 余数 1	商 -2, 余数 2	商 -1, 余数 1	商 2, 余数 2
余数符号与被除数相同（截断除法）	商 1, 余数 1	商 -1, 余数 -1	商 -1, 余数 1	商 1, 余数 -1
余数符号与除数相同（向下整除）	商 1, 余数 1	商 -2, 余数 2	商 -2, 余数 -2	商 1, 余数 -1
余数绝对值最小化（误差最小化）	商 1, 余数 1	商 -1, 余数 -1	商 -1, 余数 1	商 1, 余数 -1

注：这个表格无法展示第二个和第四个策略的区别，把4改成2再进行测试，就可以观察到两者确实是不同的策略。

这几种策略都有各自的出发点和合理性，因此哪一个也无法占据绝对的主导地位：

余数符号非负（数论）：保证余数符号非负的想法是最容易理解的，它是为了保持与数论中的定义一致。
余数符号与被除数相同（截断除法）：在很多编程语言中实现的除法实际上都是截断除法 $q = \text{trunc}(\frac{a}{b})$，这里的截断是选取可表示的数集中最靠近 $0$ 的数，即向 0 截断。这个做法导致我们始终会选择绝对值偏小的整数作为商，即选择靠近0的整除点来近似被除数，此时余数的符号就会和被除数保持一致。
余数符号与除数相同（向下整除）：如果我们将整除定义为始终向下截断：$q = \left[\frac{a}{b}\right]$，这会导致$q \le \frac{a}{b}$，此时余数的符号会和除数保持一致。
余数绝对值最小化（误差最小化）：这个可以保证 $b,q$ 作为 $a$ 的最优近似，保证 $b,q$ 是可表示的数集中最靠近 $a$ 的那个数。

值得注意的是，截断除法有一个独特的好处，在整除运算中可以直接将分子分母的负号提出来，满足
$$
(-a)/b = a/(-b) = -(a/b)
$$

这可能是它被编程语言普通采用的原因之一。

各种语言的实现

C/C++
- 对于两个整数之间的整除/和取余%，执行的策略都是截断除法，即保证余数符号与被除数相同，这也包括将其打包到一起的div函数。
- 对于浮点数之间则有两类实现：
  - fmod和std::fmod：保证余数符号与被除数相同；
  - remainder和std::remainder：保证余数的绝对值最小化。
MATLAB：常用的两个函数为
- mod：返回的余数保证符号与除数相同；
- rem：返回的余数保证符号与被除数相同。
Python：
- 整除//对应的策略是向下取整，%保证保证余数符号和除数一致（从而与//相匹配。
- 在math和numpy模块中提供了很多工具：
  - numpy.mod：行为与内置的%一致，即保证余数符号和除数一致；（numpy.remainder是它的别名）
  - math.fmod，numpy.fmod：调用的是C语言中的fmod实现，即保持余数符号和被除数一致。
Mathematica: 内置的 Mod函数采用的策略为：保证余数符号与除数一致。

例如考虑在不同语言中将浮点数 x 周期性平移到 $[0,T]$ 区间中，可以保证除数 $T > 0$，但是被除数却无法保证，我们希望得到的余数也是非负的，即保证余数的符号与除数相同。

在C/C++中，如果使用fmod，因为fmod只能保证余数符号和被除数一致，我们必须加一个补丁来保证范围

1
2
3

// [0,T]
x = std::fmod(x, T);
if (x < 0) {x += T;}

在MATLAB中，应该使用mod函数，例如

1 2	% [0,T] x = mod(x, T);

在Python中，应该使用numpy.mod函数，例如

1 2	# [0,T] x = np.mod(x, T)

小结

数学中的整除和取余的概念在计算机中并没有被完美地继承，而是根据具体情况采用了不同的实现，看起来比较混乱

截断除法可以保证余数的符号与被除数一致，并且运算中可以将符号提出，便于底层实现，因此在C/C++中被默认采用；
保证余数的符号与除数一致的做法更适合于科学计算，一些高级语言例如Python，MATLAB和Mathematica默认选择了这种策略，但是它们也为其它策略（主要是C/C++所提供的截断除法）提供了对应的接口。

浮点数

浮点数模型

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

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

其中：

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

一些非法的值：

$E$ 全为 1，$M$ 全为 0，代表正负无穷大 Inf（取决于符号位）；
$E$ 全为 1，$M$ 不全为 0，代表 NaN。

这部分参考长双精度类型 | Microsoft Learn。

单精度浮点数（32位浮点数）的长度为 4 个字节，在内存中具体分配为

其中

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

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

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

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

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

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

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

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

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

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

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

数值计算通常默认使用64位的浮点数，除非有特殊需求：

例如机器学习中对运算的精度其实并不关注，但是对运算效率非常敏感，可能使用32位的甚至更低的例如16位或8位的浮点数模型；
某些极端情况下，64位的浮点数精度无法达到要求，必须使用128位的浮点数。

由浮点数模型可知，对于固定位数的浮点数，所有可表示的数在实轴上并不是均匀分布的，而是关于 $0$ 对称，在 $0$ 附近的分布最密集，远离 $0$ 的位置则逐渐稀疏。只考虑 64 位浮点数模型，具体有：

在 $1$ 附近的两个相邻浮点数差值为 $2^{-52}$，约为 $10^{-16}$；
在 $x$ 附近的两个相邻浮点数差值约为 $|x| \times 10^{-16}$。

浮点数误差

浮点数与数学意义中的实数存在显著的差异，主要是浮点数的存储和运算规则导致的，例如浮点数存在上下界，浮点数运算不满足严格的加法和乘法的结合律和分配律等。
例如，对某些数值运算进行数学上等价的变形，可能导致结果在小数点后12位及以后产生细微的差异。

这种差异在绝大多数情况下只会产生机器精度的误差，并且不会随着运算的累计导致影响呈指数级增长（这还涉及到数值算法的稳定性）。
但是在某些情况下会非常明显，产生非常不合理的结果，下面举几个例子，只讨论最广泛使用的64位浮点数。（例子参考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++中，下面的两个整数定义语句是不等价的，第一个语句存在double到int的隐式类型转换

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

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

#include <cstdint>
#include <iostream>

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

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

    return 0;
}

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

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

N1 = 9e18 + 1 - 9e18 # 0.0

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

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

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>

严格来说，M_PI 属于 POSIX / glibc 扩展，在 Linux 下直接可用，但是 MSVC 默认不提供（虽然它也定义了），需要通过额外的选项来开启：

1 2	#define _USE_MATH_DEFINES #include <math.h>

对于 C++，cmath 作为 math,h 的封装，使用也基本类似。

基于宏的写法需要的位数实在太长了，而且 C++ 并不建议使用宏来定义常量，
我们可以利用反三角函数给出定义，它们的返回值也是足够精确的，例如

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;

Inf 和 NaN

最后记录一下关于几个特殊的浮点数值（nan，+inf，-inf），相关的运算规则。

Inf的产生通常有如下几种情景：

非零数除以零；
超过浮点数表示范围。

例如

result = 1.0 / 0.0; // +Inf
result = -1.0 / 0.0; // -Inf

result = 1e308 * 1e10; // +Inf
result = -1e308 * 1e10; // -Inf

NaN的产生通常有如下几种情景：

分子和分母同时为零，结果为NaN；
正无穷大与正无穷大相减，或负无穷大与负无穷大相减，结果为NaN；
不合法的数学运算，例如对负数计算平方根或对数等。

例如

1 2	result = 0.0 / 0.0; // NaN result = sqrt(-1.0); // NaN

与这些特殊值相关的运算规则太多了，懒得抄，不过大部分都是符合数学直觉的。
但是有一条值得注意：NaN代表的是非法状态，因此规定它与任何浮点数（包括自身）的比较结果都为假，即使和自身比较仍然是假，例如下面的这个函数可能返回false。

1
2
3

bool equal_self(double x){
    return (x == x);
}

浮点数优化问题

运算顺序问题

由于浮点数运算本身的运算性质非常糟糕，没有结合律和分配律，因此编译器如果尝试将对浮点数运算进行优化，例如大量重复的浮点数运算并行化，就可能出现与非并行版本结果不一致的情况，这类现象属于 浮点归约顺序变化 所带来的舍入误差。

例如 gcc 提供了 -ffast-math 选项，这个选项假定浮点数是满足结合律和分配律的，因此会启用一系列激进的浮点数优化，包括允许浮点数的舍入误差、忽略对 NaN 和无穷的处理、以及对运算重新排序等。
开启这个优化选项有时会显著提高计算速度，但是这可能导致浮点数结果不再严格符合 IEEE 标准。

在科学计算中，大量的操作都会涉及到对数组求和，但是浮点数的非结合特性导致了：分段求和再整体求和与直接整体求和的结果，可能不一样！当然，这种差异在大多数情况下并不显著。
所有对浮点数运算采用并行加速的计算框架，在特定运算中都会存在这类问题。

考虑下面的例子

import torch

x = torch.randn(100).abs()

z1 = x.sum().item()
z2 = x[:50].sum().item() + x[50:].sum().item()
print(f"{z1 - z2:.12f}") # != 0

实际运算中，z1 和 z2 可能并不会严格相等，也就验证了分组求和和整体求和的结果不一致。

重复索引问题

除此之外，还存在另一类问题：当程序包含对重复索引的写入更新时，多个更新可能同时作用于同一个目标位置。
此时问题不再只是“求和顺序不同”，而是涉及并发写入冲突、原子更新顺序不确定、以及高级索引写回语义等复杂因素，此时同样的输入在两次运行中也可能得到不同结果。

考虑下面的例子

import torch

x = torch.randn(9990).abs()
indices = torch.randint(0, 9990, (9000,))

def f(x, indices):
    y = x.clone()
    y[indices] += x[:9000]
    return y

y1 = f(x, indices)
y2 = f(x, indices)
print((y1 - y2).abs().max()) # != 0

注意：与前一个例子不同，这个例子必须要足够大的规模才能体现，而且结果可能非常大。

实际上，PyTorch 官方文档对这类问题的说法是：涉及到重复索引时，具体的行为未定义。在实际使用中应该改用专门的接口，例如

import torch

x = torch.randn(9990).abs()
indices = torch.randint(0, 9990, (9000,))

def f_safe(x, indices):
    y = x.clone()
    y.index_add_(0, indices, x[:9000])
    return y

y1 = f_safe(x, indices)
y2 = f_safe(x, indices)
print((y1 - y2).abs().max()) # != 0

补充

浮点数运算的特殊函数

在 C 语言标准库以及 MATLAB 中，可能会为了某些运算提供专门的函数，例如除了 $\ln x$，还会专门提供 $\ln(1+x)$ 的函数，甚至还有专门计算 $\sqrt{x^2+y^2}$ 的函数，这看起来有点奇怪，在逻辑上似乎是非必要的，但是这种设计有如下的两类原因：

计算精度不同：这些初等函数在不同取值范围内的数值计算误差不同，就像在不同位置进行泰勒展开一样，如果对某个区域的精度要求很高，可以使用对应的特化版本。
防止中间过程溢出：即使输入和精确输出都在可表示范围内，但是某些算法在中间过程中仍然有溢出的风向，通过一些额外的保护操作（例如归一化）可以避免这种情况。

例如专门计算 $\sqrt{x^2+y^2}$ 的hypot函数可以实现为如下算法

如果$x=y=0$，返回$0$；
否则，如果$|x| < |y|$，计算 $|y| \sqrt{1+(y/x)^2}$；
否则，$|x| \ge |y|$ 并且 $x \neq 0$，计算计算 $|x| \sqrt{1+(x/y)^2}$。

为了算法的健壮性，显然还需要对inf和nan等进行检查。

具体可以参考 John D. Cook 关于这些主题的两篇博文：
Math library functions that seem unnecessary，What’s so hard about finding a hypotenuse?。

累加误差导致的极小时间步

在数值计算中，浮点数的数值误差可能导致数值格式的时间推进的某个时间步仅有机器精度量级，这可能导致数值格式调阶。

考虑下面的例子

def update(tnow, dt):
    print(f"tnow = {tnow:.15f} dt = {dt:.12f}({dt:.6e})")


T = 1.0
tnow = 0.0
dt = 0.0125

while tnow < T:
    if tnow + dt >= T:
        dt = T - tnow

    update(tnow, dt)
    tnow = tnow + dt

运行结果为

tnow = 0.000000000000000 dt = 0.012500000000(1.250000e-02)

...

tnow = 0.987499999999998 dt = 0.012500000000(1.250000e-02)
tnow = 0.999999999999998 dt = 0.000000000000(1.554312e-15)

虽然这里的终止时刻恰好是时间步长的整数倍，但是实际的浮点数累和产生了误差，这导致最后一个时间步长只有机器精度量级。
这个极小的 $\Delta t$ 在某些情况下会造成显著影响：

对于多步法，这表明最后几步可能不是时间等距的，而且相邻时间步的差异非常大，多步法的格式通常假定时间步等距（不等距的格式在具体形式上会更加复杂），直接应用必然会掉阶。
对于单步格式，通常并不会出现问题
- 因为大部分单步格式都满足如下形式
  $$
  U_{n+1} = U_n + O(\Delta t)
  $$
  
  加上一个机器精度量级的时间迭代几乎没有影响。
- 但是存在例外情况：考虑线性常系数对流方程 $u_t + a u_x = 0$ 的 Lax-Friedrich 格式
  $$
  v_{j}^{n+1} = \frac{v_{j+1}^n+v_{j-1}^n}{2} - \frac{a}2 \frac{\Delta t}{\Delta x} (v_{j+1}^n-v_{j-1}^n)
  $$
  
  这是一个有条件相容的有限差分格式，局部截断误差为 $O(\Delta t + \frac{\Delta x^2}{\Delta t})$，只有在比例 $\frac{\Delta t}{\Delta x}$ 固定时才会具有一阶精度。在空间网格不变的情况下，如果存在极小的时间步 $\Delta t$，就必然影响 Lax-Friedrich 格式的精度。（同理其它有条件相容的格式也可能有同样的问题）

可以在时间循环中加上额外的检查来避免出现只有机器精度量级的最后时间步

eps = 1e-12
t = 0.0
while t < T - eps:
    dt = get_dt(data)

    if t + dt >= T - eps:
        dt = T - t

    update(data, t, dt)
    t += dt