0%

Recall to ALU/FPU 的硬件实现

继续有想法就多写点硬一些的内容。

关于计算机系统中的运算部件硬件实现,我脑子里面能想到的最近一次接触这些的时间可能还要追溯到在学校里面学数电那会了,而且大多数东西都忘得差不多了,刚好重新理一理并且顺便看一下浮点运算部件到底是怎么做的。

目前我们工作中并不会涉及到太深入到芯片设计方面的内容,因此本篇也只是对基础做一些整理,能帮助在脑子里面形成一个对不同运算部件的大体概念就不错了,充其量用这些在 MC 里面搞搞红石电路啥的。

篇尾还想找找 Nv 的硬件 spec 稍微对比分析下看看。

这里推荐一个各种不同运算算法在电路上的模拟网站:

通过上面的一些例子的模拟可以帮助更好地理解整个过程是怎么实现的。

整数运算 - ALU

先看看基础的二进制整数运算部分。

半加器(HA)

可以算是所有运算中最基础的部件了,功能是将两个输入 bit 加在一起,输出结果(Sum)和进位(Carry)。

a b Sum Carry
0 0 0 0
0 1 1 0
1 0 1 0
1 1 0 1

可以很容注意到,Sum 和 Carry 的输出可以用两个门电路来得到:

  • Sum = a XOR b
  • Carry = a AND b

Half Adder

全加器(FA)

多位二进制加法运算中考虑到每一位都有可能向前进位,因此对 3 个 bit 位进行加法操作的运算才是更常用到的,用两个 HA 串联在一起可以共同实现:

a b Last Carry Sum Carry
0 0 0 0 0
0 1 0 1 0
1 0 0 1 0
1 1 0 0 1
0 0 1 1 0
0 1 1 0 1
1 0 1 0 1
1 1 1 1 1

门电路实现:

  • s0 = a XOR b
  • c0 = a AND b
  • s1 = s0 XOR c
  • c1 = s0 AND c
  • Sum = s1
  • Carry = c0 OR c1

Full Adder

多位加法器

Ripple Carry Adder

逻辑上,通过任意数量的 FA 串联就能实现任意长度 bit 位数的加法运算,实际中需要结合硬件的物理情况进行考虑设计电路,例如 cmos 工艺、电路上的信号延迟等等。

如果直接用以上的基础实现进行位数规模扩展,整个芯片的效率会很低,因此在此基础上还有例如超前进位加法器、并行前缀进位加法器等等更加高效的实现。

乘法器

单 bit 位的乘法器只是一个 AND 的逻辑:

a b Prod
0 0 0
1 0 0
0 1 0
1 1 1

考虑多 bit 位的情况,运算过程其实与 10 进制下的乘法竖式计算是一致的,a、b 中的每一位对应相乘时需要附带一个偏移(其实也就是个加权一维卷积过程):

a1 a0
b1 b0
a1b0 a0b0
a1b1 a0b1
a1b1 a1b0 + a0b1 a0b0

$$\begin{aligned}
A &= a_1 * 2^1 + a_0 * 2^0 \\
B &= b_1 * 2^1 + b_0 * 2^0 \\
A*B &= a_1b_1 * 2^2 + (a_1b_0 + a_0b_1) * 2^1 + a_0b_0 * 2^0
\end{aligned}$$

因此我们可以得到一个通过不断移位和累加就能够完成的乘法器设计思路:

Sequential Multiplier

与多位加法器类似,如果直接按照最简单的硬件逻辑来实现电路的话的效率会很低,也非常不利于进行更多 bit 位的扩展,因此也衍生出了 Booth、Wallace树、Dadda、SRT 等等优化算法。

除法器

考虑 10 进制下对每一位进行试除和累积余数的竖式计算方法,我们可以得到一个类似的除法器设计思路:

Divider

恢复余数除法

这里来推一遍 10 / 3,在 8 位二进制下除法的过程:

Source A = 10 = (00001010)2
D = 3 = (00000011)2
Init D 初始时左移 A 的 bits 变成 00110000
Step 1 左移 A -> 00010100,A = A - D < 0 则该位结果为 0,A 保留余数 00010100 + 0
Step 2 左移 A -> 00101000,A = A - D < 0 则该位结果为 0,A 保留余数 00101000 + 0
Step 3 左移 A -> 01010000,A = A - D = 00100000 >= 0 则该位结果为 1,A 保留余数 00100000 + 1 -> 00100001
Step 4 左移 A -> 01000010,A = A - D =00010010 >= 0 则该位结果为 1,A 保留余数 00010010 + 1 -> 00010011

最终得到的结果为:
R = A[4:7] = (0001)2 = 1
Q = A[0:3] = (0011)2 = 3

由于初始计算时 D 本身已经偏移了 A 的 bits 数,因此 A - D 的正负性不回受后面补上的除完结果影响,因此这里为了简化一个寄存器,结果的 Q 直接更新在 A 中了。最终计算截断前 4 位为余数,后4位为结果。如果要继续往后面推算小数部分,则 Q 应该独立更新。

为什么这个算法要叫做恢复余数算法呢?
注意到上述计算中判断 A = A - D < 0 的这一步,事实上这里需要先发生减法运算,对结果进行 test 判断小于 0 之后,再把 D 加回去,恢复出 A,因而被称为恢复余数算法。

不恢复余数除法

那么 A - D < 0 的时候能不能不用恢复到 A 再往下走?答案是可以的,相比于不恢复余数法整体上算法复杂一些,但是在硬件实现上更有优势。
这里每产生一个商的比特位只需要一次加或减操作,核心是在减法操作后不需要进行余数恢复,这样就使得执行的速度更快了。

继续以 10 / 3,在 8 位带符号二进制下除法的过程举例:

Source A = 10 = (000001010)2
D = 3 = (000000011)2
-D = -3 = (100001101)2
Init D 初始时左移 A 的 bits 变成 000110000
-D 左移变成 111010000
Step 1 左移 A -> 000010100,A 是正数,A = A - D = 111100100 < 0 则该位结果为 0,A 保留余数 111100100 + 0
Step 2 左移 A -> 111001000,A 是负数,A = A + D = 111111000 < 0 则该位结果为 0,A 保留余数 111111000 + 0
Step 3 左移 A -> 111110000,A 是负数,A = A + D = 000100000 >=0 则该位结果为 1,A 保留余数 000100000+1-> 000100001
Step 4 左移 A -> 001000010,A 是正数,A = A - D = 000010010 >=0 则该位结果为1,A保留余数 000010010+1 -> 000010011

最终得到的结果为:
R = A[4:7] = (0001)2 = 1
Q = A[0:3] = (0011)2 = 3

再往后还有 SRT、高基除法等等更多性能更好或者更适合电路实现的变种和改进算法,这里不再展开。

浮点数运算 - FPU

浮点数的表示

首先看一下将小数部分直接进行二进制转换的表示法:

$$
(00001100.00001100)_2 = 2^3 + 2^2 + 2^{-5} + 2^{-6} = 12.046875
$$

用 2 的幂次直接来指代小数部分可以提供的范围和精度都非常有限,尤其截断之后产生的误差会相当大,因而一般通常大家讲到定点数时一般都用来指代整数了,比较少有场景直接这么用。

对于以下 x、y 两个数,数值范围上差别很大,但如果写成科学计数法的方式来表示,则写法上非常相似:

$$\begin{aligned}
x &= (00000000.00001001)_2=(1.001) * 2^{-6} \\
y &= (10010000.00000000)_2=(1.001) * 2^7
\end{aligned}$$

并且我们可以发现,用这种方式可以同时兼顾较大的数值范围和较高的精度。

在 IEEE 754 标准中,将浮点数的表示为:1 个符号位 + e 个指数位 + f 个尾数位的形式。

以 float32 为例,规定为 1/8/23 的格式。指数位有 8 位,是一个 0~255 的整数,用 Exp 移位之后的结果来表示,其中 0 和 255 这两个值被用来作为特殊标记,尾数有 23 位:

E + Bias Frac
表示 0 0 0
表示低于精度下限的非正规数
$0.f * 2^E$
0 非 0
Inf 255 0
NaN 255 非 0
常规浮点数
$1.f * 2^E$
[1, 254]
E -> [-126, 127]
任意数

因此 float32 能够表示的绝对值最小的规格化数为($2 ^ {-126}$),绝对值最大值小于($2^{128}$)。如果算上非规格化数,能表示的最小绝对值还能进一步到 $2^{-126} * 2^{-23} = 2^{-149}$。

Sign Exp Frac Max Value Min Value Normal Min Value Subnormal
float16 1 5 10 $2^{15} * \sum^{0}_{i=-10}2^i = 65504$ $2^{-14} = 6.10 * 10^{−5}$ $2^{-14} * 2^{-10} = 5.96 * 10^{−8}$
float32 1 8 23 $< 2 ^ {128} \approx 3.4 * 10^{38}$ $2^{-126} \approx 1.18 * 10^{−38}$ $2^{-126} * 2^{-23} \approx 1.4 * 10^{−45}$
float64 1 11 52 $< 2 ^ {1024} \approx 1.80 * 10^{308}$ $2^{-1022} \approx 2.23 * 10^{−308}$ $2^{-1022} * 2^{-52} \approx 4.94 * 10^{−324}$

浮点数的加减法

浮点数的加减法比整数要麻烦不少,因为两个操作数的指数往往不一样,尾数不能直接相加。整个过程大致可以拆成 对阶 -> 尾数加减 -> 规格化 -> 舍入 这几步(以 $x = m_x * 2^{e_x}$、$y = m_y * 2^{e_y}$ 为例):

  1. 对阶:比较两个数的指数,把指数较小的那个数的尾数右移 $|e_x - e_y|$ 位,让两者的指数对齐到较大的那个。右移会把一些低位挤出尾数的表示范围,这些被挤出去的位并不能直接丢掉,需要用后面提到的 Guard/Round/Sticky 位暂存下来用于舍入。
  2. 尾数加减:根据两个数的符号,对对齐后的尾数做加法或减法。
  3. 规格化:加减之后的结果不一定还是 $1.f$ 的规格形式。
    • 如果相加产生了进位溢出(结果落到 $[2, 4)$),需要把尾数右移 1 位、指数加 1;
    • 如果是相近的两个数相减,会出现大量前导 0(即所谓的 subtractive cancellation / 大量抵消),此时需要把尾数左移 $n$ 位、指数减 $n$,这一步需要一个前导零预测/计数器(Leading Zero Anticipator / Counter)来快速定位。
  4. 舍入:根据暂存的 G/R/S 位和当前的舍入模式决定最终尾数(见下一节)。舍入本身也有可能再次产生进位溢出,因此可能还要再补一次规格化。

举一个简单的二进制例子,计算 $1.5 + 0.5$:

$$\begin{aligned}
x &= 1.5 = (1.1)_2 * 2^0 \\
y &= 0.5 = (1.0)_2 * 2^{-1} \\
\text{对阶后} \quad y &\to (0.1)_2 * 2^0 \\
\text{相加} \quad x + y &= (10.0)_2 * 2^0 \\
\text{规格化} &\to (1.0)_2 * 2^1 = 2.0
\end{aligned}$$

可以看到,即便是最简单的加法,硬件上也需要移位器、加法器、前导零计数器和舍入逻辑一起配合,比整数加法复杂了不少。

FP Add/Sub

舍入

上一节多次提到舍入,这里单独展开。由于尾数位数有限,运算的中间结果往往无法被精确表示,需要按规则舍入回可表示的浮点数。IEEE 754 一共定义了 5 种舍入模式:

  • 就近舍入、向偶数取整(roundTiesToEven):舍入到最接近的可表示值,如果恰好落在两个值正中间,则选择尾数最低位为偶数的那个。这是默认模式,好处是能避免舍入误差的系统性偏移。
  • 就近舍入、远离零取整(roundTiesToAway):正中间时向绝对值更大的方向舍入,主要用于十进制。
  • 向 $+\infty$ 舍入(roundTowardPositive,向上取整)
  • 向 $-\infty$ 舍入(roundTowardNegative,向下取整)
  • 向 $0$ 舍入(roundTowardZero,直接截断)

为了在硬件里正确地实现舍入,仅保留尾数是不够的,还需要在尾数最低位(LSB)之后额外维护 3 个标志位:

  • Guard(G):LSB 后的第 1 位。
  • Round(R):LSB 后的第 2 位。
  • Sticky(S):LSB 后剩余所有位的”或”,只要还有任何一个 1,S 就是 1。

在对阶右移、乘法产生双倍位宽结果等场景中,那些被移出精度范围的信息就被压缩进了这 3 个位里。以默认的向偶数取整模式为例,是否要给尾数进位(+1)的判断如下($L$ 为当前保留尾数的最低位):

G R S 与 0.5 ULP 的关系 操作
0 x x < 0.5 ULP 截断(不进位)
1 0 0 = 0.5 ULP(正中间) 向偶:$L = 1$ 时才进位
1 0 1 > 0.5 ULP 进位
1 1 x > 0.5 ULP 进位

也就是说进位条件可以写成 $G \cdot (R + S + L)$。有了 Sticky 位,硬件就不需要保留被移出的全部低位,只用一个”或”累积起来的标志位就能保证舍入的正确性。

浮点数的乘法

相比加减法,浮点乘法的流程反而更直接,因为不需要对阶:

  1. 符号:结果符号 = 两个操作数符号位异或,$s = s_x \oplus s_y$。
  2. 指数:把两个指数相加。注意由于两个指数各自都带了一个 Bias 偏移,直接相加会多出一个 Bias,需要减回来:$e = e_x + e_y - \text{Bias}$。
  3. 尾数:把两个带隐含位的尾数 $1.f_x$ 与 $1.f_y$ 相乘。对 float32 来说就是两个 24 bit 数相乘得到一个 48 bit 的结果——这里正好可以复用前面 ALU 部分讲到的整数乘法器(Booth、Wallace 树等)。
  4. 规格化 + 舍入:两个 $[1, 2)$ 的尾数相乘,结果落在 $[1, 4)$,如果 $\ge 2$ 就右移 1 位、指数加 1,然后按舍入规则处理低位。
  5. 最后再处理上溢/下溢(转成 Inf 或非规格化数)以及 NaN 等特殊值。

FP Multiplication

浮点数的除法

除法的符号和指数处理与乘法对称:符号仍是异或,指数变成相减 $e = e_x - e_y + \text{Bias}$,核心难点在尾数相除 $1.f_x / 1.f_y$。硬件上主要有两大类实现思路:

  • 数字递归(digit recurrence):每次迭代求出商的一位,本质上就是前面整数除法里讲的不恢复余数法(SRT 是它的高基改进版,通过查表一次求出多位商)。特点是电路面积小、还能顺带得到精确余数,但延迟随位数线性增长,比较慢。

  • 乘法迭代(multiplicative):先求出除数的倒数 $1/y$,再乘以被除数。常见的有:

    • Newton-Raphson:迭代式 $x_{n+1} = x_n (2 - y \cdot x_n)$,从查找表取一个初始估计值后,每迭代一次有效位数翻倍(二次收敛)。
    • Goldschmidt:同样是二次收敛,但把两个乘法拆开使其可以并行流水,更适合硬件实现,GPU 上就常用这一类。

    这类方法迭代次数少,但每次迭代都要用到乘法器,且最后往往还需要一步修正/舍入才能得到正确舍入的结果。

FP Division

平方根 $\sqrt{x}$ 的实现思路与除法类似,既可以用数字递归,也可以先用 Newton-Raphson 求 $1/\sqrt{x}$ 再做修正。

融合乘加(FMA)

FMA(Fused Multiply-Add)是现代 CPU/GPU 里非常关键的一条指令,它一次计算 $a \times b + c$,并且只在最后做一次舍入 $\text{rn}(a \times b + c)$。而如果拆成一次乘法加一次加法,就会有两次舍入 $\text{rn}(\text{rn}(a \times b) + c)$。IEEE 754-2008 把 FMA 正式纳入了标准。

只做一次舍入带来的好处不只是省一次舍入误差,更重要的是它在乘法阶段保留了双倍位宽的完整乘积,因此在后续做减法即使发生大量抵消也不会损失精度。用 NVIDIA 白皮书里的例子,计算 $x^2 - 1$,取 $x = 1.0008$、小数点后保留 4 位十进制:

$$\begin{aligned}
x^2 - 1 &= 1.60064 * 10^{-4} \quad (\text{精确值}) \\
\text{rn}(x^2 - 1) &= 1.6006 * 10^{-4} \quad (\text{FMA,单次舍入}) \\
\text{rn}(\text{rn}(x^2) - 1) &= 1.6000 * 10^{-4} \quad (\text{先乘后加,两次舍入})
\end{aligned}$$

可以看到 FMA 的结果 $1.6006$ 明显更接近精确值,而先乘后加因为中间把 $x^2 = 1.00160064$ 舍入成了 $1.0016$,把后面对结果有贡献的位给丢掉了。

正因为这些特性,FMA 成了很多运算的基础构件:前面除法里的 Newton-Raphson 迭代、向量点积、多项式的 Horner 求值等,都能用一串 FMA 高效且更精确地完成。这也是为什么 GPU 里的 FP32/FP64 计算单元基本都是以 FMA 为核心来设计的。


GPU 上的浮点运算硬件

前面整理的都是通用的运算部件原理,最后按照开头说的,结合 NVIDIA 的硬件 spec 来看看这些东西在 GPU 上到底长什么样。

IEEE 754 合规性的演进

早期的 NVIDIA GPU 出于面积和性能考虑,在浮点合规性上是有不少妥协的:

  • Tesla 架构(CC 1.x):FP32 不支持非规格化数(直接 flush-to-zero),除法和平方根是基于 SFU 求倒数的近似实现、并非正确舍入,乘加也只是截断的 MAD 而不是真正的 FMA。直到 CC 1.3 才补上了完整符合 IEEE 754 的 FP64(含 FMA)。
  • Fermi 架构(CC 2.0,2010):第一次做到 FP32 和 FP64 全面符合 IEEE 754-2008——加减乘除、开方、FMA 都是正确舍入,支持全速的非规格化数和全部舍入模式。此后这一基线基本保持了下来。

精确通路 vs 快速通路

GPU 上一条浮点通路往往同时提供”精确”和”快速”两种选择:

  • 计算核心(CUDA Core) 里以 FMA 为核心,负责正确舍入的加减乘除。以 FP32 除法为例,实际做法是先用 SFU 取一个倒数近似作为种子,再用 Newton-Raphson 迭代把精度提升到 24 bit 并修正,从而默认得到正确舍入的结果。
  • 特殊功能单元(SFU,对应指令 MUFU.*) 用分段二次插值直接给出 $1/x$、$1/\sqrt{x}$、$2^x$、$\log_2 x$、$\sin$、$\cos$ 等的近似值,精度大约 22~24 bit(几个 ULP),但吞吐极高。CUDA 里的 rsqrtf__sinf__expf__fdividef 这些 intrinsic 走的就是这条路。

开启 -use_fast_math 后,编译器会把除法、超越函数等换成上面这些 SFU 近似,并且把非规格化数直接 flush-to-zero——速度更快,但不再是正确舍入、也不再严格符合 IEEE 754。是否使用取决于业务对精度和可复现性的要求。

低精度格式与 TensorCore

深度学习的兴起让 GPU 引入了一系列低精度浮点格式,核心思路都是在”范围”和”精度”之间做不同的取舍:

格式 Sign/Exp/Frac 特点与用途
FP16 1/5/10 范围偏小、精度尚可,推理常用
BF16 1/8/7 指数位与 FP32 一致,牺牲精度换范围,训练友好
TF32 1/8/10 Ampere TensorCore 的输入格式,共 19 个有效位
FP8 (E4M3) 1/4/3 Hopper,偏重精度
FP8 (E5M2) 1/5/2 Hopper,偏重动态范围
FP6 (E3M2) 1/3/2 Blackwell
FP4 (E2M1) 1/2/1 Blackwell,极致压缩用于大模型推理

TensorCore 做的其实就是矩阵版的 FMA:$D = A \times B + C$,通常用低精度做乘法、再用更高精度(一般是 FP32)做累加,从而在保证一定数值稳定性的前提下大幅提升算力密度。

一个容易踩坑的点:结合律

最后值得一提的是,浮点加法不满足结合律,$(a + b) + c$ 和 $a + (b + c)$ 的结果可能不一样(NVIDIA 白皮书里就有这样的例子,两者都符合 IEEE 754,但结果就是不同)。这在 GPU 上尤其常见:并行归约(reduction)时求和顺序不固定,同一份数据多跑几次、或者换一块卡,结果就可能在最低位上出现差异。如果业务需要严格的逐位可复现,就得固定运算顺序、并控制编译器对 FMA 的合并行为。


至此,从半加器一路到 GPU 上的 TensorCore,运算部件的大体脉络就算理清楚了。