Skip to content

功能单元:浮点运算

硬件描述 1 — 4.75亿美元的Bug

1994年,Thomas Nicely教授在计算布朗常数时发现Intel Pentium处理器的浮点除法单元存在错误——SRT商选择表中5个条目被遗漏,导致某些除法运算的结果在第9至第10位十进制有效数字处出错。Intel最初试图淡化此问题,声称"普通用户270亿年才会遇到一次"。但公众反应剧烈,Intel最终宣布无条件更换所有有缺陷的处理器,预提4.75亿美元准备金。这是处理器历史上最昂贵的硬件缺陷。FDIV Bug的直接后果是:浮点单元的形式化验证从此成为行业标准——Intel、AMD、ARM均在浮点除法器和平方根单元上投入了大量形式验证资源,确保商选择表的每一个条目都经过数学证明。

设计提示

统一视角。处理器设计的本质是在有限的晶体管预算和功耗约束下,通过投机和并行的层层叠加来逼近指令吞吐率的理论上限。FMA(融合乘加)正是功能单元级"并行"的体现——它将乘法和加法合并为单一操作,只执行一次舍入而非两次,在精度和延迟上都优于分离的乘法+加法。在AI工作负载中,FMA是最核心的指令——矩阵乘法的每个元素计算都是一次乘加操作。Apple M4的浮点单元每周期可以执行4次双精度FMA(32 FLOPS/cycle),这背后是整数ALU作为浮点尾数运算基础(回顾第 30.0 章的加法器和乘法器设计)的直接延伸。浮点单元虽然只占核心面积的5%\sim10%(参见第 3.0 章中面积分析),但其设计复杂度和验证难度却不成比例地高。

第 30.0 章中,我们详细探讨了整数功能单元的设计——从加法器、移位器到乘法器和除法器。整数运算的操作数格式是直观的:一个NN位二进制补码数直接表示一个整数值,运算规则明确且高效。然而,科学计算、图形渲染、机器学习和信号处理等领域的工作负载大量使用浮点运算(Floating-Point Arithmetic),浮点数的表示方式和运算规则与整数截然不同——符号位、指数、尾数的三段式编码引入了对阶、规格化、舍入等整数运算中不存在的步骤,使得浮点功能单元的设计复杂度远高于整数单元。

本章首先回顾IEEE 754浮点标准的核心内容,然后逐一剖析浮点加法器、浮点乘法器、融合乘加单元(FMA)以及浮点除法与平方根单元的微架构设计。最后讨论浮点异常处理的硬件实现。本章的重点不在于浮点算术的数学理论,而在于如何将这些数学运算高效地映射到流水线化的硬件结构上——这正是处理器设计师面临的核心挑战。

在现代超标量处理器中,浮点单元通常被组织为独立于整数单元的执行域(domain),拥有自己的物理寄存器文件、发射队列和执行端口。表表 31.1概览了现代处理器中浮点功能单元的典型配置。

运算延迟吞吐量流水线典型面积占比
浮点加法 (FADD)3\sim4周期1/周期全流水线\sim5%
浮点乘法 (FMUL)4\sim5周期1/周期全流水线\sim15%
融合乘加 (FMA)4\sim5周期1/周期全流水线\sim25%
浮点除法 (FDIV)13\sim21周期1/4\sim8周期非全流水线\sim8%
浮点平方根 (FSQRT)15\sim22周期1/4\sim8周期非全流水线共享FDIV

现代处理器浮点功能单元的典型配置

IEEE 754浮点标准

IEEE 754标准(最新版本为IEEE 754-2019)定义了浮点数的表示格式、运算语义、舍入规则和异常处理机制。几乎所有现代处理器的浮点单元都严格遵循这一标准,确保不同硬件平台上的浮点计算结果完全一致——这种跨平台一致性对于科学计算的可重复性和软件的可移植性至关重要。理解IEEE 754标准是设计浮点功能单元的前提。

IEEE 754标准的核心设计哲学是:在有限的位宽下,尽可能保证每次运算的结果都是"最好的可能值"——即运算的精确结果在目标格式可表示范围内的最接近值(正确舍入,correctly rounded)。这个看似简单的要求对硬件实现有深远的影响:为了确定"最接近的可表示值",硬件需要计算出比目标精度更多位数的中间结果,然后通过舍入逻辑来选择最终结果。

浮点数的表示

一个IEEE 754浮点数由三个字段组成:

v=(1)S×1.M×2Ebias v = (-1)^{S} \times 1.M \times 2^{E - \text{bias}}

其中SS符号位(Sign),EE指数(Exponent),MM尾数(Mantissa,也称Significand或Fraction)。对于规格化数(normalized number),尾数的整数部分隐含为1(即"hidden bit"),因此实际有效数字为1.M1.M

IEEE 754定义了多种精度格式,其中最常用的是单精度(binary32)和双精度(binary64)。图图 31.1展示了这两种格式的位域划分。


S 指数 E (8位) 尾数 M (23位)

单精度浮点格式(binary32,32位)


S 指数 E (11位) 尾数 M (52位)

双精度浮点格式(binary64,64位)
IEEE 754单精度和双精度浮点格式的位域划分

表 31.2总结了各精度格式的关键参数。

格式总位宽指数位尾数位偏移量十进制精度
半精度 (binary16)1651015\sim3.3位
单精度 (binary32)32823127\sim7.2位
双精度 (binary64)6411521023\sim15.9位
四精度 (binary128)1281511216383\sim34.0位

IEEE 754常用浮点格式的参数

指数字段使用偏移编码(biased representation):实际指数值 = EbiasE - \text{bias}。以单精度为例,指数字段为8位(EE的范围为0\sim255),偏移量为127,因此实际指数范围为126-126+127+127E=0E=0E=255E=255被保留用于表示特殊值)。这种偏移编码的优点是使得浮点数的大小比较可以像无符号整数一样通过简单的字典序比较来完成——符号位相同时,指数越大数值越大,指数相同时尾数越大数值越大。

为什么偏移量是127而非128?

一个自然的问题是:单精度的8位指数字段有256个编码值(02550 \sim 255),去掉保留的全0(E=0E=0,用于零和非规格化数)和全1(E=255E=255,用于无穷和NaN),剩余254个有效编码值(12541 \sim 254),对应实际指数126-126+127+127。偏移量bias=127\text{bias} = 127使得这254个有效编码值的中点恰好位于E=127E = 127处,对应实际指数00

为什么不选bias=128\text{bias} = 128使得正负指数对称?如果bias=128\text{bias} = 128,则有效指数范围变为1128=1271 - 128 = -127254128=+126254 - 128 = +126。此时负指数范围比正指数范围多1——最大指数为+126+126而非+127+127。这意味着最大可表示值变小了(从3.4×1038\approx 3.4 \times 10^{38}下降到1.7×1038\approx 1.7 \times 10^{38}),而最小正规格化数却改善不大(从1.18×1038\approx 1.18 \times 10^{-38}5.88×1039\approx 5.88 \times 10^{-39})。

IEEE 754选择bias=2k11\text{bias} = 2^{k-1} - 1(而非2k12^{k-1})的根本原因在于最大化动态范围的上界。在科学计算和工程应用中,大数值(如天文距离、粒子物理能量)的出现频率不低于极小数值。偏移量为2k112^{k-1} - 1使得正指数比负指数多一个值(+127+127对比126-126),从而使最大值多了一倍。此外,这个选择使得倒数关系更加平衡:最大规格化数(2128\approx 2^{128})的倒数(2128\approx 2^{-128})仍然在非规格化数的表示范围内,不会直接下溢到零。

从硬件设计的角度来看,bias=2k11=0111k1\text{bias} = 2^{k-1} - 1 = 0\underbrace{11\ldots1}_{k-1}的二进制表示是除最高位外全1,这使得偏移量减法可以利用简单的位操作技巧来优化。

设计提示

偏移量的选择还影响了浮点比较器的设计。由于bias=127\text{bias} = 127,实际指数为0的数编码为E=127=011111112E = 127 = 01111111_2,正好是编码值范围的中点。对于两个同号正数的比较,可以直接将指数和尾数字段拼接后作为无符号整数比较——更大的编码值对应更大的浮点数。如果偏移量选为128,这个性质同样成立,但bias=127\text{bias} = 127的选择在实际指数=0= 0时的编码值(E=127E=127)恰好使得指数字段的最高位为0——这对某些分支预测和指数通路的优化有微弱的面积收益。

为什么使用隐含位"1"?

尾数字段存储的是小数部分MM,隐含的整数位"1"不需要存储。这意味着一个23位尾数字段实际提供了24位的精度(1+231 + 23位)。

隐含位设计的核心动机是免费获得一位额外精度。对于规格化浮点数,其有效数字的最高位必然为1(否则可以通过左移尾数并减小指数来使最高位变为1——这正是规格化的定义)。既然这一位始终为1,存储它就是浪费。通过隐含这个"1",23位的尾数字段提供了24位的精度——相当于在不增加存储成本的情况下多了1位有效数字,精度提升约213dB2^{-1} \approx 3\text{dB}

这个设计的代价是:(1)零不能用规格化形式表示(因为隐含位为1意味着最小的规格化值也大于零),必须用特殊编码(指数和尾数全零);(2)非规格化数的隐含位为0,打破了统一性,增加了硬件复杂度;(3)浮点运算单元在内部需要将隐含位显式地拼接到尾数前面,形成完整的有效数字进行计算。这个解包(unpack)操作是浮点流水线第一级的必要步骤:

verilog
// 解包: 从IEEE 754编码提取三元组 (sign, exp, mantissa_with_hidden_bit)
logic        sign;
logic [7:0]  biased_exp;
logic [23:0] full_mantissa;  // 24位: 隐含位 + 23位尾数

assign sign       = fp_in[31];
assign biased_exp = fp_in[30:23];

// 隐含位: 指数非零时为1(规格化数), 指数为零时为0(非规格化数/零)
assign full_mantissa = {|biased_exp, fp_in[22:0]};
// |biased_exp 是指数字段的OR归约: 指数非全零时为1

值得注意的是,并非所有浮点格式都使用隐含位。Intel的x87扩展精度格式(80位)使用显式整数位(explicit integer bit)——尾数字段为64位,其中最高位是显式的整数位。这种设计简化了硬件(不需要拼接隐含位),但浪费了一位存储空间。x87格式的这一特殊性还导致了一个有趣的问题:存在一类"非规范"(unnormal)编码——指数非零但整数位为0的值——在IEEE 754标准中没有定义,x87硬件需要额外的逻辑来处理这些非规范值。

浮点数轴上的分布

理解浮点编码设计的另一个角度是考察浮点数在数轴上的分布密度。在指数为ee的区间[2e,2e+1)[2^e, 2^{e+1})内,相邻两个可表示的浮点数之间的间距为:

ULP(e)=2ep+1 \text{ULP}(e) = 2^{e - p + 1}

其中pp为有效位数(单精度p=24p=24,双精度p=53p=53)。ULP(Unit in the Last Place)随指数ee指数增长——在[1,2)[1, 2)区间内,相邻浮点数的间距为2231.19×1072^{-23} \approx 1.19 \times 10^{-7}(单精度);而在[223,224)[2^{23}, 2^{24})区间内,间距增大到11——此时单精度浮点数只能精确表示整数值。

这种"离零越远,间距越大"的非均匀分布正是浮点数的核心特征。它意味着:

  • 浮点运算的相对精度是近似恒定的(2p\approx 2^{-p}),而绝对精度随数值大小变化。

  • 接近零的区域有最密集的可表示值,而远离零的区域(大数和极小的正数)间距极大。

  • 非规格化数的存在进一步填充了零附近的区域,使得在[0,21bias)[0, 2^{1-\text{bias}})范围内也有均匀间距的可表示值(间距为21biasp+12^{1-\text{bias}-p+1}),实现渐进下溢。

从硬件角度来看,ULP的概念直接关系到舍入逻辑的设计——"向上舍入1个ULP"在硬件中对应的操作是对尾数的最低有效位执行"+1"增量,这是一个简单的增量器操作。但要注意,增量可能导致进位传播——当尾数为全1时,增量溢出到指数字段,需要将指数加1并将尾数清零。这种"尾数溢出"的处理在舍入硬件中必须考虑。

为了更直观地理解浮点编码,考虑以下具体示例。将十进制数6.625-6.625表示为IEEE 754单精度浮点数:

  1. 符号位S=1S = 1(负数)。

  2. 将绝对值6.625转换为二进制:6.62510=110.1012=1.10101×226.625_{10} = 110.101_2 = 1.10101 \times 2^2

  3. 指数E=2+127=129=100000012E = 2 + 127 = 129 = 10000001_2

  4. 尾数M=101010000000000000000002M = 10101000000000000000000_2(隐含位"1"不存储)。

  5. 最终编码:1,即0xC0D40000

这个编码/解码过程在硬件中的实现虽然简单,但必须在每次浮点操作的输入端和输出端完成——输入端将IEEE 754编码"解包"为符号、指数和带隐含位的尾数三元组,输出端将运算结果"打包"回IEEE 754编码格式。解包和打包逻辑虽然简单,但占据了浮点流水线的第一级和最后一级的一部分关键路径。

浮点数的动态范围与精度

IEEE 754格式的设计在动态范围和精度之间做了精心的平衡。单精度格式的动态范围为103810^{-38}103810^{38}(约76个数量级),而精度约为7位十进制有效数字。双精度格式的动态范围扩展到1030810^{-308}1030810^{308}(约616个数量级),精度约为16位十进制有效数字。

从硬件角度来看,这种"指数+尾数"的分离编码方式意味着浮点数的运算需要对指数和尾数分别进行不同的处理——加法需要对阶(指数操作)后再做尾数加减(尾数操作),乘法需要指数相加(指数操作)且尾数相乘(尾数操作)。这种"分离处理、合并结果"的模式贯穿了所有浮点运算单元的设计。

在现代处理器中,浮点寄存器文件存储的数据宽度通常远大于单个标量浮点数的宽度。例如x86的XMM寄存器为128位,YMM寄存器为256位,ZMM寄存器为512位——它们可以存储多个单精度或双精度浮点数用于SIMD并行运算。RISC-V的向量寄存器(V扩展)宽度由VLEN参数决定,可以是128位到65536位。浮点功能单元的数据通路宽度需要与最宽的寄存器格式匹配,或通过拆分为多个较窄的操作来处理。

设计提示

在硬件实现中,浮点数据通路的宽度通常大于格式定义的位宽。例如单精度的内部尾数通路宽度通常为26\sim28位(24位有效位 + 保护位guard + 舍入位round + 粘滞位sticky),双精度则为55\sim57位。这些额外的精度位(guard/round/sticky,简称G/R/S位)是实现正确舍入所必需的,将在31.2.4 节中详细讨论。

舍入模式

浮点运算的一个本质特征是:精确结果通常需要比目标格式更多的位数来表示,因此必须进行舍入(rounding)。IEEE 754标准定义了四种舍入模式:

  1. 就近舍入到偶数(Round to Nearest, ties to Even,简称RNE):将结果舍入到最接近的可表示值。当结果恰好位于两个可表示值的中间时,选择最低有效位为0(偶数)的那个。RNE是默认的舍入模式,因为它在统计上是无偏的——中间值有一半的概率向上舍入,一半向下舍入。

  2. 向零舍入(Round toward Zero,简称RZ):将结果截断到绝对值不超过精确结果的最接近可表示值。对于正数,RZ等同于向下取整;对于负数,RZ等同于向上取整。RZ在整数与浮点转换中经常使用。

  3. 向正无穷舍入(Round toward Positive Infinity,简称RP):将结果舍入到不小于精确结果的最小可表示值。

  4. 向负无穷舍入(Round toward Negative Infinity,简称RN):将结果舍入到不大于精确结果的最大可表示值。

RP和RN模式在区间算术(interval arithmetic)中有重要应用:通过分别使用RP和RN对同一计算进行两次求值,可以获得结果的上界和下界,从而确定计算误差的范围。

舍入模式的硬件控制

舍入模式通过浮点控制/状态寄存器(如x86的MXCSR寄存器、RISC-V的fcsr寄存器)进行设置。在大多数程序中,舍入模式始终保持为默认的RNE模式。从硬件设计的角度来看,四种舍入模式的支持增加了舍入逻辑的复杂度——需要一个4选1的舍入决策多路选择器,根据舍入模式位和精度位来决定最终的舍入动作。然而,由于舍入逻辑位于浮点流水线的最后阶段,且舍入模式通常不频繁切换,这部分额外的硬件开销是可以接受的。

RISC-V的浮点ISA在舍入模式处理上有一个值得注意的设计选择:每条浮点运算指令的编码中包含一个3位的指令级舍入模式字段(rm字段),允许每条指令独立指定舍入模式而不需要修改全局的fcsr寄存器。当rm = 111(动态模式)时,使用fcsr中的全局舍入模式;否则使用指令编码中指定的静态舍入模式。

这种指令级舍入模式对硬件的影响是:舍入模式信号不再是一个全局的静态配置,而是需要随每条指令一起在流水线中传播。具体而言,rm字段在译码时被提取并写入发射队列的entry中,随指令一起流经流水线,最终在舍入级作为舍入决策的输入。这增加了流水线寄存器的宽度(额外的2\sim3位/指令),但面积开销极小。

指令级舍入模式的优势在于:不需要在两次舍入模式不同的运算之间插入fcsr写指令(这会引入序列化点),直接在指令编码中指定即可。这对区间算术和某些数值库的实现非常有利——同一段代码中可以交替使用RP和RN模式而不需要修改控制寄存器。x86的浮点指令(SSE/AVX)不提供指令级舍入模式(但AVX-512的EVEX编码增加了SAE——Suppress All Exceptions选项,可以指定RNE模式),ARM的A64浮点指令也不支持指令级舍入。

在硬件实现中,舍入决策依赖于三个额外的精度位:保护位(Guard bit, G)、舍入位(Round bit, R)和粘滞位(Sticky bit, S)。G位是精确结果中紧跟在最低有效位之后的第一位,R位是第二位,S位是所有更低位的OR归约——只要有任何更低位为1,S就为1。

为什么需要这三个位而非更多或更少?答案在于IEEE 754的舍入规则的信息需求。为了执行正确的舍入,硬件需要知道精确结果相对于两个相邻可表示值的位置关系——精确结果是更接近较小值、更接近较大值、还是恰好在两者中间。这个位置关系完全由G/R/S三个位决定:

  • G=0G = 0:精确结果在"截断值"和"截断值 + 0.5 ULP"之间的下半区,更接近截断值。

  • G=1,R=0,S=0G = 1, R = 0, S = 0:精确结果恰好等于"截断值 + 0.5 ULP",是正中间值(tie)。

  • G=1G = 1RS=1R \vee S = 1:精确结果在"截断值 + 0.5 ULP"之上的上半区,更接近增量值(截断值 + 1 ULP)。

只有这三种位置关系,因此三个精度位就足够完全确定舍入方向。如果只用G和R两个位(不用S),则无法区分"精确结果恰好在中间"和"精确结果略大于中间"——前者在RNE模式下需要向偶数舍入,后者需要向上舍入。S位正是用来做这个区分的。如果增加更多的精度位(如第三位、第四位),它们的信息可以被OR归约到S位中——不需要单独保留。因此G/R/S三个位是信息论意义上的最小充分集

从硬件实现的角度来看,G和R位是精确结果中特定位置的单个位,可以直接从移位器或乘法器的输出中提取。S位则是一个OR归约操作的结果——需要将所有比R位更低的位OR在一起。这个OR归约是一个树形结构,延迟为O(logn)O(\log n),但如前所述可以与桶形移位器集成,不增加额外的关键路径延迟。

表 31.3给出了RNE模式下的舍入决策逻辑。

GRS舍入动作说明
0×\times×\times截断精确结果更接近较小值
100向偶数舍入恰好在中间,看L位
101进位(+1)精确结果更接近较大值
11×\times进位(+1)精确结果更接近较大值

RNE舍入模式下的舍入决策(L=最低有效位)

特殊值

IEEE 754通过指数字段的两个极端值(全0和全1)来编码特殊值,这些特殊值的存在对浮点硬件的设计产生了深远的影响。

正负无穷(±\pm\infty

当指数字段为全1且尾数字段为全0时,表示正无穷或负无穷(由符号位决定)。无穷的产生有两种途径:(1)溢出——当运算结果超出可表示的最大值时;(2)非零数除以零。无穷参与后续运算时遵循数学上的直觉:+x=\infty + x = \infty×x=±\infty \times x = \pm\inftyx0x \neq 0)等。

NaN(Not a Number)

当指数字段为全1且尾数字段不为全0时,表示NaN。NaN分为两类:

  • 信号NaN(Signaling NaN, sNaN):尾数最高位为0。sNaN的存在意味着"无效操作"——当sNaN参与几乎任何运算时,都会触发Invalid Operation异常。sNaN主要用于调试:将未初始化的浮点变量设置为sNaN,任何对其的使用都会立即触发异常。

  • 安静NaN(Quiet NaN, qNaN):尾数最高位为1。qNaN在运算中会悄无声息地传播——任何与qNaN的运算都产生qNaN,不触发异常。这允许程序在不使用异常处理的情况下检测无效结果。

NaN是唯一不等于自身的值(NaNNaN\text{NaN} \neq \text{NaN}),这一特性在硬件比较器中需要特殊处理。当两个操作数都是NaN时,比较结果不是"大于"也不是"小于",而是"无序"(unordered)——这要求浮点比较器的输出至少有四种状态(大于、等于、小于、无序),而非整数比较器的三种状态。

NaN还涉及NaN传播(NaN propagation)规则:当一个操作的输入包含qNaN时,结果应当是该qNaN(如果有多个NaN输入,IEEE 754推荐传播第一个NaN,但具体行为由实现定义)。NaN传播逻辑通常在浮点流水线的特殊值快速路径中实现——当输入检查阶段检测到NaN时,跳过正常的运算流程,直接将NaN传递到输出。

正负零(±0\pm 0

当指数字段和尾数字段都为全0时,表示零(符号位区分+0+00-0)。在比较运算中,+0=0+0 = -0,但在某些运算中它们的行为不同:1/(+0)=+1 / (+0) = +\infty1/(0)=1 / (-0) = -\infty。双符号零的存在增加了浮点比较器的设计复杂度。

非规格化数(Denormal / Subnormal)

当指数字段为全0但尾数字段不为全0时,表示非规格化数。非规格化数的隐含位为0(而非规格化数的1),实际值为:

v=(1)S×0.M×21biasv = (-1)^{S} \times 0.M \times 2^{1 - \text{bias}}

特殊值的运算规则

IEEE 754标准详细定义了特殊值参与运算时的结果。表表 31.4列出了常见运算中涉及特殊值的结果规则。

运算结果
+\infty + \infty++\infty
+()\infty + (-\infty)NaN (Invalid)
×0\infty \times 0NaN (Invalid)
×x\infty \times x (x0x \neq 0, xx \neq NaN)±\pm\infty
x/0x / 0 (x0x \neq 0, xx \neq NaN)±\pm\infty (DivByZero)
0/00 / 0NaN (Invalid)
x\sqrt{-x} (x>0x > 0)NaN (Invalid)
+\sqrt{+\infty}++\infty
0\sqrt{-0}0-0

部分涉及特殊值的运算结果

这些规则必须在浮点流水线的特殊值快速路径中精确实现。快速路径的逻辑虽然不复杂(主要是一系列条件判断和MUX),但需要覆盖所有特殊值组合——仅对加法器而言,就需要处理两个操作数分别为规格化数、非规格化数、零、无穷和NaN五种类型的5×5=255 \times 5 = 25种组合。遗漏任何一种组合都可能导致严重的正确性问题。浮点单元的形式化验证中,特殊值路径的覆盖是最重要的验证目标之一。

对于FMA(三个操作数),特殊值组合更为庞大:53=1255^3 = 125种输入分类组合。虽然其中许多组合可以归并处理(例如"任一操作数为NaN\to结果为NaN"覆盖了大量组合),但完整的分类表仍然有数十条规则。表表 31.5列出了FMA中一些容易被忽视的特殊值规则。

操作结果说明
0×+c0 \times \infty + cNaN (Invalid)0×0 \times \infty无意义
0×+qNaN0 \times \infty + \text{qNaN}qNaNNaN传播优先
a×b+a \times b + \infty (abab有限)++\infty\infty主导
×b+()\infty \times b + (-\infty) (b>0b > 0)NaN (Invalid)\infty - \infty无意义
×b+()\infty \times b + (-\infty) (b<0b < 0)-\infty+()=-\infty + (-\infty) = -\infty
a×b+0a \times b + 0 (a,ba, b规格化)a×ba \times bFMA退化为纯乘法

FMA中容易被忽视的特殊值运算规则

浮点比较操作

浮点数的比较操作在硬件上看似简单,但IEEE 754的特殊值使其变得微妙。浮点比较器的设计需要处理以下规则:

  • +0=0+0 = -0:正零和负零在比较时相等,但位模式不同。比较器需要检测零值(指数和尾数都为0)并特殊处理。

  • NaN不参与正常比较:NaNNaN\text{NaN} \neq \text{NaN},且NaN与任何值的比较结果都是"无序"(unordered)。比较器的输出需要增加"无序"状态。

  • 规格化数的快速比较:对于两个同号的规格化正数,由于IEEE 754的编码特性(指数使用偏移编码,尾数使用隐含位),可以直接将整个32位或64位的位模式视为无符号整数进行比较。这是IEEE 754编码设计的一个精巧之处,极大简化了硬件比较器的实现。

浮点比较器通常与浮点加法器共享部分硬件(如指数比较器),并且通常在浮点执行端口上执行,延迟为2\sim3个周期。

浮点比较器的硬件实现流程如下:

  1. NaN检测:检查两个操作数是否为NaN(指数全1且尾数非零)。如果任一操作数为NaN,比较结果为"无序"(unordered),同时如果是sNaN则触发Invalid Operation异常。NaN检测的延迟仅为1\sim2个门延迟(指数的AND归约 + 尾数的NOR归约),可以与后续步骤并行。

  2. 符号比较:如果两个操作数符号不同(且都不是零),则正数大于负数,比较立即完成。特殊情况:如果两个操作数都是零(检测到指数和尾数都为0),则无论符号如何都相等(+0=0+0 = -0)。

  3. 幅度比较:对于同号的两个操作数,比较指数和尾数的拼接值。由于IEEE 754的偏移指数编码,对于正数,指数-尾数拼接后的无符号整数越大表示浮点值越大。对于负数,关系相反(无符号整数越大表示浮点值越小,即绝对值越大)。因此幅度比较可以使用一个标准的无符号整数比较器(减法器),根据符号位翻转比较方向。

  4. 结果编码:浮点比较的输出是四种状态之一:大于(GT)、等于(EQ)、小于(LT)、无序(UO)。这四种状态用2位编码,其中UO状态使得浮点比较结果不能直接用于条件分支——需要先检查UO标志,再根据GT/EQ/LT做分支决策。

RISC-V的浮点比较指令(FEQ.SFLT.SFLE.S)将比较结果直接写入整数寄存器(0或1),这涉及浮点域到整数域的跨域结果传输。x86的浮点比较指令(COMISSUCOMISD)将结果写入EFLAGS寄存器,同样涉及跨域传输。这种跨域传输的延迟通常为1\sim2个额外周期。

设计提示

浮点比较操作中的FMIN/FMAX指令(取两个浮点数中较大/较小的一个)在向量化代码中非常常见——例如ReLU激活函数可以实现为FMAX(x,0)\text{FMAX}(x, 0)。FMIN/FMAX的硬件实现基本上就是一个浮点比较器加上一个2选1 MUX,但IEEE 754-2019对NaN输入的处理规则(如果一个输入是NaN,返回另一个非NaN输入)增加了控制逻辑的复杂度。RISC-V的FMIN.S/FMAX.S遵循IEEE 754-2019的语义,而较早的x86 MINPS/MAXPS则遵循较简单的"如果任一输入为NaN,返回第二个操作数"的规则——两者在NaN处理上不兼容,移植代码时需要注意。

IEEE 754-2019标准还引入了totalOrder谓词,提供了一种包含NaN和±0\pm 0在内的全序比较。在totalOrder下,NaN<<<0<+0<<+<+NaN-\text{NaN} < -\infty < \ldots < -0 < +0 < \ldots < +\infty < +\text{NaN}。这种全序比较在排序算法中很有用(可以避免NaN导致的排序异常),但在硬件中较少直接支持——通常通过软件组合多条比较指令来实现。

非规格化数的存在填补了零与最小规格化数之间的"间隙",实现了渐进下溢(gradual underflow)——当运算结果逐渐趋近于零时,精度逐渐降低(因为有效位在减少),而不是突然跳到零。

为了量化渐进下溢的意义,考虑没有非规格化数的假想格式。最小正规格化数为m=1.0×21261.18×1038m = 1.0 \times 2^{-126} \approx 1.18 \times 10^{-38}(单精度)。如果没有渐进下溢,任何幅度小于mm的运算结果都会被直接冲刷为零——从mm00之间的整个数轴上没有任何可表示的值。这会导致一个严重的数值问题:aba - b可能产生零,但aba \neq b。例如a=1.0001×2126a = 1.0001 \times 2^{-126}b=1.0000×2126b = 1.0000 \times 2^{-126},两者的差0.0001×2126=1.0×21300.0001 \times 2^{-126} = 1.0 \times 2^{-130}小于mm,在无渐进下溢的系统中会被冲刷为零——导致"ab=0a - b = 0aba \neq b"的逻辑矛盾。有了非规格化数,0.0001×21260.0001 \times 2^{-126}可以用非规格化编码精确表示(虽然只有4位有效精度),避免了这个矛盾。

Kahan(IEEE 754标准的主要设计者之一)将渐进下溢描述为"在不增加存储成本的情况下,显著提高了数值程序的鲁棒性"。然而,非规格化数的硬件处理代价极高:由于隐含位的变化,常规的对阶和规格化逻辑需要增加特殊路径。许多处理器选择将非规格化数的处理交给微码辅助或操作系统陷阱,这将在31.6.2 节中详细讨论。

表 31.6总结了IEEE 754特殊值的编码规则。

类别符号指数尾数表示的值
正零0全0全0+0+0
负零1全0全00-0
非规格化数SS全00\neq 0(1)S×0.M×21bias(-1)^S \times 0.M \times 2^{1-\text{bias}}
规格化数SS0,\neq 0, \neq全1任意(1)S×1.M×2Ebias(-1)^S \times 1.M \times 2^{E-\text{bias}}
正无穷0全1全0++\infty
负无穷1全1全0-\infty
sNaNSS全1最高位=0,0\neq 0信号NaN
qNaNSS全1最高位=1安静NaN

IEEE 754特殊值的编码规则

在硬件实现中,特殊值的检测需要在浮点流水线的第一级完成。这个检测逻辑包括:(1)指数字段的全0/全1检测(使用AND/NOR门的归约操作);(2)尾数字段的全0检测(类似的归约操作);(3)尾数最高位的检查(区分sNaN和qNaN)。检测结果作为控制信号伴随数据一起流过流水线,在每个阶段引导特殊值的快速路径处理。

特殊值快速路径的微架构实现

特殊值快速路径(special value fast path)是浮点流水线中与正常运算路径并行的一条简化路径。当输入分类检测到特殊值时,快速路径直接生成结果,跳过耗时的正常运算(对阶、尾数运算、规格化等)。

快速路径的具体实现取决于运算类型。以浮点加法为例,快速路径需要处理的输入组合和对应结果如下:

  • 任一输入为qNaN:结果为该qNaN(NaN传播)。如果两个输入都是qNaN,传播第一个操作数的qNaN(或由实现定义)。

  • 任一输入为sNaN:结果为对应的qNaN(sNaN的最高尾数位置1变为qNaN),同时触发Invalid Operation异常。

  • +\infty + \infty(同号):结果为\infty

  • +()\infty + (-\infty):结果为qNaN(Invalid Operation异常)。

  • +x\infty + xxx为有限数):结果为\infty

  • 0+00 + 0:结果为+0+0(在RN模式下为0-0,其他模式下为+0+0)。

  • 0+x0 + xx0x \neq 0:结果为xx(直接传递,不需要运算)。

快速路径的硬件实现是一个基于输入分类标志的多路选择器树——根据(a_class,b_class,op)(a\_\text{class}, b\_\text{class}, \text{op})的组合选择对应的结果。输入分类标志在流水线第一级生成(延迟\sim2个门延迟),快速路径结果在第一级或第二级就可以确定。但为了与正常路径保持流水线同步,快速路径的结果通常也经过所有流水级(在空闲的流水线寄存器中传递),最终在最后一级通过MUX与正常路径的结果合并输出。

快速路径的存在不仅保证了特殊值的正确处理,还可以提供性能优化:当检测到快速路径激活时,正常路径的硬件可以被时钟门控,减少无意义的信号翻转和功耗。这种"检测到特殊值后门控正常路径"的优化在功耗敏感的设计中是有价值的——尽管特殊值出现的概率很低,门控逻辑的面积开销也很小(主要是分类标志到时钟门控信号的几个门延迟)。

浮点分类指令

RISC-V定义了FCLASS.S/FCLASS.D指令,将浮点数的分类结果直接写入整数寄存器。FCLASS返回一个10位的独热编码(one-hot encoding),分别对应负无穷、负规格化数、负非规格化数、0-0+0+0、正非规格化数、正规格化数、正无穷、sNaN和qNaN十种分类。

FCLASS指令的硬件实现完全复用了特殊值检测电路——指数全0/全1检测、尾数全0检测和尾数最高位检测的组合就足以确定所有分类。FCLASS的延迟通常为1\sim2个周期(比较和多路选择器),不需要经过浮点运算流水线。这条指令对于软件实现的浮点异常处理和数值分析工具非常有用——可以在不触发异常的情况下检查浮点值的分类。

x86没有直接等价的FCLASS指令——需要通过多条比较指令和位操作来实现类似的分类功能,效率较低。ARM的A64指令集也不提供直接的分类指令。RISC-V的FCLASS指令是其浮点ISA设计中的一个精巧之处。

案例研究 1 — x87与SSE的浮点模型差异

Intel x87浮点单元使用80位扩展精度(64位尾数 + 15位指数 + 1位符号)作为内部工作格式,所有运算都以80位精度进行,只在存储到内存时才舍入到目标精度。这种"扩展精度"模型会导致双重舍入(double rounding)问题:先在80位精度舍入一次,存储时又在目标精度舍入一次,两次舍入的组合可能产生与直接在目标精度运算不同的结果。

SSE/AVX浮点单元则严格按照目标精度进行运算和舍入——单精度操作产生单精度结果,双精度操作产生双精度结果——完全消除了双重舍入问题。这也是现代编译器默认使用SSE而非x87进行浮点运算的主要原因之一。从微架构设计的角度看,SSE模型简化了浮点流水线的设计,因为数据通路的宽度是确定的,不需要x87那样的动态精度转换逻辑。

浮点加法器

浮点加法(和减法)是最基本的浮点运算,但其硬件实现远比整数加法复杂。对于整数加法,两个操作数可以直接送入加法器;而对于浮点加法,必须首先将两个操作数的指数对齐(对阶),才能对尾数进行有意义的加法操作。这个"对阶"步骤是浮点加法复杂度的根源——它引入了移位器、前导零检测器、粘滞位计算等整数加法中不需要的硬件结构,使得浮点加法器的面积通常是同位宽整数加法器的3\sim5倍。

从延迟角度来看,浮点加法的典型延迟(3\sim4个周期)也显著高于整数加法(1个周期)。这个延迟差异在程序性能中有直接影响:浮点加法依赖链(如归约求和sum=sum+a[i]\text{sum} = \text{sum} + a[i])的吞吐量瓶颈是加法延迟,而非执行端口带宽。这也是为什么高性能数值计算中广泛使用循环展开和多路累加来隐藏浮点加法延迟。

对阶与尾数运算与规格化

浮点加法的基本流程可以分为三个主要步骤:

第一步:对阶(Exponent Alignment)。比较两个操作数的指数,将指数较小的操作数的尾数向右移位,使两个操作数的指数相等。移位量等于两个指数的差值ΔE=EaEb\Delta E = |E_a - E_b|。右移时被移出的位不能简单丢弃——它们被用于计算粘滞位(sticky bit),以确保后续舍入的正确性。

例如计算1.0×25+1.1×221.0 \times 2^5 + 1.1 \times 2^2时,指数差ΔE=52=3\Delta E = 5 - 2 = 3,需要将第二个操作数的尾数右移3位:1.1×22=0.0011×251.1 \times 2^2 = 0.0011 \times 2^5。然后对齐后的两个尾数相加:1.0000+0.0011=1.00111.0000 + 0.0011 = 1.0011,结果为1.0011×251.0011 \times 2^5

第二步:尾数运算。在对阶完成后,对两个对齐后的尾数进行加法或减法(取决于操作类型和操作数的符号)。这一步使用的是普通的定点加法器,但数据通路宽度要包含guard/round/sticky位。需要注意的是,浮点加法和减法在硬件层面没有本质区别——浮点减法ABA - B可以视为A+(B)A + (-B),通过翻转BB的符号位来实现。因此硬件中只需要一个"浮点加减法器",根据操作类型和两个操作数的符号来决定尾数运算是加法还是减法。这个决策逻辑被称为有效运算判断(effective operation determination)。

第三步:规格化与舍入。尾数运算的结果可能不是规格化形式(即最高有效位可能不是1)。规格化分为两种情况:

  • 右移规格化:当尾数加法产生了进位溢出(如1.111+1.000=10.1111.111\ldots + 1.000\ldots = 10.111\ldots),需要将结果右移1位,指数加1。这种情况只需要固定的1位右移。

  • 左移规格化:当尾数减法导致了大规模消去(如1.00011.0000=0.00011.0001 - 1.0000 = 0.0001),需要将结果左移多位直到最高有效位为1,指数相应减小。左移量可以从1位到nn位不等,需要前导零检测器来确定。

规格化完成后,对结果进行舍入操作。如果舍入导致了新的进位溢出(如1.11111.1111\ldots被向上舍入到10.000010.0000\ldots),还需要再执行一次1位的右移——这称为后规格化(post-normalization)。后规格化使得舍入与规格化之间形成潜在的迭代关系,增加了硬件设计的复杂度。

性能分析 1 — 浮点加法的延迟分析

一个朴素的单路径浮点加法器的关键路径包括:

  • 指数比较:一个nn位减法器,延迟约为 O(logn)O(\log n)

  • 尾数对阶移位:一个最多53位(双精度)的桶形移位器,延迟约为 O(logn)O(\log n)

  • 尾数加法:一个约55位的加法器,延迟约为 O(logn)O(\log n)

  • 前导零检测 + 规格化移位:延迟约为 O(logn)O(\log n)

  • 舍入加法:一个约53位的增量器,延迟约为 O(logn)O(\log n)

以上各步骤串行执行时,总延迟约为 5×O(logn)5 \times O(\log n)。对于双精度(n55n \approx 55),在先进工艺下这大约相当于3\sim4个时钟周期。为了将浮点加法流水线化到单周期吞吐量(每周期启动一次新的加法),需要将上述步骤分配到多个流水级中——典型的浮点加法器使用3\sim5级流水线。

完整的浮点加法流水线推导示例

为了将以上三步流程具体化,我们通过一个详细的数值示例来推导整个浮点加法的硬件执行过程。

问题:计算双精度浮点数A=1.5=1.12×20A = 1.5 = 1.1_2 \times 2^0B=0.0625=1.02×24B = 0.0625 = 1.0_2 \times 2^{-4}的和。

第一步(流水线第1级):解包与指数比较

  • AA的解包:SA=0S_A = 0EA=0+1023=1023E_A = 0 + 1023 = 1023(偏移编码),MA=1.10000M_A = 1.1000\ldots0(53位含隐含位)。

  • BB的解包:SB=0S_B = 0EB=4+1023=1019E_B = -4 + 1023 = 1019MB=1.00000M_B = 1.0000\ldots0

  • 指数差:ΔE=EAEB=10231019=4\Delta E = E_A - E_B = 1023 - 1019 = 4

  • 路径选择:ΔE=4>1|\Delta E| = 4 > 1,选择Far Path。

  • 有效运算:两个操作数同号且操作为加法\to有效运算为加法。

第二步(流水线第2级):对阶与尾数运算

  • 对阶移位:将BB的尾数右移ΔE=4\Delta E = 4位:

    MB=0.000100000M_B' = 0.000\underline{1}0000\ldots0,被移出的4位为10001000,粘滞位S=1000=1S = 1 \vee 0 \vee 0 \vee 0 = 1

    但实际上被移出的位是原始尾数的低4位:00010001(只有最低位为1——不对,让我们重新计算)。

    原始MB=1.00000M_B = 1.0000\ldots0,右移4位:0.0001  000000.0001\;0000\ldots0。被移出的位:第1\sim4位为00010001,其中第4位的"1"是隐含位移出的。

    更精确地:MB=1.0000052M_B = 1.\underbrace{0000\ldots0}_{52\text{位}}右移4位变为0.000100000480.0001\underbrace{0000\ldots0}_{48\text{位}},扩展为55位(含G/R/S位)后:MB=0.0001  00000    G=0,R=0,S=0M_B' = 0.0001\;0000\ldots0\;|\;G=0, R=0, S=0

  • 尾数加法:MA+MB=1.10000+0.0001  00000=1.1001  00000M_A + M_B' = 1.1000\ldots0 + 0.0001\;0000\ldots0 = 1.1001\;0000\ldots0

第三步(流水线第3级):规格化与舍入

  • 规格化:结果1.1001  000001.1001\;0000\ldots0已经是规格化形式(最高位为1),不需要移位。

  • 舍入:G=0,R=0,S=0G=0, R=0, S=0,精确结果没有被截断的位,不需要舍入。结果是精确的。

  • 指数:Eresult=EA=1023E_{\text{result}} = E_A = 1023(偏移编码),实际指数为00

  • 最终结果:1.10012×20=1.5625101.1001_2 \times 2^0 = 1.5625_{10},即1.5+0.0625=1.56251.5 + 0.0625 = 1.5625。正确。

这个简单的例子展示了Far Path的正常流程。如果我们改为计算1.51.4375=0.06251.5 - 1.4375 = 0.0625ΔE=0|\Delta E| = 0,走Near Path),则需要LZA检测前导零并执行大规模左移——过程更为复杂,但原理与前述一致。

上述三步流程在硬件中需要仔细处理许多细节。例如,在指数比较的同时还需要确定哪个操作数的指数更大——因为需要将指数较小的那个操作数进行右移。这涉及到一个"交换"操作:如果Eb>EaE_b > E_a,则需要交换AABB的角色。交换逻辑可以通过2选1 MUX来实现,但增加了数据通路的延迟。一种优化是在指数比较的同时预计算两种交换情况的结果,然后根据比较结果选择正确的一个——这是一种典型的推测并行(speculative parallelism)技巧。

另一个重要的实现细节是补码处理。当有效运算为减法时(两个操作数的有效符号不同),需要计算补码来执行减法。补码计算(取反加一)本身需要加法器,这与尾数加法器形成了潜在的资源竞争。双路径设计中,Near Path通常配置独立的补码/减法逻辑,Far Path的减法则集成在尾数加法器中。

减法的补码处理有两种实现策略:

策略一:二的补码减法。 计算ABA - B时,将BB取反后与AA相加,并在最低位注入进位Cin=1C_{\text{in}} = 1(实现A+B+1=ABA + \overline{B} + 1 = A - B)。这种方式只需要一个加法器,通过控制BB的取反和进位输入来切换加法/减法模式。这是最常用的实现方式,因为它不需要额外的硬件——加法器本身就可以通过XOR门来实现BB的条件取反。

策略二:独立减法器。 使用一个专门的减法器来计算ABA - B(或BAB - A),避免补码的"+1"操作。这种方式需要额外的硬件,但在Near Path中可能有优势——因为Near Path的操作数幅度接近,ABA - B的结果可能为负(B>AB > A),此时需要对结果取补码。如果使用策略一,可能需要两次补码操作(一次对BB取反,一次对负结果取反),而策略二可以通过检测结果符号来避免二次补码。

实际设计中,Near Path通常采用"条件交换 + 减法"的方式:先根据尾数大小比较来交换AABB(确保大数减小数),然后执行减法得到非负结果。这个尾数大小比较可以与指数比较并行执行,不增加关键路径。但对于ΔE=0|\Delta E| = 0的情况(两个操作数指数完全相同),尾数比较是必须的——它确定了最终结果的符号和哪个操作数应该被"减去"。

浮点加法器的操作数交换

在浮点加法的第一步(对阶)中,需要将指数较小的操作数的尾数向右移位。但在硬件接收到两个操作数时,并不知道哪个的指数更大——需要先比较指数,然后根据比较结果交换操作数。这个"比较-交换"操作在数据通路上引入了一个额外的2选1 MUX延迟。

一种常见的优化是推测交换:假设操作数AA的指数更大(EAEBE_A \geq E_B),在比较完成之前就开始按此假设执行后续操作。如果比较结果证实假设正确,后续操作继续;如果假设错误(EA<EBE_A < E_B),则从已经预计算的"交换版本"中选择正确的结果。这种方式需要同时计算两种交换情况的中间结果,但消除了比较-交换MUX在关键路径上的位置。

具体实现中,推测交换的面积开销取决于需要并行计算的逻辑深度。在双路径加法器中,指数比较通常在流水线的第一级完成,其结果在第二级之前可用,因此交换MUX可以放在第一级和第二级之间——不需要推测交换。推测交换主要用于单路径设计或需要极低延迟的场景。

然而,朴素的串行实现存在一个根本性的效率问题:对阶移位和规格化移位不可能同时发生大幅度的移动。这可以严格地证明。

大规模对阶与大规模消去的互斥性

定理:对于两个规格化浮点数的减法ABA - BA>B>0A > B > 0),如果对阶移位量ΔE=EAEB2\Delta E = E_A - E_B \geq 2,则尾数减法后的前导零数量最多为1。

证明:设AABB对阶后的有效数字分别为A=1.mAA' = 1.m_AB=0.000ΔE11.mBB' = 0.0\underbrace{0\ldots0}_{\Delta E - 1}1.m_B。由于ΔE2\Delta E \geq 2BB'至少右移了2位,因此B<0.012=0.25B' < 0.01_2 = 0.25。而A1.0A' \geq 1.0,因此AB1.00.25=0.75=0.112A' - B' \geq 1.0 - 0.25 = 0.75 = 0.11_20.1120.11_2的前导零为0,因此ABA' - B'的前导零最多为0位。如果考虑ABA' - B'可能略小于1(当A=1.0A' = 1.0B>0B' > 0),则AB(0.75,1.0)A' - B' \in (0.75, 1.0),前导零仍为0。唯一可能出现1位前导零的情况是进位借位的边界效应,但即使如此也不超过1位。\square

推论:大规模消去(k>1k > 1位的前导零)只可能在ΔE1|\Delta E| \leq 1时发生。

这个定理揭示了浮点加法中的一个根本性结构:对阶移位和规格化移位是互斥的——不可能同时需要大型右移位器(对阶)和大型左移位器(规格化)。在单路径设计中,这两个大型移位器串行在关键路径上,但实际上任何一次加法操作至多使用其中一个。这种"设备冲突但无实际并发"的模式正是双路径设计的出发点——将两种互斥的场景分配到两条独立的硬件路径上,使每条路径只需要一个大型移位器而非两个。

从时序的角度来看,单路径设计的关键路径长度为T=T对阶+T加法+TLZD+T规格化+T舍入T_{\text{单}} = T_{\text{对阶}} + T_{\text{加法}} + T_{\text{LZD}} + T_{\text{规格化}} + T_{\text{舍入}},而双路径设计的关键路径为T=max(TNear,TFar)T_{\text{双}} = \max(T_{\text{Near}}, T_{\text{Far}}),其中:

TNear=T加法+TLZA+T规格化TFar=T对阶+T加法+T舍入\begin{aligned} T_{\text{Near}} &= T_{\text{加法}} + T_{\text{LZA}} + T_{\text{规格化}} \\ T_{\text{Far}} &= T_{\text{对阶}} + T_{\text{加法}} + T_{\text{舍入}} \end{aligned}

注意Near Path不需要大型对阶移位器(ΔE1|\Delta E| \leq 1),Far Path不需要LZD和大型规格化移位器(前导零1\leq 1)。TT_{\text{双}}显著短于TT_{\text{单}},典型的改善约为30%\sim40%的关键路径延迟。

这个观察直接引出了双路径加法器的设计思想。

双路径加法器

双路径浮点加法器(Dual-Path Floating-Point Adder)是现代高性能处理器中最广泛采用的浮点加法器架构,最早由Oberman和Flynn(1997年)系统化提出。其核心思想是:将浮点加法分为两种情况,分别用不同的硬件路径处理,最后通过多路选择器选择正确的结果。

Near Path(近路径,ΔE1|\Delta E| \leq 1

当两个操作数的指数差不超过1时,走Near Path。此路径的特点是:

  • 对阶移位量最多为1位,无需大型移位器。

  • 尾数减法可能导致大规模消去(如1.00000.1111=0.000011.0000\ldots - 0.1111\ldots = 0.0000\ldots 1),需要前导零检测器(LZD)和大型左移位器进行规格化。

  • 由于大规模消去丢失了大量精度位,舍入的精度损失相对于消去本身可以忽略,舍入逻辑可以简化。

Far Path(远路径,ΔE>1|\Delta E| > 1

当两个操作数的指数差大于1时,走Far Path。此路径的特点是:

  • 对阶需要大型右移位器。

  • 尾数加减法的结果最多需要1位的右移(处理进位溢出)或1位的左移(处理减法的微小消去),无需大型左移位器。

  • 需要完整的舍入逻辑(包括guard/round/sticky位的处理)。

图 31.2展示了双路径浮点加法器的整体结构。

双路径浮点加法器架构。Near Path处理$|\Delta E| \leq 1$的情况(需要大规模规格化),Far Path处理$|\Delta E| > 1$的情况(需要大规模对阶移位)。
双路径浮点加法器架构。Near Path处理$|\Delta E| \leq 1$的情况(需要大规模规格化),Far Path处理$|\Delta E| > 1$的情况(需要大规模对阶移位)。

双路径设计的关键优势在于每条路径都消除了对方路径中最复杂的部件——Near Path不需要大型右移位器,Far Path不需要前导零检测器和大型左移位器。由于两条路径并行执行,总延迟由较慢的那条路径决定,而每条路径的延迟都比单路径设计短得多。

Far Path中的粘滞位计算

在Far Path中,对阶右移可能将大量的尾数位移出有效范围。这些被移出的位的信息必须通过粘滞位(sticky bit)来保留。粘滞位的定义是:被移出的所有位的OR归约——只要有任何被移出的位为1,粘滞位就为1。

粘滞位虽然只有1位,但它承载的信息量是关键的——它区分了"精确结果等于截断值"(S=0S=0,即被移出的位全为零)和"精确结果大于截断值"(S=1S=1,即被移出的位中至少有一个为1)。这个1位信息在RNE舍入模式的"ties to even"决策中是决定性的:当G=1,R=0G=1, R=0时,如果S=0S=0则舍入到偶数(可能不进位),如果S=1S=1则一定进位。因此粘滞位的正确计算对IEEE 754的正确舍入至关重要。

朴素的粘滞位计算需要一个nn位的OR归约树,与对阶移位串行。一种优化方法是将粘滞位的计算与移位器的层级结构集成:在桶形移位器的每一级,将"被移出"的位进行OR归约,各级的OR结果最终再OR在一起。这种集成式设计使粘滞位计算不增加额外的延迟。

具体实现中,桶形移位器(barrel shifter)由log2n\lceil \log_2 n \rceil级MUX组成,第kk级根据移位量的第kk位决定是否移位2k2^k位。在每一级中,"被移出"的位可以通过一个2k2^k位的OR门归约为1位粘滞标志sks_k。最终的粘滞位为所有级别的粘滞标志的OR:S=s0s1sK1S = s_0 \vee s_1 \vee \ldots \vee s_{K-1}

以8位右移4位为例(桶形移位器为3级):

  • 第0级(移1位):移位量第0位=0,不移位,s0=0s_0 = 0

  • 第1级(移2位):移位量第1位=0,不移位,s1=0s_1 = 0

  • 第2级(移4位):移位量第2位=1,移4位。被移出的低4位通过4输入OR门归约得到s2s_2

  • 最终:S=s0s1s2=s2S = s_0 \vee s_1 \vee s_2 = s_2

这种集成式粘滞位计算的关键优势是:OR归约逻辑与MUX移位逻辑完全并行——在每一级中,MUX选择移位或不移位的同时,OR门计算该级的粘滞标志。因此粘滞位计算不在桶形移位器的关键路径上,也不增加额外的流水级。

Far Path中的有效运算判断

Far Path还需要处理一个重要的控制逻辑问题:有效运算判断(effective operation determination)。浮点加法指令和浮点减法指令在硬件层面的区别仅在于是否翻转第二个操作数的符号位。有效运算(实际对尾数执行加法还是减法)取决于操作码和两个操作数的符号的组合:

  • FADD(A+BA + B):如果AABB同号,有效运算为加法;如果异号,有效运算为减法。

  • FSUB(ABA - B):如果AABB同号,有效运算为减法;如果异号,有效运算为加法。

在Far Path中,有效减法(effective subtraction)的结果只可能需要0或1位的规格化移位(因为ΔE>1|\Delta E| > 1保证了不会大规模消去)。但即使只是1位的规格化,也需要在Far Path中配置一个简单的条件左移逻辑。有效加法的结果则可能溢出1位,需要1位的右移——这与减法的左移方向相反。硬件中通常使用同一个1位双向移位器来处理这两种情况。

Near Path中的简化舍入

Near Path的一个重要特性是:当大规模消去发生时(左移量2\geq 2位),所有的guard/round/sticky位都变成了有效位的一部分,因此不需要舍入——结果已经是精确的(或仅需极简单的处理)。只有当左移量为0或1时才可能需要舍入,此时舍入逻辑可以大幅简化。这一观察使得Near Path的硬件实现显著轻于Far Path。

为什么大规模消去后不需要舍入?考虑两个pp位精度的操作数AABB执行减法ABA - BΔE1|\Delta E| \leq 1。精确的减法结果最多有p+1p + 1位有效位(含进位/借位)。如果结果有k2k \geq 2个前导零,则规格化左移kk位后,结果的有效位数为p+1kp1<pp + 1 - k \leq p - 1 < p——即结果的有效位数少于目标精度pp,没有需要舍掉的低位。结果是精确的。

仅当左移量为0或1时,结果的有效位数可能等于或超过pp,此时G/R/S位可能不全为零,需要舍入。但即使需要舍入,Near Path的舍入逻辑也比Far Path简单:

  • 左移0位时,只有G位可能为1(R和S位都在消去中被抵消了),舍入决策只需考虑G和L位。

  • 左移1位时,G位来自原始结果的最低有效位之后第一位,舍入同样简化。

这种简化使得Near Path的舍入逻辑面积约为Far Path的20%\sim30%,延迟也更短。

设计权衡 1 — 单路径vs双路径浮点加法器

  • 单路径设计:使用一条通用路径处理所有情况。硬件面积较小(只需一套移位器和加法器),但关键路径最长(大型右移位器 + 加法器 + LZD + 大型左移位器 + 舍入 串行在一起)。适用于面积和功耗敏感的低性能设计。

  • 双路径设计:面积增加约40%\sim60%(两条独立的数据通路),但关键路径显著缩短,可以实现更高的时钟频率或更少的流水线级数。现代高性能处理器几乎都采用双路径设计。

  • 三路径设计:学术研究中提出的进一步优化,将Far Path分为"远加法"和"远减法"两条路径。硬件面积进一步增加,但延迟收益有限,工业界较少采用。

以下代码展示了简化的双路径浮点加法器的SystemVerilog实现框架。代码仅展示近路径(Near Path)和远路径(Far Path)的选择逻辑和关键数据通路,省略了完整的舍入和异常处理。

verilog
module fp_add_dual_path #(
    parameter EXP_W  = 8,   // exponent width
    parameter MAN_W  = 23   // mantissa width (excl hidden bit)
)(
    input  logic [EXP_W+MAN_W:0] a, b,   // IEEE 754 packed
    input  logic                 sub,     // 1 = subtract
    output logic [EXP_W+MAN_W:0] result
);
    localparam FULL = MAN_W + 1;  // including hidden bit

    // Unpack
    logic             sa, sb;
    logic [EXP_W-1:0] ea, eb;
    logic [FULL-1:0]  ma, mb;
    assign {sa, ea, ma[MAN_W-1:0]} = a;
    assign {sb, eb, mb[MAN_W-1:0]} = b;
    assign ma[MAN_W] = |ea;  // hidden bit (1 for normalized)
    assign mb[MAN_W] = |eb;

    // Effective operation
    logic eff_sub;
    assign eff_sub = sa ^ sb ^ sub;

    // Exponent difference
    logic signed [EXP_W:0] exp_diff;
    assign exp_diff = $signed({1'b0,ea}) - $signed({1'b0,eb});

    // Path selection: Near if |dE|<=1 AND effective subtract
    logic use_near;
    assign use_near = eff_sub
                    & (exp_diff >= -1) & (exp_diff <= 1);

    // === Near Path (|dE| <= 1, eff subtract) ===
    logic [FULL+1:0] near_result;
    logic [EXP_W-1:0] near_exp;
    always_comb begin
        logic [FULL:0] m_big, m_small, diff;
        if (exp_diff >= 0) begin
            m_big   = {1'b0, ma};
            m_small = {1'b0, mb} >> exp_diff;
            near_exp = ea;
        end else begin
            m_big   = {1'b0, mb};
            m_small = {1'b0, ma};
            near_exp = eb;
        end
        diff = m_big - m_small;
        // LZA: count leading zeros (simplified)
        near_result = diff;
    end

    // === Far Path (|dE| > 1 or eff add) ===
    logic [FULL+1:0] far_result;
    logic [EXP_W-1:0] far_exp;
    always_comb begin
        logic [FULL+1:0] m_aligned;
        if (exp_diff >= 0) begin
            m_aligned = {1'b0, mb, 1'b0} >> exp_diff;
            far_exp   = ea;
        end else begin
            m_aligned = {1'b0, ma, 1'b0}>>(-exp_diff);
            far_exp   = eb;
        end
        if (eff_sub)
            far_result = {1'b0, ma, 1'b0} - m_aligned;
        else
            far_result = {1'b0, ma, 1'b0} + m_aligned;
    end

    // Path MUX
    logic [FULL+1:0] mux_result;
    logic [EXP_W-1:0] mux_exp;
    assign mux_result = use_near ? near_result : far_result;
    assign mux_exp    = use_near ? near_exp    : far_exp;

    // Simplified pack (no rounding / normalize)
    assign result = {sa, mux_exp, mux_result[MAN_W-1:0]};
endmodule

上述代码中,use_near信号是路径选择的关键:当且仅当有效运算为减法指数差ΔE1|\Delta E| \leq 1时,选择近路径。近路径中需要LZA(前导零预测)和大规模左移规格化,远路径中需要大规模右移对阶。两条路径并行执行,由最终的MUX选择结果。实际硬件中还需要完整的舍入逻辑(Guard/Round/Sticky位处理)和异常检测,此处为展示架构而省略。

前导零检测器

在Near Path中,尾数减法的结果可能有大量的前导零。例如两个非常接近的数相减(大规模消去),结果可能是0.0000001×0.00000\ldots 01\times的形式。为了将结果规格化(使最高有效位为1),需要知道前导零的数量,然后左移相应的位数。前导零检测器(Leading Zero Detector, LZD)或前导一检测器(Leading One Detector, LOD)正是用于完成这一任务的电路。

朴素的前导零检测方法是从最高位开始逐位扫描——这是一个O(n)O(n)的串行过程,对于53位的双精度尾数来说太慢了。高效的LZD使用类似于优先编码器的树形结构,将检测时间降低到O(logn)O(\log n)

一种经典的LZD设计基于递归分治。将nn位数分为高半部分和低半部分,分别检测各自的前导零数量。如果高半部分全为零,则前导零总数 = n/2n/2 + 低半部分的前导零数量(此处需注意高低半的含义);否则前导零总数 = 高半部分的前导零数量。这个递归过程形成一棵二叉树,深度为log2n\lceil \log_2 n \rceil,每一层的延迟为一个2选1MUX加一个OR门。

图 31.3展示了一个8位LZD的树形结构。

8位前导零检测器(LZD)的递归树形结构。每一层将检测范围翻倍,总深度为$\log_2 n = 3$层。
8位前导零检测器(LZD)的递归树形结构。每一层将检测范围翻倍,总深度为$\log_2 n = 3$层。

每一层合并节点的逻辑如下:如果高半部分的"全零"标志(VHV_H)为1(即高半全为零),则该层的前导零数量 = n/2n/2 + 低半部分的前导零数量;否则 = 高半部分的前导零数量。"全零"标志向上传播:当前层的全零标志 = VHVLV_H \wedge V_L(高半和低半都全零)。这种结构的延迟为每层一个MUX延迟加一个AND门延迟。

LZD的电路实现细节

树形LZD的每个叶节点处理2位输入,输出1位前导零计数和1位"全零"标志。以2位LZD为例,输入为(b1,b0)(b_1, b_0)b1b_1是高位),输出为前导零计数cc和全零标志vv

c=b1b0(高位为0且低位为1时,前导零 = 1)v=b1b0(两位都为0时全零)\begin{aligned} c &= \overline{b_1} \wedge b_0 \qquad \text{(高位为0且低位为1时,前导零 = 1)} \\ v &= \overline{b_1} \wedge \overline{b_0} \qquad \text{(两位都为0时全零)} \end{aligned}

b1=1b_1 = 1时,c=0c = 0(无前导零);当b1=0,b0=1b_1 = 0, b_0 = 1时,c=1c = 1(1个前导零);当b1=b0=0b_1 = b_0 = 0时,v=1v = 1(全零),cc的值无意义(由上层节点通过vv来处理)。

合并节点的逻辑更加关键。设高半部分的输出为(cH,vH)(c_H, v_H),低半部分为(cL,vL)(c_L, v_L),合并后的输出为(c,v)(c, v)

\begin{aligned} v &= v_H \wedge v_L \\ c &= v_H \;?\; \{1, c_L\} : \{0, c_H\} \end{aligned}$$ 式[式 (31.6)](/part5/ch31#eq:lzd-merge-c)中,如果高半全零($v_H = 1$),则前导零计数 = 半宽度(最高位为1)+ 低半的计数$c_L$;否则计数 = 0(最高位)+ 高半的计数$c_H$。这个MUX选择操作是每层的关键路径。 ::: info 硬件描述 2 — LZD合并节点的RTL实现 ```verilog module lzd_merge #(parameter HALF_W = 4, CNT_W = 2) ( input logic [CNT_W-1:0] cnt_hi, cnt_lo, input logic allzero_hi, allzero_lo, output logic [CNT_W:0] cnt_out, // 位宽 +1 output logic allzero_out ); assign allzero_out = allzero_hi & allzero_lo; // 若高半全零, 计数 = HALF_W + cnt_lo; 否则 = cnt_hi assign cnt_out = allzero_hi ? {1'b1, cnt_lo} // 最高位=1表示+HALF_W : {1'b0, cnt_hi}; endmodule ``` ::: 对于$n$位的LZD,树的深度为$\lceil \log_2 n \rceil$层,每层的延迟为1个MUX延迟 + 1个AND门延迟。以双精度55位尾数为例,LZD的深度为$\lceil \log_2 55 \rceil = 6$层,总延迟约为6个MUX延迟。在先进工艺下,一个MUX延迟约为1$\sim$2个FO4门延迟,因此55位LZD的总延迟约为6$\sim$12个FO4门延迟,大约相当于一个中等复杂度加法器的延迟。 #### LZA:前导零预测 {#sec:31-lza-detail} 一个关键的优化是<strong>前导零预测</strong>(Leading Zero Anticipation, LZA)。传统的LZD在尾数减法<strong>完成之后</strong>才能开始工作,这使得LZD和减法器形成串行依赖——减法器和LZD的延迟之和构成关键路径。LZA则是在减法<strong>进行的同时</strong>,直接根据两个操作数(减法之前)预测结果的前导零位置。 LZA的基本原理如下:对于两个操作数$A$和$B$的减法$A - B$,对每一位$i$定义三个位级信号: - $T_i = A_i \wedge \overline{B_i}$:<strong>传播位</strong>(Transmit),该位$A$为1、$B$为0,差的该位"倾向于"为1 - $Z_i = \overline{A_i} \wedge B_i$:<strong>生成位</strong>(Zero/Generate),该位$A$为0、$B$为1,差的该位"倾向于"需要借位 - $G_i = A_i \odot B_i = \overline{A_i \oplus B_i}$:<strong>相等位</strong>(Equal),两位相同,差的该位为0(不考虑借位) 前导零的位置大致对应于$T/Z$模式的转变点。更精确地说,Schmookler和Nowka(2001年)提出的经典LZA算法定义了一个位向量$f$: $$ f_i = \begin{cases} T_i \oplus Z_{i-1} & \text{if } G_i = 0 \\ G_i \oplus G_{i-1} & \text{if } G_i = 1 \end{cases}

实际上,最广泛使用的简化版LZA直接计算:

fi=TiZi1fi=(TiTi1Zi1)(ZiZi1Ti1) f_i = T_i \oplus Z_{i-1} \quad \text{或} \quad f_i = (T_i \wedge \overline{T_{i-1}} \wedge \overline{Z_{i-1}}) \vee (Z_i \wedge \overline{Z_{i-1}} \wedge \overline{T_{i-1}})

ff向量中第一个为1的位置(即ff的前导零数量)就是减法结果中前导零的预测值。这个预测可能有±1\pm 1位的误差,这是因为借位链的传播可能使实际前导零比预测值多或少1位。

ff向量的计算只涉及AABB的对应位之间的局部逻辑(不需要借位传播),因此可以在O(1)O(1)的门延迟内完成。随后对ff向量进行前导零计数(使用上述的树形LZD),延迟为O(logn)O(\log n)。关键在于:ff向量的计算和尾数减法可以完全并行执行,因为两者的输入相同(AABB),但计算路径完全独立。

考虑以下具体示例来理解LZA的工作原理:

AA:1.0 1 1 0 1 0 0 0
BB:1.0 1 1 0 0 1 1 1
T/Z/GT/Z/G:GGGGTZZZ
ABA - B:0.0 0 0 0 0 0 0 1

LZA通过分析T/Z/GT/Z/G序列来预测前导零位置。在此例中,从最高位开始都是GG(两位相同),直到出现第一个TTZZ的位置(第4位),LZA预测前导零约为4\sim5位。实际减法结果的前导零为6位(考虑借位传播),LZA的误差为+1+1位,通过后续的1位修正移位来补偿。

LZA误差的修正

LZA的±1\pm 1位误差需要在规格化移位之后进行修正。实际的修正策略有两种:

策略一:后置1位移位器。 先按LZA的预测值进行规格化左移,然后检查结果的最高位是否为1。如果不是(说明预测偏小1位),再做一次1位左移;如果最高两位都为1(说明预测偏大1位,极罕见情况),则做一次1位右移。这种方法的延迟为大型移位器 + 1个MUX延迟。

策略二:并行双版本。 同时计算两个版本的规格化结果——一个按LZA预测值移位,另一个按预测值+1+1移位——然后根据减法结果的最高位选择正确的版本。这种方法以额外的移位器面积换取消除了修正MUX在关键路径上的串行位置,但移位器面积近乎翻倍。

大多数工业设计采用策略一,因为1位修正移位器的延迟很短(约1个MUX延迟),而策略二的面积开销过大。

LZD与LZA的对比

表 31.7总结了LZD和LZA在硬件实现中的关键差异。

特性LZDLZA
输入减法结果(1个数)减法的两个操作数
时序必须在减法完成之后与减法并行执行
精度精确可能有±1\pm 1位误差
电路结构优先编码器树位级逻辑 + 优先编码器树
门延迟O(logn)O(\log n)O(1)+O(logn)O(1) + O(\log n)
面积较小较大(约多30%\sim50%)
应用位置低性能设计 / 独立使用双路径加法器的Near Path

前导零检测(LZD)与前导零预测(LZA)的对比

LZD/LZA的面积开销

对于双精度浮点加法器,LZA电路处理的位宽约为55位。使用树形结构的LZA占用的面积约为同位宽加法器的30%\sim40%,这是一个不可忽视的硬件开销。然而,与LZA并行化带来的性能收益相比(节省一个加法器延迟的关键路径),这个面积开销是完全值得的。

LZA的面积主要来自两部分:(1)ff向量的逐位计算逻辑——nn位的ff向量需要nn个小型组合逻辑单元,面积为O(n)O(n);(2)对ff向量的优先编码器树——与LZD的树结构相同,面积为O(n)O(n)。总面积约为2×O(n)2 \times O(n),而LZD仅需O(n)O(n)。额外的30%\sim50%面积来自ff向量计算逻辑。

设计提示

在双路径浮点加法器中,LZA与尾数减法器的并行化是Near Path性能的关键。在典型的3级流水线设计中,第一级完成指数比较和对阶;第二级同时启动尾数减法和LZA;第三级根据LZA结果执行规格化移位和舍入。如果LZA和减法器是串行的,则第二和第三级会合并成一个很长的级,要么增加周期时间,要么需要拆分为更多流水级。

值得注意的是,LZA只在Near Path(ΔE1|\Delta E| \leq 1)中使用。Far Path不需要LZA,因为Far Path的规格化移位最多为1位——这是双路径设计的又一个优势:不仅消除了每条路径中不需要的大型移位器,也消除了不需要的LZA电路。

舍入逻辑

舍入是浮点加法器中最容易出错也最难验证的部分。朴素的舍入实现是在规格化完成后,根据G/R/S位和舍入模式决定是否对结果执行"+1"操作(增量舍入)。然而,这个"+1"操作本身是一个加法,其延迟(约O(logn)O(\log n))被串行地加到了已经很长的关键路径上。

并行舍入

高性能浮点加法器通常采用并行舍入(parallel rounding)技术:在规格化移位的同时,同时计算"截断结果"和"截断结果+1"两个版本,然后根据舍入决策选择其中之一。这样,舍入的"+1"加法与规格化移位并行执行,不增加关键路径延迟。

具体实现中,增量器(+1电路)可以利用尾数加法器的空闲端口,或使用一个独立的快速增量器(利用进位链很短的特点——增量器的进位传播通常在低几位就终止了)。

为了更清晰地说明并行舍入的原理,考虑以下示例。假设规格化后的尾数(含G位)为:

Result = 1.01101001…01, G=1, R=0, S=1

并行计算两个版本:

  • 截断版本:1.01101001…0 (直接丢弃G/R/S位以下的部分)

  • 增量版本:1.01101001…1 (截断结果 + 1 ULP)

在RNE模式下,由于G=1G=1RS=1R|S=1(精确结果更接近较大值),舍入决策为"进位",选择增量版本作为最终结果。两个版本的计算与G/R/S位的舍入决策判断完全并行执行,最终通过一个2选1 MUX输出正确的版本。

舍入引起的指数调整(后规格化)

舍入操作的一个微妙后果是:向上舍入(+1 ULP)可能导致尾数溢出。例如,双精度的尾数为1.111111.1111\ldots1(52个1),加1后变为10.0000010.0000\ldots0——尾数超出了[1,2)[1, 2)的规格化范围。此时需要后规格化(post-normalization):将尾数右移1位(10.0001.00010.000\to 1.000),指数加1。

后规格化本身很简单(固定的1位右移和指数+1),但它引入了一个潜在的问题:指数加1可能导致指数溢出。如果舍入前的指数已经是最大值(E=254E = 254对于单精度),加1后E=255E = 255——这是±\pm\infty的编码。因此,后规格化逻辑需要检测指数溢出并在必要时将结果设置为±\pm\infty、产生Overflow异常标志。

这个"舍入\to后规格化\to可能的指数溢出\to可能的异常"链条说明了浮点运算中各步骤之间的级联效应——每个步骤的输出可能影响后续步骤的行为,增加了硬件设计的复杂性和验证的难度。在流水线实现中,后规格化通常与舍入在同一级完成,指数溢出检测则在下一级(或同一级的后半部分)完成。

舍入逻辑的硬件实现

将前述的所有舍入相关逻辑整合在一起,一个完整的浮点舍入单元需要包含以下组件:

  1. G/R/S位提取:从规格化后的尾数中提取Guard、Round和Sticky位。

  2. 舍入模式译码:根据CSR中的舍入模式位(2位编码4种模式)和结果符号位,生成舍入控制信号。

  3. 舍入决策逻辑:组合G/R/S/L位和舍入控制信号,产生"是否加1"的决策信号RupR_{\text{up}}

  4. 条件增量器:根据RupR_{\text{up}}信号,对尾数执行条件+1操作。

  5. 后规格化检测:检查增量后的尾数是否溢出,如果溢出则右移1位、指数加1。

  6. 指数溢出/下溢检测:检查最终指数是否在有效范围内。

  7. Inexact标志生成Inexact=GRS\text{Inexact} = G \vee R \vee S(只要有任何非零的被舍弃位,结果就是不精确的)。

表 31.8给出了四种舍入模式的完整决策逻辑。

舍入模式RupR_{\text{up}}说明
RNEG(RSL)G \wedge (R \vee S \vee L)G=1且(R或S为1,或L为1即中间值时选偶数)
RZ00永远截断(不加1)
RPsign(GRS)\overline{\text{sign}} \wedge (G \vee R \vee S)正数时有非零被舍弃位则加1
RNsign(GRS)\text{sign} \wedge (G \vee R \vee S)负数时有非零被舍弃位则加1

四种IEEE 754舍入模式的舍入决策逻辑

RupR_{\text{up}}信号的生成只需要几个逻辑门,延迟极短。条件增量器的延迟也通常很短——因为"+1"的进位传播在统计上很快终止(随机数据中,进位平均只传播\sim2位即终止)。但在最坏情况下(尾数为全1),进位需要传播整个尾数位宽,增量器的延迟与全尺寸加法器相当。

舍入预测

更激进的优化是舍入预测(rounding prediction):在尾数加法进行的同时,预测舍入结果。这需要提前知道G/R/S位的值,而G/R/S位取决于对阶移位——在Far Path中,对阶移位在尾数加法之前完成,因此G/R/S位在加法开始时已经可用。利用这一点,舍入决策可以在尾数加法完成的同一个周期内做出,然后在下一个周期执行条件加一。

性能分析 2 — 现代浮点加法器的典型延迟

以双精度双路径浮点加法器为例,在先进工艺(3\sim5 nm)下的典型流水线划分为:

流水级操作延迟
第1级指数比较、路径选择、对阶移位(Far)1个周期
第2级尾数加/减法 + LZA(Near)1个周期
第3级规格化移位 + 并行舍入 + 结果选择1个周期

总延迟为3个周期,吞吐量为1次加法/周期。这是现代高性能处理器中双精度浮点加法的典型配置。Intel Skylake到Golden Cove系列和AMD Zen系列的浮点加法延迟均为3\sim4个周期。

浮点乘法器

浮点乘法在概念上比浮点加法更简单:两个浮点数相乘时,指数相加、尾数相乘、符号异或——不存在对阶问题,也不存在大规模消去和规格化的问题。然而,尾数乘法的硬件实现是浮点乘法器中最主要的面积和延迟开销——一个53×\times53位的乘法器阵列(双精度)包含27个部分积(Booth编码后)和多层CSA压缩树,其晶体管数量和布线复杂度都非常可观。

指数加法

浮点乘法的指数处理非常直观。给定两个操作数A=(1)SA×1.MA×2EAbiasA = (-1)^{S_A} \times 1.M_A \times 2^{E_A - \text{bias}}B=(1)SB×1.MB×2EBbiasB = (-1)^{S_B} \times 1.M_B \times 2^{E_B - \text{bias}},乘积的指数为:

Eresult=(EAbias)+(EBbias)+bias=EA+EBbias E_{\text{result}} = (E_A - \text{bias}) + (E_B - \text{bias}) + \text{bias} = E_A + E_B - \text{bias}

在硬件实现中,这是一个简单的加法后减去偏移量的操作。对于单精度,指数为8位,偏移量为127;对于双精度,指数为11位,偏移量为1023。指数加法器的延迟远小于尾数乘法器,因此不在关键路径上。

需要注意的是,指数加法的结果可能溢出(超出最大指数,产生±\pm\infty)或下溢(低于最小指数,产生非规格化数或零)。溢出和下溢的检测可以通过扩展指数位宽(增加1\sim2位)来实现——在扩展位宽下进行加法,然后检查结果是否超出有效范围。具体而言,对于双精度,11位指数扩展为13位来进行加法,这样可以检测到超过+2046+2046(溢出)或低于00(下溢)的情况。扩展后的指数加法器可以共享前面讨论的整数ALU中的加法器设计技术(如Kogge-Stone或Han-Carlson),但由于指数位宽较窄(13位),即使使用简单的行波进位加法器也不会成为瓶颈。

符号位的计算最为简单:Sresult=SASBS_{\text{result}} = S_A \oplus S_B(异或)。

指数的特殊情况处理

指数加法阶段还需要处理多种特殊情况:

  • 任一操作数为±\pm\infty(指数全1,尾数全0):结果为±\pm\infty(除非另一个操作数为零,此时结果为NaN)。

  • 任一操作数为零(指数全0,尾数全0):结果为±0\pm 0

  • 任一操作数为NaN:结果传播NaN。

  • 任一操作数为非规格化数:需要修正指数(非规格化数的实际指数为1bias1 - \text{bias}而非0bias0 - \text{bias})。

这些特殊情况通常在流水线的第一级通过输入分类逻辑检测,检测到特殊情况后绕过正常的运算流程,直接在快速路径中产生结果。快速路径的存在意味着特殊值的处理延迟通常只有1\sim2个周期(远小于正常运算的4\sim5个周期),但更重要的是它确保了特殊值不会在正常数据通路中引起错误行为。

尾数乘法

尾数乘法是浮点乘法器的核心,其硬件实现可以直接复用第 30.0 章中介绍的整数乘法器技术。

对于单精度,两个24位(含隐含位)的尾数相乘,产生一个48位的乘积。对于双精度,两个53位的尾数相乘,产生一个106位的乘积。在实际的浮点乘法器中,通常使用以下技术组合来实现高效的尾数乘法:

  1. Booth编码:使用基-4 Booth编码将部分积的数量减半。对于53位双精度尾数,原本需要53个部分积,Booth编码后减少到27个。

  2. Wallace树/Dadda树压缩:使用CSA(Carry-Save Adder)树将27个部分积逐层压缩为2个(一个和向量和一个进位向量)。

  3. 最终加法:使用一个快速加法器(如Kogge-Stone加法器)将最终的两个向量相加,得到精确的乘积。

设计提示

在许多处理器设计中,浮点乘法器和整数乘法器共享同一个乘法器阵列。这是因为尾数乘法本质上就是无符号整数乘法——24位×\times24位(单精度)或53位×\times53位(双精度)的无符号乘法。通过在乘法器阵列的输入端添加多路选择器,可以在整数乘法和浮点尾数乘法之间复用同一套Booth编码器和Wallace树。这种共享显著节省了芯片面积——一个53位×\times53位的乘法器阵列(含Booth编码和Wallace树)在先进工艺下仍然是一个面积可观的结构。Apple和ARM的处理器设计中广泛采用了这种共享策略。

浮点乘法的规格化过程比加法简单得多。两个规格化数的乘积的形式为1.x×1.y=1.xy1.x \times 1.y = 1.xy\ldots1x.y1x.y\ldots(当乘积2\geq 2时),因此规格化最多只需要1位的右移。这意味着浮点乘法器不需要前导零检测器和大型左移位器——这是浮点乘法在某些方面比浮点加法更"简单"的原因。

乘积截断与舍入

虽然尾数乘法器产生的是完整的106位(双精度)乘积,但最终结果只需要53位有效位加上G/R/S位。因此,乘积的低位部分可以在早期阶段就进行截断和OR归约(计算粘滞位),而不需要将完整的106位乘积传递到流水线的最后阶段。这种早期截断技术可以减少后续加法器和规格化移位器的宽度,节省面积和功耗。

具体而言,乘积的高54位(含可能的1位规格化右移)被保留用于精确的结果计算,低52位通过OR树归约为粘滞位。这个OR树可以集成到Wallace树的最后几层中,利用CSA压缩过程中的中间结果来计算,进一步避免额外的延迟。

混合精度乘法

现代处理器的浮点乘法器通常需要同时支持半精度(16位)、单精度(32位)和双精度(64位)的乘法。一种常见的实现方式是设计一个双精度宽度的乘法器阵列(53位×\times53位),然后通过分区(partitioning)来支持较低精度的运算。例如,一个53×\times53的乘法器可以同时执行两个24×\times24的单精度乘法(利用乘法器阵列的上半部和下半部),或者四个11×\times11的半精度乘法。这种SIMD分区乘法器在向量浮点单元中广泛使用。

分区乘法器的实现需要在乘法器阵列中插入分区屏障(partition barrier),阻止进位跨越分区边界传播。具体来说,在Booth编码阶段,分区屏障确保高位分区的Booth编码不依赖低位分区的操作数位;在CSA压缩树中,屏障阻断跨分区的进位传播;在最终CPA中,屏障将进位链截断。

分区屏障的控制信号来自精度模式选择——一个2位的模式信号控制屏障的激活/去激活:

  • 双精度模式:所有屏障去激活,整个53×\times53阵列作为一个统一的乘法器工作。

  • 单精度模式:中间位置的屏障激活,将阵列分为两个24×\times24的乘法器。

  • 半精度模式:更多的屏障激活,将阵列分为四个11×\times11的乘法器。

分区屏障的硬件开销很小——每个屏障只需要一个AND门来门控进位信号。但分区的存在对布局布线有影响:分区边界处的布线密度会降低(因为跨分区的信号被截断),可能影响整体布线的均匀性。

硬件描述 3 — 分区乘法器中CSA的分区逻辑

verilog
// 全加器(FA)的分区版本
// partition_en=1时, 截断进位输出(分区模式)
// partition_en=0时, 正常传播进位(双精度模式)
module partitioned_fa (
  input  logic a, b, cin,
  input  logic partition_en,  // 分区使能
  output logic sum, cout
);
  assign sum  = a ^ b ^ cin;
  // 分区使能时, 截断进位输出
  assign cout = partition_en ? 1'b0 : ((a & b) | (a & cin) | (b & cin));
endmodule

在实际设计中,分区屏障不是简单地将进位设为零,而是通过更精细的逻辑来确保分区后各子乘法器的结果正确。这包括处理Booth编码中跨分区的符号扩展和部分积的对齐。

乘法器的流水线化

为了实现高吞吐量,浮点乘法器通常被流水线化为4\sim5级。表表 31.9展示了一个典型的双精度浮点乘法器流水线。

流水级名称操作
1Booth编码Booth编码生成部分积;指数加法;符号异或
2部分积压缩(上半)Wallace/Dadda树的前几层CSA压缩
3部分积压缩(下半)Wallace/Dadda树的后几层CSA压缩,生成最终的和/进位向量
4最终加法 + 规格化CPA加法器完成最终相加;0/1位规格化移位
5舍入根据G/R/S位和舍入模式完成舍入;异常检测

双精度浮点乘法器的典型流水线划分

在某些设计中,第4级和第5级可以合并——通过在最终加法的同时预计算舍入结果,将规格化、最终加法和舍入压缩到单个流水级中。Intel的浮点乘法器采用的就是这种优化,将总延迟控制在4个周期。

一个值得注意的优化是舍入信息的早期提取。在Wallace/Dadda树压缩的最后几层,虽然最终的精确和尚未计算出来(还需要CPA),但乘积的低位(用于计算G/R/S位)已经可以部分确定。具体来说,乘积的最低若干位(低于最终截断点的位)在CSA压缩过程中可以提前"收割"——将这些低位通过OR归约提前计算粘滞位,而不需要等待CPA完成。这种粘滞位早期提取(early sticky extraction)技术可以使舍入决策在CPA完成之前就部分完成,进一步压缩最终CPA+舍入的总延迟。

Booth编码对浮点乘法的影响

基-4 Booth编码(也称Booth-2编码)是浮点尾数乘法中最常用的部分积缩减方法。对于53位的双精度尾数,原始方法需要53个部分积,而Booth-2编码将部分积数量减少到53/2+1=27\lceil 53/2 \rceil + 1 = 27个。每个Booth编码的部分积是被乘数的00±1\pm 1±2\pm 2倍——±1\pm 1倍通过直接取值或取反实现,±2\pm 2倍通过左移1位实现。

Booth编码对浮点乘法有一个特殊的影响:由于浮点尾数的最高位始终为1(规格化数),最高几位的Booth编码模式是受限的。例如,两个规格化数的尾数1.MA1.M_A1.MB1.M_BMAM_AMBM_B的最高位为1(隐含位),这意味着最高位的Booth编码总是非零的——可以利用这一特性来简化最高位部分积的生成逻辑。这种隐含位感知的Booth优化在某些设计中可以节省几个门的延迟和面积。

处理器延迟(周期)吞吐量(周期/次)
Intel Skylake\simGolden Cove40.5
AMD Zen 2\simZen 53\sim40.5\sim1
Apple Firestorm/Everest40.5
ARM Cortex-X43\sim41

现代处理器浮点乘法延迟(双精度FMUL)

吞吐量为0.5个周期/次意味着处理器每周期可以启动两次浮点乘法——这通过配置两个独立的浮点乘法流水线来实现。在大多数现代处理器中,这两个乘法流水线实际上是FMA单元的一部分——浮点乘法被编码为a×b+0a \times b + 0的FMA操作,通过将加数cc设为零来执行。

性能分析 3 — 浮点乘法器的面积分布

对于一个双精度浮点乘法器(53×\times53位),各子模块的典型面积分布如下:

  • Booth编码器和部分积生成:\sim15%

  • Wallace/Dadda CSA压缩树:\sim40%

  • 最终CPA(Carry-Propagate Adder):\sim15%

  • 规格化、舍入和异常逻辑:\sim10%

  • 流水线寄存器:\sim20%

CSA压缩树是面积的主要贡献者,因为它包含了27个部分积的多层3:2和4:2压缩器,每层的宽度约为106位。流水线寄存器的面积也不可忽视——在4\sim5级流水线中,每级都需要存储约106位的中间和/进位向量,总的流水线寄存器位数约为106×2×4=848106 \times 2 \times 4 = 848位。

乘法器流水线的切割点选择

Wallace树/Dadda树的CSA层级结构为流水线切割提供了天然的分界点。一个27个部分积的Dadda树需要约5层CSA,每层的延迟约为一个全加器(FA)的延迟。流水线寄存器可以在CSA层级之间插入,将压缩过程分为多个流水级。

选择切割点时需要考虑以下因素:

  • 延迟均衡:每个流水级的延迟应当尽量接近,避免某一级成为瓶颈限制时钟频率。

  • 寄存器开销:CSA中间结果的位宽与部分积相当(约106位),流水线寄存器的面积和功耗不可忽视。在CSA的每一层之间都插入寄存器是不划算的——通常将2\sim3层CSA合并为一个流水级。

  • 最终加法器的匹配:最后一级CPA的延迟通常大于单层CSA的延迟,因此CPA单独占一个流水级是合理的。

设计提示

流水线切割点的选择也影响了旁路网络(bypass network)的设计。如果乘法器的结果需要被后续指令立即使用,旁路网络需要在流水线的最后一级(舍入后)提供结果——这意味着乘法结果的旁路延迟等于乘法器的总流水线延迟(4\sim5个周期)。对于乘法结果作为加法器输入的情况(如fadd f3, f1, f2依赖于前一条fmul f2, f4, f5),加法器需要等待乘法器的全部延迟才能获得操作数——这就是FMA指令的另一个优势所在,它将乘法和加法融合在一个流水线中,消除了乘法到加法的旁路延迟。

在面向SIMD向量运算的浮点乘法器中,流水线化还需要考虑操作数宽度的复用。例如一个支持256位AVX的FMA单元,在执行四个并行的双精度乘法时需要四个53×\times53的乘法器同时工作。这四个乘法器可以共享Booth编码的控制逻辑,但各自维护独立的Wallace树和最终加法器。在某些设计中,这四个乘法器被进一步组织为一个大的212×53212 \times 53的乘法器阵列(拼接四组部分积),通过统一的Wallace树进行压缩——这种方式减少了独立Wallace树之间的布线拥挤,但增加了整体的布局布线复杂度。

融合乘加单元FMA

融合乘加(Fused Multiply-Add, FMA)是现代浮点单元中最重要的运算单元,也是衡量处理器浮点性能的核心指标。FMA将乘法和加法融合为单一操作,在一次流水线通过中完成a×b+ca \times b + c的计算,不仅提高了运算吞吐量(一条指令完成两个浮点操作),更重要的是提供了比分开执行乘加更高的数值精度。FMA的概念最早由IBM在1990年代的POWER架构中引入,现已成为所有主流ISA的标准配置。

FMA的数学定义

FMA执行的运算为:

FMA(a,b,c)=a×b+c \text{FMA}(a, b, c) = a \times b + c

关键在于FMA只执行一次舍入——先计算a×ba \times b的完整精确乘积(不舍入),然后与cc精确相加,最后对完整的精确结果执行一次舍入。这与先独立执行乘法(舍入一次)再执行加法(舍入一次)的分开计算方式有本质区别。

单次舍入 vs 两次舍入的精度分析

FMA的"单次舍入"特性在IEEE 754-2008标准中被明确要求:FMA的结果必须等价于先计算a×ba \times b的精确无限精度乘积,然后与cc精确相加,最后对无限精度的结果执行一次目标精度的舍入。这个定义决定了FMA硬件的内部精度必须足以容纳完整的乘积——对于双精度,这意味着至少106位的内部尾数宽度。

为什么单次舍入比两次舍入更精确?根本原因在于舍入误差的累积。设fl(x)\text{fl}(x)表示将精确值xx舍入到浮点格式,ϵ\epsilon为机器精度(单精度ϵ=224\epsilon = 2^{-24},双精度ϵ=253\epsilon = 2^{-53})。对于分开执行的乘加:

fl(a×b)=a×b×(1+δ1),δ1ϵfl(fl(a×b)+c)=(fl(a×b)+c)(1+δ2),δ2ϵ\begin{aligned} \text{fl}(a \times b) &= a \times b \times (1 + \delta_1), \quad |\delta_1| \leq \epsilon \\ \text{fl}(\text{fl}(a \times b) + c) &= (\text{fl}(a \times b) + c)(1 + \delta_2), \quad |\delta_2| \leq \epsilon \end{aligned}

组合这两个误差,分开计算的总误差界为:

fl(fl(a×b)+c)(ab+c)abϵ+ab+cϵ+O(ϵ2)|\text{fl}(\text{fl}(a \times b) + c) - (ab + c)| \leq |ab| \cdot \epsilon + |ab + c| \cdot \epsilon + O(\epsilon^2)

而FMA的单次舍入误差界仅为:

FMA(a,b,c)(ab+c)ab+cϵ|\text{FMA}(a, b, c) - (ab + c)| \leq |ab + c| \cdot \epsilon

关键的差异在于abϵ|ab| \cdot \epsilon项——当abab+c|ab| \gg |ab + c|时(即发生灾难性消去),分开计算的误差可以远大于FMA。具体地,如果ab+c0ab + c \approx 0ab|ab|很大,分开计算的相对误差abϵ/ab+c|ab| \cdot \epsilon / |ab + c|可以趋向无穷大,而FMA的相对误差始终有界于ϵ\epsilon

性能分析 4 — FMA精度优势的量化示例:灾难性消去

考虑以下双精度计算,展示灾难性消去(catastrophic cancellation)的具体影响。

a=1+230a = 1 + 2^{-30}, b=1+230b = 1 + 2^{-30}, c=(1+229)c = -(1 + 2^{-29})

精确结果

a×b=(1+230)2=1+229+260a×b+c=(1+229+260)(1+229)=260\begin{aligned} a \times b &= (1 + 2^{-30})^2 = 1 + 2^{-29} + 2^{-60} \\ a \times b + c &= (1 + 2^{-29} + 2^{-60}) - (1 + 2^{-29}) = 2^{-60} \end{aligned}

分开计算(假设双精度,53位尾数):

  • fl(a×b)=fl(1+229+260)\text{fl}(a \times b) = \text{fl}(1 + 2^{-29} + 2^{-60})。由于2602^{-60}位于有效位的第60位,超出了53位精度,被舍入丢弃。因此fl(a×b)=1+229\text{fl}(a \times b) = 1 + 2^{-29}

  • fl(fl(a×b)+c)=fl((1+229)(1+229))=0\text{fl}(\text{fl}(a \times b) + c) = \text{fl}((1 + 2^{-29}) - (1 + 2^{-29})) = 0

分开计算的结果为00,而精确结果为2608.67×10192^{-60} \approx 8.67 \times 10^{-19}相对误差为100%——结果完全错误。

FMA计算

  • FMA保留a×ba \times b的精确乘积(106位),即1+229+2601 + 2^{-29} + 2^{-60},不进行中间舍入。

  • 精确乘积与c=(1+229)c = -(1 + 2^{-29})相加:1+229+2601229=2601 + 2^{-29} + 2^{-60} - 1 - 2^{-29} = 2^{-60}

  • 最终对2602^{-60}进行一次舍入——由于2602^{-60}可以精确表示为双精度浮点数,舍入误差为零

FMA的结果为精确值2602^{-60}

这个例子清晰地展示了分开计算的两次舍入如何在消去场景中导致精度的完全丧失。第一次舍入(乘法)丢弃了2602^{-60}这个在减法后成为唯一有效信息的位。FMA的单次舍入保留了完整的乘积,使得减法后仍能获得正确的结果。

在实际应用中,灾难性消去出现在许多关键算法中:

  • 二次方程ax2+bx+c=0ax^2 + bx + c = 0的求根公式中,当b24acb^2 \approx 4ac

  • 矩阵行列式的计算中,当矩阵接近奇异时

  • 有限差分法中相邻采样点函数值的相减

  • Newton-Raphson迭代中残差r=aq×br = a - q \times b的计算

FMA的存在使得这些算法的数值稳定性得到了硬件级的保障。

FMA的精度优势使其成为以下算法的核心:

  • 矩阵乘法(GEMM):大量的乘累加操作直接映射到FMA。在BLAS Level-3的DGEMM实现中,计算核心的99%以上的浮点操作都是FMA。

  • 多项式求值:Horner法 p(x)=((anx+an1)x+an2)x+p(x) = ((a_n x + a_{n-1})x + a_{n-2})x + \ldots 的每一步都是一个FMA。FMA的单次舍入保证了Horner求值链的误差不超过理论下界。

  • Newton-Raphson迭代:除法和平方根的软件实现依赖FMA的精度保证来实现正确的收敛。FMA的关键作用在于计算残差r=aq×br = a - q \times b时不引入额外的舍入误差。

  • 补偿求和(Compensated Summation):Kahan求和算法利用FMA的精确乘加来计算误差补偿项。

  • 双精度模拟(Double-Double Arithmetic):使用两个双精度值来模拟约30位十进制精度的运算,其核心算法(如Dekker乘法和Knuth二和)可以利用FMA来简化和加速。

性能分析 5 — 五步推导:FMA vs 分离乘加的精度差异量化

本推导为上述灾难性消去示例提供一个通用的误差界分析框架。

第1步:定义变量。aa, bb, cc为双精度浮点数,ϵ=253\epsilon = 2^{-53}(机器精度),精确结果r=ab+cr = ab + c。定义消去度κ=ab/ab+c\kappa = |ab| / |ab + c|——当ababcc近似抵消时,κ1\kappa \gg 1

第2步:分离乘加的误差界。fl(ab)=ab(1+δ1)\text{fl}(ab) = ab(1+\delta_1)fl(fl(ab)+c)=(fl(ab)+c)(1+δ2)\text{fl}(\text{fl}(ab)+c) = (\text{fl}(ab)+c)(1+\delta_2)δiϵ|\delta_i| \leq \epsilon。展开得绝对误差Esepabϵ+ab+cϵ+O(ϵ2)|E_{\text{sep}}| \leq |ab|\epsilon + |ab+c|\epsilon + O(\epsilon^2)。相对误差:

Eseprκϵ+ϵ(κ+1)ϵ\frac{|E_{\text{sep}}|}{|r|} \leq \kappa \cdot \epsilon + \epsilon \approx (\kappa + 1)\epsilon

第3步:FMA的误差界。FMA保留abab的精确乘积(106位),只在最终结果处执行一次舍入:FMA(a,b,c)rrϵ|\text{FMA}(a,b,c) - r| \leq |r| \cdot \epsilon。相对误差恒为ϵ\epsilonκ\kappa无关

第4步:精度差异倍数。分离乘加与FMA的相对误差之比为(κ+1)ϵ/ϵ=κ+1(\kappa+1)\epsilon / \epsilon = \kappa + 1。当κ=1015\kappa = 10^{15}(双精度可表示的典型大消去度)时,分离乘加的相对误差可以是FMA的101510^{15}倍——即15个十进制有效位全部丢失

第5步:代入具体数值验证。a=b=1+230a = b = 1 + 2^{-30}c=(1+229)c = -(1+2^{-29})。精确结果r=260r = 2^{-60}κ=ab/r=(1+229+260)/2602601.15×1018\kappa = |ab|/|r| = (1+2^{-29}+2^{-60})/2^{-60} \approx 2^{60} \approx 1.15 \times 10^{18}。分离乘加相对误差260×253=27=1281\approx 2^{60} \times 2^{-53} = 2^7 = 128 \gg 1,实际结果为0(100%相对误差)。FMA相对误差 = 0(精确结果可表示)。

从ISA的角度来看,FMA指令是所有主流指令集的标配:x86通过FMA3扩展(VFMADD132PD等)支持FMA,ARM从ARMv8.0开始在A64指令集中包含FMADD,RISC-V在F/D浮点扩展中定义了FMADD.S/FMADD.D等四种FMA变体(FMADD, FMSUB, FNMADD, FNMSUB),对应不同的符号组合。

值得注意的是,x86的FMA3编码方式使用了三种操作数排列(132、213、231),分别对应不同的源/目的寄存器映射。例如VFMADD132PD xmm1, xmm2, xmm3计算的是xmm1×xmm3+xmm2\text{xmm1} \times \text{xmm3} + \text{xmm2},结果存入xmm1。这三种排列的存在是由于x86的破坏性操作数约定(结果覆盖一个源操作数),编译器需要根据操作数的活性(liveness)来选择最优的排列以避免额外的寄存器移动。RISC-V的FMA指令则使用独立的目标寄存器(rdrs1rs2rs3\texttt{rd} \neq \texttt{rs1} \neq \texttt{rs2} \neq \texttt{rs3}),不存在这种编码复杂性。

FMA的微架构设计

FMA的微架构设计是浮点单元中最复杂的部分。其核心挑战在于:必须将乘法的精确结果(双倍位宽)与第三个操作数cc进行对阶和加法,而cc的指数可能与乘积的指数相差很大。

图 31.4展示了FMA单元的典型微架构。

FMA单元的典型5级流水线微架构:乘法的部分积压缩与加数$c$的对阶并行执行,然后将$c$注入CSA树进行最终合并。
FMA单元的典型5级流水线微架构:乘法的部分积压缩与加数$c$的对阶并行执行,然后将$c$注入CSA树进行最终合并。

FMA微架构的几个关键设计决策:

加数注入

FMA设计中最巧妙的部分是将加数cc"注入"(inject)到乘法器的CSA压缩树中。在传统的分离设计中,乘法先产生完整的乘积,然后再送入加法器与cc相加——这需要两次CPA(Carry-Propagate Adder)操作:一次在乘法器内部(将CSA树的最终和/进位向量相加),一次在加法器中(将乘积与cc相加)。而在FMA中,对阶移位后的cc被作为CSA树的一个额外输入,与乘法的部分积一起压缩。这意味着乘法和加法在CSA树中被融合——乘法的部分积压缩和cc的加入在同一个硬件结构中完成,最终只需要一次CPA操作

加数注入的具体实现方式取决于cc在CSA树中的注入位置。最常见的策略是在CSA树的最后几层注入cc

  • 末端注入(end-around injection):将cc作为CSA树的最后一层额外输入。此时CSA树的输出从2个向量变为3个向量(和、进位、加数cc),需要额外一层CSA来压缩回2个向量。这种方式对CSA树的修改最小,但增加了1层CSA延迟。

  • 中间注入(mid-tree injection):将cc在CSA树的中间某层注入,替代该层的一个部分积输入。这种方式不增加CSA层级数,但需要修改树的拓扑结构。

  • 独立注入(separate injection):cc不注入CSA树,而是在CSA树输出后通过一个额外的CSA与树的和/进位向量合并。这种方式最灵活,但需要额外的CSA层级。

大多数工业设计采用末端注入或独立注入方式,因为它们对原有乘法器CSA树的修改最小,有利于复用已有的乘法器设计。

加数cc的符号处理

当有效运算为减法(a×bca \times b - c(a×b)+c-(a \times b) + c)时,需要对cc或乘积取补码后再注入CSA树。直接对一个160位宽的数取二进制补码(取反加一)的开销很大。一种常用的优化是使用一的补码(one’s complement)取反代替二的补码——将cc的所有位取反后注入CSA树,然后在最终CPA加法中在最低位加入一个"+1+1"修正。这个+1+1可以通过CPA的进位输入来实现,不增加额外的加法器。

宽数据通路

FMA的内部数据通路必须足够宽以容纳精确的乘积和对阶后的cc。对于双精度FMA,乘积为53×53=10653 \times 53 = 106位,加上cc的对阶可能需要额外的位来覆盖指数差范围,内部数据通路的宽度通常为160\sim170位。这个宽度远超双精度的53位有效位,是FMA面积开销大的主要原因。

具体的位宽分析如下。设双精度尾数为53位(含隐含位),则:

  • 乘积宽度:53×2=10653 \times 2 = 106位(精确乘积)。

  • cc的对阶范围:cc可能比乘积大很多(cc的指数远大于Ea+EbbiasE_a + E_b - \text{bias}),也可能比乘积小很多。当cc远大于乘积时,cc的有效位对齐在乘积的高位之上;当cc远小于乘积时,cc被移到乘积的低位之下(贡献为粘滞位)。

  • 需要覆盖的总宽度:乘积的106位 + cc可能高于乘积的位(最多53位)+ 额外的guard位106+53+3=162\approx 106 + 53 + 3 = 162位。

  • 实际实现通常使用162/8×8=168\lceil 162 / 8 \rceil \times 8 = 168位或直接使用170位,方便字节对齐和布线。

性能分析 6 — FMA内部数据通路宽度的面积影响

FMA内部数据通路的宽度直接影响CSA树、CPA和移位器的面积。对于双精度FMA:

  • CSA树的面积W×Npp\propto W \times N_{\text{pp}},其中W170W \approx 170位为通路宽度,Npp28N_{\text{pp}} \approx 28为部分积数量(27个乘法部分积 + 1个注入的cc)。

  • CPA的面积W\propto W,约170位。

  • 规格化移位器的面积W×log2W\propto W \times \log_2 W,约170×81360170 \times 8 \approx 1360个MUX单元。

相比独立的双精度加法器(数据通路约55位)和乘法器(CSA树约106位),FMA的数据通路宽了约1.6\sim3倍。这使得FMA单元的面积通常是独立加法器和乘法器面积之和的1.5\sim2倍——并非简单的加法器+乘法器,而是一个深度融合的结构。

对阶移位器的挑战

FMA中的对阶移位器面临着比普通加法器更大的挑战。在普通加法器中,对阶移位的范围等于指数差,最大约为2111=20472^{11} - 1 = 2047位(双精度),但实际有效的移位范围被限制在尾数位宽内(约53位),超出这个范围的移位意味着被移位的操作数完全变为零(只贡献粘滞位)。

在FMA中,加数cc需要与a×ba \times b的乘积对齐。乘积的指数范围是Ea+EbbiasE_a + E_b - \text{bias},而cc的指数是EcE_c,因此对阶移位量为(Ea+Ebbias)Ec|(E_a + E_b - \text{bias}) - E_c|。由于乘积的内部表示为106位(双精度),cc的对阶移位范围需要覆盖整个106位的乘积宽度加上额外的保护位——这使得FMA的对阶移位器宽度约为160\sim170位,远大于普通加法器的53位对阶移位器。

对阶方向也更复杂。在普通加法器中,总是将指数较小的操作数向右移位以对齐指数较大的操作数。在FMA中,cc可能比乘积大也可能比乘积小:

  • Ec>Ea+EbbiasE_c > E_a + E_b - \text{bias}时(cc比乘积大),cc的有效位在乘积之上。此时可以将乘积向右移位,但乘积的精确位已经分布在106位中,移位意味着丢弃低位(通过粘滞位保留信息)。

  • Ec<Ea+EbbiasE_c < E_a + E_b - \text{bias}时(cc比乘积小),cc的有效位在乘积之下。此时将cc向右移位以对齐乘积。

  • EcEa+EbbiasE_c \approx E_a + E_b - \text{bias}时,cc与乘积幅度接近,可能发生大规模消去。

实际设计中通常固定将cc进行移位(无论方向),因为乘积的106位已经在CSA树中以部分积形式存在,不适合再移位。cc的对阶移位器宽度为170位,使用7级桶形移位器(log2170=8\lceil \log_2 170 \rceil = 8,但可以优化为7级)。

大规模消去场景

FMA最困难的场景是当乘积a×ba \times b与加数cc的符号相反且幅度非常接近时。此时加法的结果可能发生大规模消去,需要大量的前导零检测和规格化左移。在FMA中,这个问题比普通加法更严重,因为消去可能发生在106位宽的乘积范围内,规格化移位量可能高达106位——需要一个极宽的左移位器。这是FMA硬件面积和延迟开销的主要来源之一。

大规模消去场景的一个具体例子:设a=1.0a = 1.0b=1.0+252b = 1.0 + 2^{-52}(双精度中与1.01.0最近的数),c=1.0c = -1.0。则精确乘积a×b=1.0+252a \times b = 1.0 + 2^{-52},与c=1.0c = -1.0相加后得到2522^{-52}——消去了52位有效位。在FMA的170位内部通路中,这意味着结果的前52位都是零,需要左移52位才能规格化。

FMA的LZA/LZD需要处理170位宽的输入——这比普通加法器的55位LZA宽约3倍,树的深度增加1\sim2层(log2170=8\lceil \log_2 170 \rceil = 8 vs log255=6\lceil \log_2 55 \rceil = 6),面积和延迟都有不可忽视的增加。

双路径FMA

与浮点加法器类似,FMA也可以采用双路径设计:当乘积的指数与cc的指数接近时走Near Path,远时走Far Path。但由于FMA的复杂性,工业界的实践表明单路径FMA加上精心优化的流水线通常是更好的选择——双路径FMA的面积开销过大(两条170位宽的独立数据通路),而且复杂度使验证变得困难。

学术界提出的一种折中方案是"1.5路径"FMA——共享乘法器和CSA树,但在CPA之后分为两条路径:一条用于大规模消去的规格化(宽左移位器),一条用于正常情况的舍入(简单的0/1位移位 + 舍入逻辑)。这种设计减少了面积开销(乘法部分不重复),同时仍能缩短正常情况下的关键路径。

设计提示

FMA的验证是浮点单元中最具挑战性的任务。FMA有三个输入操作数(aabbcc),每个操作数可以是规格化数、非规格化数、零、无穷或NaN五种类型,加上四种舍入模式和四种操作变体(FMADD/FMSUB/FNMADD/FNMSUB),总的输入空间为53×4×4=20005^3 \times 4 \times 4 = 2000种分类组合——每种组合都需要在RTL仿真和形式化验证中覆盖。工业界的FMA验证通常使用数百万条随机测试向量加上数千条手工编写的边界测试用例,验证周期长达数月。

FMA流水线

现代高性能处理器中的FMA单元通常采用4\sim5级流水线,延迟为4\sim5个周期。表表 31.11总结了几款处理器的FMA延迟。

处理器延迟(周期)吞吐量(周期/次)
Intel Skylake\simGolden Cove40.5
AMD Zen 2\simZen 540.5
Apple Firestorm/Everest40.5
ARM Cortex-X441
RISC-V 香山昆明湖51

现代处理器FMA延迟(双精度FMADD)

Intel和AMD在Skylake/Zen 2之后的设计中都配置了两个FMA流水线(分配在两个浮点执行端口上),因此吞吐量为0.5个周期/次(每周期可以启动2个FMA操作)。这意味着对于双精度FMA,每周期可以完成2个乘加操作,即4个浮点操作(FLOP),对于256位AVX向量则是每周期16个双精度FLOP。

FMA的子操作模式

FMA单元不仅可以执行a×b+ca \times b + c,还可以通过将某些输入设置为特殊值来执行纯乘法和纯加法:

  • 乘法a×b+0a \times b + 0,将cc设为+0+0

  • 加法1×b+c1 \times b + c,将aa设为1.01.0

  • 减法1×b+(c)1 \times b + (-c)(1)×b+c(-1) \times b + c

  • 取反乘加(a×b)+c-(a \times b) + c(FNMADD)。

这种设计使得处理器可以用FMA单元替代独立的浮点加法器和乘法器——某些处理器确实这样做了,只配置FMA单元而不配置独立的浮点加法器。但这种选择有其代价,如下一节所讨论的。

性能分析 7 — FMA吞吐量与峰值FLOPS计算

以一个配置了两个256位FMA流水线的处理器为例(如Intel Skylake),在运行双精度FMA指令时的峰值性能为:

  • 每个FMA流水线每周期处理一条256位(4个双精度)FMA指令。

  • 每条FMA指令包含1次乘法和1次加法,即2个FLOP。

  • 两个流水线并行:2×4×2=162 \times 4 \times 2 = 16 双精度FLOP/周期。

  • 在4 GHz频率下:16×4=6416 \times 4 = 64 GFLOPS/核心。

对于512位AVX-512指令(如Intel Skylake-SP),每条FMA指令处理8个双精度值,两个流水线并行提供2×8×2=322 \times 8 \times 2 = 32 FLOP/周期,在4 GHz下达到128 GFLOPS/核心。这也是为什么FMA单元被视为HPC和深度学习处理器中最关键的浮点执行资源。

FMA与独立乘加的权衡

在处理器的浮点执行单元配置中,设计师面临一个重要的架构选择:是使用FMA单元来统一处理所有浮点运算,还是同时配置独立的加法器和乘法器?

设计权衡 2 — FMA统一 vs 独立加法器+乘法器

方案一:仅使用FMA单元

  • 优点:硬件统一,只需要一种功能单元;FMA操作获得最优性能。

  • 缺点:纯加法操作的延迟增加——通过FMA执行加法(1×b+c1 \times b + c)需要4\sim5个周期,而独立加法器只需3个周期。对于加法密集的代码(如归约操作),这1\sim2个周期的额外延迟直接影响依赖链上的性能。

  • 代表:Intel Haswell\simIce Lake的配置——FMA单元同时服务于乘法和加法。

方案二:FMA + 独立加法器

  • 优点:纯加法操作获得最低延迟(3周期);FMA操作也获得最优性能;更多的执行端口提供更高的吞吐量。

  • 缺点:额外的独立加法器增加了面积和功耗。

  • 代表:AMD Zen系列——两个FMA单元 + 两个独立的浮点加法器。

方案三:独立加法器 + 独立乘法器(无FMA)

  • 优点:加法和乘法各获得最优延迟;硬件设计相对简单。

  • 缺点:无法执行FMA指令,对FMA密集的工作负载(如GEMM、深度学习)性能显著受损;不符合IEEE 754-2008对FMA的精度要求。

  • 代表:早期处理器(2010年前)和某些嵌入式处理器。

案例研究 2 — Intel浮点单元配置的演进

Intel浮点执行单元的配置演进很好地展示了FMA与独立加法器之间的架构权衡:

  • Sandy Bridge (2011):1个128位FP ADD + 1个128位FP MUL。不支持FMA指令。加法延迟3周期,乘法延迟5周期。

  • Haswell (2013):2个256位FMA单元(Port 0和Port 1)。所有浮点运算(包括纯加法和纯乘法)都通过FMA单元执行。纯加法的延迟从3周期增加到了FMA的统一延迟(3\sim5周期,取决于具体实现)。

  • Skylake (2015):2个256位FMA单元 + 1个独立的FP ADD单元(Port 1)。Intel认识到纯加法延迟增加对某些工作负载的负面影响,重新增加了独立加法器。纯加法回到4周期延迟。

  • Alder Lake/Golden Cove (2021):继续保持2个FMA + 独立ADD的配置,并在某些新微架构中进一步优化了FMA延迟。

这一演进轨迹说明了一个重要的微架构设计教训:统一化虽然简化了硬件设计,但不一定带来最优的整体性能。当特定操作的使用频率足够高时,为其配置专用的低延迟路径是值得的。

设计提示

在RISC-V浮点执行单元的设计中,浮点延迟配置的选择需要根据目标市场来决定:

  • 嵌入式IoT(如RISC-V RV32F):单精度FMA延迟3\sim4周期即可。面积预算有限,不需要独立加法器。除法可以使用简单的SRT基数-2,延迟30+周期但面积极小。

  • 移动/客户端(如RISC-V RV64FD):双精度FMA延迟4\sim5周期,独立加法器延迟3周期。除法使用SRT基数-4或基数-8,延迟15\sim20周期。

  • HPC/服务器(如RISC-V RV64FDHV):双精度FMA延迟4周期,配置2个FMA端口。除法使用混合SRT+NR方法,延迟13\sim15周期。向量单元支持BF16/FP16混合精度。

RISC-V的ISA模块化特性使得这种差异化的浮点配置在指令集层面就得到了支持——不同的目标选择不同的浮点扩展(F、D、Zfh、V等),硬件实现的复杂度随之调整。

设计提示

在RISC-V处理器设计中,浮点执行单元的配置选择尤为灵活,因为RISC-V的ISA模块化设计允许实现者根据目标应用选择不同的浮点扩展。对于嵌入式应用,可能只实现F扩展(单精度),使用面积较小的单精度FMA单元。对于HPC应用,需要同时实现F+D+V扩展(单精度+双精度+向量),配置多个宽数据通路的FMA单元。香山处理器的浮点单元配置为2个FMA端口+1个独立的FADD端口+1个FDIV/FSQRT端口,兼顾了FMA密集型和加法密集型的工作负载。

浮点除法与平方根

浮点除法和平方根是所有浮点运算中延迟最长的操作。与加法和乘法不同,除法和平方根没有已知的O(logn)O(\log n)时间的精确算法——它们的本质复杂度使得硬件实现必须依赖迭代方法。现代处理器使用两大类迭代算法:减法类算法(如SRT除法,参见第 30.0 章)和乘法类算法(如Newton-Raphson和Goldschmidt)。本节重点讨论乘法类算法在浮点除法和平方根中的应用。

浮点除法和平方根的一个关键特征是其使用频率远低于加法和乘法。在典型的科学计算和机器学习工作负载中,除法操作占浮点指令总数的比例通常低于5%,平方根更低(通常低于1%)。这一统计特性对硬件设计有重要影响:由于使用频率低,为除法和平方根配置全流水线化的高吞吐量执行单元是不划算的——其硬件资源的利用率太低。这就是为什么大多数处理器选择使用非全流水线的迭代单元来实现除法和平方根,将面积预算节省给更常用的加法器和乘法器/FMA单元。

Newton-Raphson迭代法

Newton-Raphson(牛顿-拉弗森)方法通过迭代逼近来计算除法和平方根,每次迭代将精度翻倍(二次收敛)。

除法的Newton-Raphson实现

计算a/ba/b可以转化为求bb的倒数r=1/br = 1/b,然后计算a×ra \times r。求1/b1/b等价于求函数f(x)=1/xbf(x) = 1/x - b的零点。Newton-Raphson迭代公式为:

xi+1=xi(2bxi) x_{i+1} = x_i \cdot (2 - b \cdot x_i)

每次迭代需要两次乘法和一次减法(2bxi2 - b \cdot x_i可以看作一次取反加法)。如果使用FMA单元,这两个操作可以高效地完成。

初始近似值x0x_0通常通过查找表(ROM)获得。查找表使用bb的高几位作为索引,返回1/b1/b的近似值。表的位宽和深度决定了初始精度:一个256项×\times8位的查找表可以提供约8位的初始精度。

更先进的初始近似方法使用分段线性近似(Piecewise Linear Approximation, PLA):将1/b1/b函数在区间[1,2)[1, 2)上分段,每段用一次线性函数ax+bax + b近似。一个64段的PLA可以提供约14位的初始精度,只需要一个小型ROM(存储斜率aa和截距bb)加上一个乘加器。与纯ROM查找表相比,PLA在更小的ROM面积下获得了更高的精度——这正是因为PLA利用了1/b1/b函数的平滑性。

某些处理器还使用二次近似(Quadratic Approximation)来获得更高的初始精度(约24位),但二次近似需要一次平方运算和两次乘加运算来计算ax2+bx+cax^2 + bx + c,硬件复杂度较高。

性能分析 8 — Newton-Raphson除法的迭代次数

Newton-Raphson具有二次收敛性:每次迭代将精度位数翻倍。从初始精度p0p_0位开始,达到目标精度pp位需要的迭代次数为:

n=log2(p/p0)n = \lceil \log_2(p / p_0) \rceil

对于双精度除法(目标精度53位):

  • 初始精度8位(256项查找表):需要 log2(53/8)=3\lceil \log_2(53/8) \rceil = 3 次迭代

  • 初始精度16位(更大的查找表或分段线性近似):需要 log2(53/16)=2\lceil \log_2(53/16) \rceil = 2 次迭代

每次迭代需要2次FMA操作,因此总共需要4\sim6次FMA操作。加上最终的a×ra \times r乘法和舍入校正,双精度Newton-Raphson除法的总延迟约为20\sim30个周期(取决于FMA延迟和调度效率)。

Newton-Raphson的一个重要实现细节是最终舍入校正。迭代收敛后的近似倒数rnr_n虽然精度已经足够,但a×rna \times r_n的结果可能不是正确舍入的(因为rnr_n本身不是精确的1/b1/b)。为了获得IEEE 754要求的正确舍入结果,需要额外的校正步骤:计算余数rem=aq×b\text{rem} = a - q \times b(其中q=a×rnq = a \times r_n是试商),然后根据余数的符号和幅度决定qq是否需要增加或减少1个ULP(Unit in the Last Place)。这个校正步骤通常需要额外的1\sim2次FMA操作。

校正的具体逻辑如下:

  1. 计算试商q=fl(a×rn)q = \text{fl}(a \times r_n)(截断到目标精度,不舍入)。

  2. 使用FMA精确计算余数rem=FMA(b,q,a)=ab×q\text{rem} = \text{FMA}(-b, q, a) = a - b \times q。由于FMA的单次舍入特性,这个余数的计算是精确的(ab×qa - b \times q的精确值可以在双倍位宽内表示,FMA内部保留了足够的精度)。

  3. 根据余数的符号和舍入模式决定最终的商:

    • 如果rem=0\text{rem} = 0qq就是精确的商,直接输出。

    • 如果rem>0\text{rem} > 0:试商偏小,在向上舍入模式下需要q+1ULPq + 1\,\text{ULP}

    • 如果rem<0\text{rem} < 0:试商偏大,在向下舍入模式下需要q1ULPq - 1\,\text{ULP}

    • RNE模式下,需要进一步比较rem|\text{rem}|0.5ULP×b0.5\,\text{ULP} \times b的大小来决定舍入方向。

这个校正过程保证了除法结果的正确舍入(correctly rounded),即结果等于将无限精度的精确商舍入到目标精度后的值。正确舍入是IEEE 754对除法的强制要求,不能用"足够接近"来替代。

初始近似查找表的设计

Newton-Raphson迭代的收敛速度取决于初始近似的精度p0p_0。更高的初始精度意味着更少的迭代次数,但更大的查找表。表表 31.12展示了不同初始精度下查找表的大小和所需的迭代次数。

初始精度查找表条目每条目位宽总ROM大小达到53位的迭代数
8位2568位256 B3
10位102410位1.25 KB3
14位1638414位28 KB2
16位 (PLA)64段32位256 B2
24位 (二次近似)256段72位2.25 KB1

Newton-Raphson除法初始近似的查找表配置

分段线性近似(PLA)是一种面积效率极高的初始近似方法。PLA的原理是将输入范围[1,2)[1, 2)分为KK段(通常K=32K = 326464),每段存储一个线性函数y=aix+biy = a_i x + b_i的系数(ai,bi)(a_i, b_i)。计算初始近似时,用输入bb的高几位选择段号,然后执行一次乘加y=aibfrac+biy = a_i \cdot b_{\text{frac}} + b_i(其中bfracb_{\text{frac}}bb的段内偏移)。PLA需要一个小型ROM(存储KK对系数)和一个乘加器(可以复用FMA单元的一部分),总面积远小于等效精度的直接查找表。

查找表的设计还需要考虑舍入策略。如果查找表返回的初始近似x0x_0的误差恒为正(即x01/bx_0 \geq 1/b),则Newton-Raphson的所有迭代结果都从同一侧收敛——这种单侧收敛简化了最终的校正逻辑。如果查找表的误差可正可负,则迭代可能从两侧振荡收敛,校正逻辑更复杂。大多数实际实现选择单侧收敛的查找表。

平方根的Newton-Raphson实现

计算a\sqrt{a}可以通过先求1/a1/\sqrt{a}(逆平方根),然后乘以aa得到a\sqrt{a}。求1/a1/\sqrt{a}等价于求函数f(x)=1/x2af(x) = 1/x^2 - a的零点。迭代公式为:

xi+1=xi(3axi2)/2=xi2(3axi2) x_{i+1} = x_i \cdot (3 - a \cdot x_i^2) / 2 = \frac{x_i}{2} \cdot (3 - a \cdot x_i^2)

每次迭代需要三次乘法(xi2x_i^2axi2a \cdot x_i^2xiresultx_i \cdot \text{result})和一次减法,同样可以通过FMA来实现。迭代次数与除法相同。

使用FMA指令来表达Newton-Raphson平方根迭代的伪代码如下:

c
// 输入: a (待求平方根的正数)
// 输出: sqrt(a)

double y0 = rsqrt_approx(a);   // 查找表初始近似 1/sqrt(a)
// 第1次迭代: y1 = y0 * (3 - a*y0*y0) / 2
double h  = 0.5 * y0;          // h = y0/2
double r  = fma(-a, y0*y0, 1.0); // r = 1 - a*y0^2
double y1 = fma(h, r, y0*0.5);   // y1 近似
// 第2次迭代 (类似)
double h1 = 0.5 * y1;
double r1 = fma(-a, y1*y1, 1.0);
double y2 = fma(h1, r1, y1*0.5);
// 最终结果
double result = a * y2;  // sqrt(a) = a * (1/sqrt(a))
// + 舍入校正步骤

注意上述代码是高度简化的——实际的硬件或微码实现需要仔细处理特殊值(a=0a = 0a<0a < 0a=a = \inftya=NaNa = \text{NaN})以及最终的舍入校正。

平方根的舍入校正

平方根的舍入校正比除法更为微妙。Newton-Raphson求得的是y1/ay \approx 1/\sqrt{a},最终平方根s=a×ys = a \times y。为了验证ss是否正确舍入,需要检查s2s^2aa的关系:

  1. 计算试平方根s=fl(a×y)s = \text{fl}(a \times y)(截断到目标精度)。

  2. 使用FMA精确计算余数r=FMA(s,s,a)=as2r = \text{FMA}(-s, s, a) = a - s^2

  3. 根据余数的符号判断:

    • 如果r=0r = 0ss是精确的平方根。

    • 如果r>0r > 0ss偏小,s2<as^2 < a,在向上舍入模式下可能需要s+1ULPs + 1\,\text{ULP}

    • 如果r<0r < 0ss偏大,s2>as^2 > a,在向下舍入模式下可能需要s1ULPs - 1\,\text{ULP}

  4. 对于RNE模式,需要进一步判断rr是否恰好在两个候选值的中间(tie),此时选择偶数结果。

这个校正过程保证了FSQRT(a)=fl(a)\text{FSQRT}(a) = \text{fl}(\sqrt{a})——即结果等于将精确平方根舍入到目标精度。平方根校正的一个额外复杂性在于:精确余数的计算as2a - s^2涉及一次平方和一次减法,如果不使用FMA,需要两次独立的运算(两次舍入),可能引入错误。FMA的存在使得余数可以在单次运算中精确计算——这是FMA对除法/平方根实现的重要贡献之一。

除法与平方根的延迟分析

表 31.13分析了基于Newton-Raphson/Goldschmidt方法的双精度除法和平方根的延迟组成。

步骤除法延迟平方根延迟
初始近似查找 (PLA)4周期4周期
第1次NR迭代 (2×\timesFMA串行)8周期8周期
第2次NR迭代 (2×\timesFMA串行)8周期8周期
最终乘法 (a×ra \times ra×ya \times y)4周期4周期
余数计算与舍入校正4\sim5周期4\sim5周期
总计\sim28\sim29周期\sim28\sim29周期

双精度除法/平方根的延迟分解(假设FMA延迟4周期,初始精度14位)

在双FMA端口的处理器上,使用Goldschmidt方法可以将迭代延迟减半(每次迭代4周期而非8周期),总延迟降至\sim16\sim18周期。结合SRT初始近似(省去PLA的延迟),混合方法可以进一步将总延迟降至\sim13\sim15周期——这与表表 31.15中现代处理器的实测数据相吻合。

Goldschmidt除法

Goldschmidt除法是Newton-Raphson的一种变体,其核心思想是:通过反复乘以一个接近1的因子,使除数收敛到1,同时被除数收敛到商。

计算q=a/bq = a / b的Goldschmidt迭代如下。令初始近似 r01/br_0 \approx 1/b,然后:

N0=a,D0=bFi=2DiNi+1=Ni×FiDi+1=Di×Fi\begin{aligned} N_0 &= a, \quad D_0 = b \\ F_i &= 2 - D_i \\ N_{i+1} &= N_i \times F_i \\ D_{i+1} &= D_i \times F_i \end{aligned}

每次迭代,DiD_i越来越接近1,NiN_i越来越接近qq。与Newton-Raphson一样,Goldschmidt也具有二次收敛性。

为了更直观地理解Goldschmidt的收敛过程,考虑计算7/37/3的示例。设初始近似r01/30.333r_0 \approx 1/3 \approx 0.333,则D0=3D_0 = 3

迭代Fi=2DiF_i = 2 - D_iNi+1=Ni×FiN_{i+1} = N_i \times F_iDi+1=Di×FiD_{i+1} = D_i \times F_i
i=0i=023×0.333=1.0012 - 3 \times 0.333 = 1.0017×0.333×1.0012.3367 \times 0.333 \times 1.001 \approx 2.3363×0.333×1.0011.0013 \times 0.333 \times 1.001 \approx 1.001
i=1i=121.001=0.9992 - 1.001 = 0.9992.336×0.9992.3342.336 \times 0.999 \approx 2.3341.001×0.9991.0001.001 \times 0.999 \approx 1.000
i=2i=21.000\approx 1.0002.333\approx 2.3331.000\approx 1.000

经过2\sim3次迭代,NiN_i收敛到7/32.3337/3 \approx 2.333DiD_i收敛到1.0001.000

Goldschmidt除法相对于Newton-Raphson的关键优势在于:每次迭代中的两个乘法Ni×FiN_i \times F_iDi×FiD_i \times F_i相互独立的——它们共享同一个乘数FiF_i,但被乘数不同。这意味着如果处理器有两个FMA单元,这两个乘法可以在同一个周期内同时发射,将每次迭代的延迟从2次FMA操作减少到1次。

Goldschmidt与Newton-Raphson的系统性对比

两种算法具有相同的数学收敛阶——都是二次收敛(每次迭代精度翻倍)。但它们在硬件实现上的差异是显著的。

依赖链结构的差异。 Newton-Raphson的迭代公式xi+1=xi(2bxi)x_{i+1} = x_i \cdot (2 - b \cdot x_i)展开为两步:

  1. ti=bxit_i = b \cdot x_i (乘法,需要xix_i的值)

  2. xi+1=xi(2ti)x_{i+1} = x_i \cdot (2 - t_i) (FMA,需要tit_ixix_i的值)

步骤2依赖步骤1的结果tit_i,因此两步是串行的。每次迭代的延迟=2×TFMA= 2 \times T_{\text{FMA}}

Goldschmidt的迭代展开为三步:

  1. Fi=2DiF_i = 2 - D_i (减法/取反加法)

  2. Ni+1=Ni×FiN_{i+1} = N_i \times F_i (乘法)

  3. Di+1=Di×FiD_{i+1} = D_i \times F_i (乘法)

步骤2和步骤3都只依赖步骤1的结果FiF_i,两者相互独立,可以在两个FMA端口上并行执行。每次迭代的延迟=Tsub+TFMATFMA= T_{\text{sub}} + T_{\text{FMA}} \approx T_{\text{FMA}}(减法延迟被FMA延迟掩盖或可以集成到FMA中)。

收敛速度分析。 设初始近似的误差为ϵ0=x01/b\epsilon_0 = |x_0 - 1/b|。Newton-Raphson的误差递推为:

ϵi+1=bϵi2\epsilon_{i+1} = b \cdot \epsilon_i^2

Goldschmidt的误差递推为:

δi+1=δi2,其中 δi=1Di\delta_{i+1} = \delta_i^2, \quad \text{其中 } \delta_i = 1 - D_i

两者都是二次收敛,但Goldschmidt的误差递推更简洁——没有乘以bb的系数,这使得Goldschmidt的实际收敛速度在b>1b > 1时略优于Newton-Raphson。

表 31.14总结了两种算法在硬件实现中的全面对比。

特性Newton-RaphsonGoldschmidt
收敛阶数二次二次
每次迭代的FMA操作数2(串行)2(可并行)+ 1减法
双FMA端口下每次迭代延迟2×TFMA2 \times T_{\text{FMA}}1×TFMA\approx 1 \times T_{\text{FMA}}
单FMA端口下每次迭代延迟2×TFMA2 \times T_{\text{FMA}}3×TFMA3 \times T_{\text{FMA}}
计算的值1/b1/b(倒数)a/ba/b(直接商)
是否需要最终乘aa是(q=a×(1/b)q = a \times (1/b)否(NiN_i直接收敛到qq
数值稳定性较好需注意DiD_i的舍入累积
硬件复用使用通用FMA使用通用FMA

Newton-Raphson与Goldschmidt除法的硬件实现对比

硬件复用的优势。 Newton-Raphson和Goldschmidt的一个共同优势是它们完全复用已有的FMA硬件——不需要任何专用的除法电路,只需要一个查找表ROM(提供初始近似)和一个微码/控制状态机(控制迭代序列)。这使得在已经配置了FMA单元的处理器中,除法的增量面积开销极小(主要是查找表和控制逻辑,通常不到FMA单元面积的5%)。

Goldschmidt的数值稳定性问题。 与Newton-Raphson相比,Goldschmidt有一个微妙的数值稳定性问题:每次迭代涉及三个乘法和一个减法,这些运算如果在有限精度下执行(即每步都进行舍入),累积的舍入误差可能导致最终结果不是正确舍入的。Newton-Raphson只计算倒数1/b1/b,最终的a×(1/b)a \times (1/b)可以精确地在一步FMA中完成并正确舍入。而Goldschmidt的NiN_i在每次迭代中都经过舍入,最终的NnN_n可能存在大于0.5 ULP的误差。

为了解决这个问题,Goldschmidt除法的实际实现通常采用扩展内部精度——在迭代过程中使用比目标精度更宽的数据通路(如双精度除法使用68\sim70位内部尾数),确保累积舍入误差不超过0.5 ULP。最终结果再截断到目标精度并执行正确舍入。

设计提示

Goldschmidt除法特别适合具有两个FMA流水线的处理器。在Intel和AMD的现代设计中,两个FMA端口可以在同一周期内分别执行Ni×FiN_i \times F_iDi×FiD_i \times F_i,使得每次Goldschmidt迭代只需要4\sim5个周期(一次FMA延迟 + 调度开销)。相比之下,Newton-Raphson的每次迭代中两个乘法是串行依赖的(xi(2bxi)x_i \cdot (2 - b \cdot x_i)中的第二个乘法依赖第一个乘法的结果),无法并行。

选择建议:在配置了双FMA端口的高性能处理器中,Goldschmidt是更好的选择(如IBM POWER系列)。在只有单个FMA端口的设计中(如某些嵌入式处理器),Newton-Raphson因其更少的操作总数而更优。在使用SRT+迭代混合方案的设计中(如Intel Skylake),Newton-Raphson的最终精化步骤更简单直接。

共享除法平方根单元

浮点除法和平方根有一个重要的共同点:它们都可以通过乘法迭代来实现。这使得处理器可以设计一个共享的除法/平方根执行单元,复用同一套乘法器和迭代控制逻辑。

在基于乘法迭代的实现中,除法和平方根的差异仅在于:

  • 初始近似:除法使用1/b1/b的查找表,平方根使用1/a1/\sqrt{a}的查找表。

  • 迭代公式:除法使用式式 (31.13)或式式 (31.15)\sim式 (31.17),平方根使用式式 (31.14)

  • 最终乘法:除法为a×(1/b)a \times (1/b),平方根为a×(1/a)a \times (1/\sqrt{a})

这些差异可以通过多路选择器和少量的控制逻辑来处理,而主要的硬件资源(乘法器、CSA树、查找表ROM)在除法和平方根之间完全共享。

Goldschmidt方法也可以用于平方根的计算。Goldschmidt平方根的迭代公式为:

Y0=a,H01/(2a)Ri=0.5Yi×HiYi+1=Yi+Yi×RiHi+1=Hi+Hi×Ri\begin{aligned} Y_0 &= a, \quad H_0 \approx 1/(2\sqrt{a}) \\ R_i &= 0.5 - Y_i \times H_i \\ Y_{i+1} &= Y_i + Y_i \times R_i \\ H_{i+1} &= H_i + H_i \times R_i \end{aligned}

其中YiY_i收敛到a\sqrt{a}HiH_i收敛到1/(2a)1/(2\sqrt{a})。与Goldschmidt除法类似,每次迭代中式式 (31.20)和式式 (31.21)相互独立的——可以在双FMA端口上并行执行。Goldschmidt平方根的收敛速度也是二次的。

共享除法/平方根单元的控制状态机通过以下方式区分两种运算:

  • 操作码(FDIV或FSQRT)控制初始近似查找表的选择(1/b1/b的表或1/a1/\sqrt{a}的表)。

  • 迭代公式中的系数不同(除法使用2Di2 - D_i,平方根使用0.5YiHi0.5 - Y_i H_i),通过MUX选择。

  • 迭代次数可能不同(平方根通常比除法多一次迭代,因为平方根需要更高的内部精度来保证正确舍入)。

  • 最终校正步骤不同(除法计算余数aqba - qb,平方根计算余数as2a - s^2)。

在另一种实现方式中,处理器不配置专用的除法/平方根硬件,而是通过微码将除法和平方根分解为一系列FMA操作。这种软件化的方式节省了专用硬件的面积,但延迟更长(因为需要排队等待FMA单元)且占用了FMA单元的带宽。在面积极度敏感的嵌入式处理器中,微码实现是一种合理的选择;而在高性能处理器中,专用的迭代硬件(可能与SRT结合)是标准配置。

SRT与乘法迭代的混合实现

许多现代处理器采用的是SRT与乘法迭代的混合方法。以Intel Skylake为例,其浮点除法的实现分为两个阶段:

  1. SRT阶段:使用高基数SRT除法器(基数-16或更高)产生一个初始商的近似值,精度约为14\sim16位。SRT阶段使用专用的商位选择查找表和迭代余数寄存器。

  2. Goldschmidt/NR精化阶段:以SRT的结果作为初始近似,执行1\sim2次Goldschmidt或Newton-Raphson迭代来提升精度到目标值(53位)。这个阶段复用FMA单元。

这种混合方法结合了SRT的低延迟初始近似(不需要查找表ROM)和乘法迭代的快速收敛,在延迟和面积之间取得了良好的平衡。

混合方法的优势可以量化如下。纯SRT基数-16方法(每迭代产生4位商)需要53/4=14\lceil 53/4 \rceil = 14次迭代,延迟约14\sim16周期。纯Newton-Raphson方法(初始精度14位,FMA延迟4周期)需要log2(53/14)=2\lceil \log_2(53/14) \rceil = 2次迭代,每次迭代8周期(双FMA串行),加上初始近似和校正,总延迟约24周期。而混合方法:SRT阶段4次基数-16迭代(\sim4\sim5周期)产生16位初始近似,然后1次Goldschmidt迭代(\sim4周期)精化到32位,再1次迭代精化到53位(\sim4周期),最终校正(\sim2周期),总延迟\sim14\sim15周期——优于两种纯方法。

混合方法的另一个优势是资源复用的灵活性。SRT阶段使用的专用硬件(商位选择表、余数寄存器、加法器)面积很小;精化阶段复用已有的FMA单元,增量面积近乎为零。如果FMA单元正忙于其他操作,精化阶段可以等待——这在除法使用频率较低的场景下是完全可以接受的。

变精度除法优化

一个细节但重要的优化是变精度除法——根据目标精度(单精度24位 vs 双精度53位)调整迭代次数。对于单精度除法,Newton-Raphson从14位初始近似只需要1次迭代即可达到24位精度(14×2=28>2414 \times 2 = 28 > 24);对于双精度则需要2次迭代。这意味着单精度除法的延迟可以显著短于双精度——在Intel Skylake上,单精度FDIV延迟为11周期,而双精度为13\sim14周期。

类似地,半精度(FP16)除法只需要11位精度,从14位初始近似甚至不需要迭代——初始近似已经足够。半精度除法的延迟可以低至6\sim8周期。这种变精度优化使得低精度除法获得了接近加法的延迟,在深度学习推理等大量使用低精度运算的场景中是有价值的。

除法/平方根的非阻塞实现

由于除法和平方根的延迟很长(13\sim21个周期),如果在执行期间阻塞整个浮点流水线,将严重影响后续指令的吞吐量。因此,现代处理器的除法/平方根单元通常是非阻塞(non-blocking)的——除法在专用的迭代单元中执行,不占用FMA流水线的执行槽位(除非使用微码实现)。其他浮点指令可以在除法执行期间继续正常发射和执行。除法完成后,结果通过公共数据总线(CDB)写回,与其他功能单元共享写回端口。

非阻塞实现的微架构细节包括:

  • 独立的迭代寄存器:除法/平方根单元拥有自己的内部状态寄存器(存储当前迭代的中间结果NiN_iDiD_ixix_i),与主流水线的流水线寄存器分离。这使得主流水线可以在除法迭代期间继续处理其他指令。

  • 写回端口仲裁:当除法完成时,它需要一个写回端口来将结果广播到CDB。如果该端口正被另一个功能单元使用(如FMA),除法的写回需要等待一个周期。在高吞吐量的处理器中,通常为除法配置一个低优先级的专用写回端口,或允许除法在所有端口繁忙时延迟1\sim2个周期写回。

  • 分步(pipelined)SRT:对于基于SRT的除法实现,SRT迭代器本身可以是部分流水线化的——允许两个除法操作在SRT单元中重叠执行(前一个除法在SRT的后半段迭代时,后一个除法可以开始其前半段迭代)。这种分步执行使得除法的吞吐量提升到\sim4周期/次(而非等待前一个除法完全完成后再开始下一个)。

  • 微码实现的资源占用:如果除法通过微码分解为FMA操作序列来实现,则除法在执行期间会占用FMA端口。这意味着在除法执行的4\sim6次FMA迭代期间,FMA端口无法服务其他FMA/FMUL/FADD指令,显著降低了浮点吞吐量。这是专用除法硬件vs微码实现的核心权衡——专用硬件不占用FMA端口,而微码实现虽然节省面积但影响吞吐量。

案例研究 3 — RISC-V处理器中的除法实现策略

RISC-V处理器的浮点除法实现策略因目标应用的不同而差异很大:

  • BOOM(Berkeley Out-of-Order Machine):使用SRT基数-4除法器,延迟约为53/2+3=3053/2 + 3 = 30周期(双精度,每次迭代产生2位商)。非阻塞,独立于FMA流水线。

  • 香山处理器:使用SRT基数-16除法器,延迟约为14\sim16周期。配置在独立的FDIV端口上,不影响FMA吞吐量。

  • Rocket(顺序RISC-V核心):使用简单的SRT基数-2除法器,延迟较长但面积极小。由于是顺序核心,阻塞执行不影响乱序调度(因为没有乱序调度)。

  • CVA6(前Ariane):通过迭代器实现除法,可配置为SRT或Newton-Raphson方式,面积和延迟可根据目标应用调整。

这种实现的多样性正是RISC-V模块化ISA设计理念的体现——ISA只定义了FDIV.S/FDIV.D指令的语义(正确舍入的除法),而不规定实现方式。

处理器FDIV延迟FSQRT延迟实现方式
Intel Skylake\simGolden Cove13\sim1418\sim19SRT + 后续迭代
AMD Zen 2\simZen 513\sim1515\sim21基数-16 SRT
Apple Firestorm/Everest10\sim1412\sim17未公开
ARM Cortex-X414\sim1716\sim21未公开

现代处理器浮点除法和平方根延迟(双精度)

需要注意的是,浮点除法和平方根通常是非全流水线(not fully pipelined)的——它们无法每周期启动一次新操作。典型的吞吐量为4\sim8个周期/次(即每4\sim8个周期才能启动一次新的除法或平方根)。这使得除法和平方根成为处理器中"最慢"的浮点运算,编译器和算法设计者通常尽量避免使用除法(例如用a×(1/b)a \times (1/b)代替a/ba/b,其中1/b1/b可以预计算)。

设计权衡 3 — SRT除法 vs 乘法迭代 vs 混合方法

  • SRT除法:每次迭代产生固定位数的商(基数-4产生2位/迭代,基数-16产生4位/迭代)。优点是不需要乘法器,使用加法器和查找表即可实现。缺点是延迟与精度成线性关系——双精度需要约53/4=1453/4 = 14次基数-16迭代。适用于面积敏感的设计。

  • 纯乘法迭代(Newton-Raphson/Goldschmidt):二次收敛使得迭代次数仅为O(loglogn)O(\log \log n)。但需要乘法器(或FMA),初始近似需要查找表ROM。延迟取决于乘法器延迟和迭代次数的乘积。适用于已经配置了高性能乘法器的设计。

  • 混合方法:先用SRT产生中等精度(14\sim16位)的初始近似,再用1\sim2次乘法迭代精化到目标精度。结合了两种方法的优点:SRT阶段不需要乘法器,迭代阶段利用已有的FMA单元。Intel和AMD的现代设计大多采用此方法。

案例研究 4 — 快速逆平方根:从软件Hack到硬件支持

1999年,Quake III Arena游戏引擎中出现了一段著名的"快速逆平方根"代码,通过将浮点数的位模式重新解释为整数来获得1/x1/\sqrt{x}的初始近似:

c
float Q_rsqrt(float number) {
    long i;
    float x2, y;
    x2 = number * 0.5F;
    y  = number;
    i  = *(long *)&y;           // 位模式 -> 整数
    i  = 0x5f3759df - (i >> 1); // "魔数"近似
    y  = *(float *)&i;          // 整数 -> 位模式
    y  = y * (1.5F - (x2*y*y)); // 一次NR迭代
    return y;
}

这段代码利用了IEEE 754浮点编码的数学特性:将浮点位模式视为整数时,log2\log_2近似为该整数值的线性函数。这个"魔数"技巧提供了约3\sim4位的初始精度,随后一次Newton-Raphson迭代将精度提升到足以满足3D图形渲染的需求(约7\sim8位)。

现代处理器已经通过专用的逆平方根估计指令(如x86的RSQRTPS和ARM的FRSQRTE)在硬件中提供了类似的功能,典型精度为11\sim12位,延迟仅为4\sim5个周期——远快于完整的FSQRT指令。这些估计指令加上1\sim2次Newton-Raphson迭代(使用FMA实现)可以达到单精度或双精度的完整精度,总延迟仍然远小于FSQRT指令。

类似地,x86的RCPPS指令提供了1/x1/x的快速近似,ARM的FRECPE指令提供了类似的功能。RISC-V目前没有标准的近似倒数/逆平方根指令,但正在讨论中的Zfa扩展可能包含此类指令。从微架构角度看,这些近似估计指令通常使用查找表ROM实现(与Newton-Raphson迭代的初始近似共享ROM),只需要一个ROM查找加上少量的线性插值逻辑,延迟和面积开销都很小。

浮点异常处理

IEEE 754标准定义了五种浮点异常,处理器的浮点单元必须能够正确地检测并报告这些异常。异常处理的硬件实现直接影响浮点运算的正确性和性能。

异常标志的产生

IEEE 754定义的五种浮点异常如下:

  1. Invalid Operation(无效操作):当运算没有数学意义时触发。典型情况包括:0/00/0\infty - \infty×0\infty \times 0x\sqrt{x}x<0x < 0)、对sNaN的任何运算。默认结果为qNaN。

  2. Division by Zero(除以零):当非零有限数除以零时触发。注意0/00/0不是Division by Zero(而是Invalid Operation)。默认结果为±\pm\infty

  3. Overflow(上溢):当运算的精确结果的幅度超过目标格式可表示的最大有限数时触发。例如单精度中,当结果>3.4028235×1038> 3.4028235 \times 10^{38}时。默认结果取决于舍入模式:RNE模式下结果为±\pm\infty,RZ模式下结果为最大有限数。

  4. Underflow(下溢):当运算的精确结果为非零但其幅度小于目标格式的最小规格化数时触发。单精度中,这意味着结果<21261.18×1038< 2^{-126} \approx 1.18 \times 10^{-38}。Underflow的检测有两种模式:

    • 舍入前检测(before rounding):如果精确结果在规格化范围之外,则触发Underflow——无论舍入后结果是否回到规格化范围。

    • 舍入后检测(after rounding):只有当舍入后的结果仍然在非规格化范围内时才触发。IEEE 754允许两种检测方式,但大多数现代处理器使用舍入后检测。

  5. Inexact(不精确):当运算的舍入结果与精确结果不同时触发。这是最常见的异常——几乎所有的除法和大多数的乘法和加法都会触发Inexact异常。Inexact通常不会引起中断,只是设置状态标志位。

在硬件实现中,异常检测逻辑分布在浮点流水线的各个阶段:

  • 输入检查(第1级):检测操作数是否为NaN、±\pm\infty±0\pm 0或非规格化数。如果输入为sNaN,立即标记Invalid Operation。如果是0/00/0/\infty/\infty,标记Invalid Operation。

  • 指数检查(中间级):在指数加法或减法后检查是否超出有效范围,预判Overflow或Underflow。

  • 舍入后检查(最后级):检查舍入是否改变了结果(Inexact);检查舍入后的指数是否溢出或下溢(确认Overflow/Underflow)。

异常标志通过流水线寄存器逐级传递,在指令提交时写入浮点状态寄存器(如x86的MXCSR或RISC-V的fcsr)。在乱序执行的处理器中,异常标志的写入必须是按序的——只有当指令成功提交(commit)时才更新状态寄存器,推测执行的指令即使检测到异常也不能提前修改状态标志。

异常检测的时序约束

异常检测的一个重要时序约束是:某些异常只有在运算完成后才能确定,而不能在运算进行中提前判定。具体来说:

  • Invalid Operation:大多数情况可以在输入分类时(第1级)立即检测——如0/00/0\infty - \infty、sNaN输入等。但某些Invalid Operation(如x\sqrt{-x}xx为负数)需要检查符号位,也可以在第1级完成。

  • Division by Zero:在除法的第1级检测(检查除数是否为零且被除数非零)。

  • Overflow:只有在规格化和指数计算完成后才能确定——需要检查最终指数是否超出最大值。这通常在最后1\sim2级完成。

  • Underflow:类似Overflow,需要在规格化后检查指数是否低于最小值。如果使用"舍入后检测"模式,还需要等到舍入完成。

  • Inexact:需要在舍入决策完成后才能确定(检查G/R/S位是否全零)。这是最后确定的异常——通常在最后一级。

这些时序约束意味着异常标志的确定是渐进的——从第1级开始逐步确定各种异常标志,到最后一级才能确定全部五种异常的状态。异常标志通过流水线寄存器逐级传递并累积,最终在流水线输出端得到完整的异常向量。

异常屏蔽与陷阱

IEEE 754标准要求每种异常都可以被屏蔽(masked)或不屏蔽(unmasked)。当异常被屏蔽时,处理器使用默认结果(如NaN、±\pm\infty等)继续执行,仅设置状态标志位。当异常不被屏蔽时,触发一个陷阱(trap)到操作系统或用户定义的异常处理程序。

在实践中,几乎所有的浮点异常都被屏蔽运行——不屏蔽的浮点异常会导致频繁的陷阱(特别是Inexact异常几乎每条指令都会触发),严重影响性能。不屏蔽的浮点异常主要用于调试和数值分析,不用于生产环境。

然而,不屏蔽异常的存在对浮点流水线的设计有重要影响:当异常不被屏蔽时,处理器不能使用默认结果继续执行后续依赖指令。在乱序处理器中,这意味着检测到不屏蔽异常的指令必须等到提交阶段才能确认是否真正触发陷阱(因为之前的推测可能被取消),这增加了异常处理的复杂度。大多数处理器通过在提交阶段检查异常标志和屏蔽位来决定是否触发陷阱。

硬件描述 4 — 浮点异常检测的流水线集成

以下伪代码展示了浮点加法器中异常检测逻辑的集成方式(简化示意):

verilog
// Stage 1: 输入分类
logic a_is_nan, a_is_inf, a_is_zero, a_is_denorm;
logic b_is_nan, b_is_inf, b_is_zero, b_is_denorm;

always_comb begin
  a_is_nan    = (a_exp == '1) && (a_mant != 0);
  a_is_inf    = (a_exp == '1) && (a_mant == 0);
  a_is_zero   = (a_exp == 0) && (a_mant == 0);
  a_is_denorm = (a_exp == 0) && (a_mant != 0);
  // ... b 同理
end

// 特殊情况快速路径
logic invalid_op, result_is_nan;
always_comb begin
  invalid_op = 0;
  // inf - inf = NaN (Invalid)
  if (a_is_inf && b_is_inf && (op == SUB)
      && (a_sign == b_sign))
    invalid_op = 1;
  // sNaN input
  if (a_is_snan || b_is_snan)
    invalid_op = 1;
end

// Stage 3: 舍入后异常检测
logic overflow, underflow, inexact;
always_comb begin
  overflow  = (result_exp >= exp_max);
  underflow = (result_exp <= 0) && !result_is_zero;
  inexact   = (guard | round_bit | sticky);
end

非规格化数的处理

非规格化数(denormal / subnormal number)的硬件处理是浮点单元设计中最令人头痛的问题之一。非规格化数的存在打破了规格化数中"隐含位为1"的假设,导致常规的浮点运算逻辑需要大量的特殊处理。

为什么非规格化数的处理困难?

要理解非规格化数为什么对硬件造成巨大困扰,需要先回顾规格化数的一个关键假设:所有规格化浮点数的有效数字最高位为1。这个假设使得浮点运算单元的设计可以大幅简化——加法器的对阶逻辑、乘法器的部分积生成、以及规格化移位器都假定输入的有效位以"1.xxx1.\text{xxx}"开头。非规格化数的有效位以"0.xxx0.\text{xxx}"开头,这意味着:

(1)有效位的前导零不确定。 规格化数的有效位前导零始终为0(因为最高位是隐含的1)。而非规格化数的有效位0.M0.M可能有0到p2p-2个前导零(其中pp为有效位总位宽)。例如,单精度非规格化数0.000000000000000000000012×21260.00000000000000000000001_2 \times 2^{-126}有22个前导零,仅有1位有效精度。这意味着非规格化数的实际精度是可变的——从1位到p1p-1位不等。

(2)需要额外的前导零计数和预规格化移位。 当一个非规格化操作数参与运算时,为了正确对阶,硬件需要首先确定该操作数的实际有效位位置。这需要对非规格化数的尾数执行前导零计数(LZC),然后左移到规格化位置。这个"输入预规格化"操作在正常的规格化数流水线中不存在——它需要一个额外的LZC电路和一个左移位器,而且这个操作位于流水线的最前端,直接增加了第一级的关键路径。

(3)指数系统的不一致性。 规格化数的实际指数 = EbiasE - \text{bias}E1E \geq 1),而非规格化数的实际指数固定为1bias1 - \text{bias}(无论EE字段的值如何,因为E=0E = 0被特殊定义)。这个1bias1 - \text{bias}而非0bias0 - \text{bias}的选择是为了使非规格化数与最小规格化数之间没有间隙(渐进下溢),但它使得指数比较逻辑需要增加一个特殊情况分支。

非规格化数对各种运算的具体影响

以浮点加法为例,非规格化数带来的额外复杂度包括:

  1. 隐含位不再为1:规格化数的有效位是1.M1.M,而非规格化数的有效位是0.M0.M。加法器的尾数拼接逻辑需要根据指数是否为零来选择隐含位的值。

  2. 指数修正:非规格化数的实际指数是1bias1 - \text{bias}(而非0bias0 - \text{bias}),这意味着指数比较逻辑需要特殊处理零指数的情况。

  3. 规格化后的非规格化:运算结果可能是非规格化数——当结果的指数下溢到零以下时,需要将尾数右移并将指数固定为零。这个"非规格化结果"的处理增加了规格化阶段的复杂度。

  4. 双重舍入风险:非规格化结果的右移可能需要额外的舍入操作,而这个舍入与常规舍入的交互需要仔细处理以避免双重舍入错误。

对于浮点乘法,非规格化数的影响同样严重:

  1. 乘积指数的下溢检测:两个规格化数的乘积可能产生非规格化结果(当指数之和足够小时)。指数通路需要检测下溢并触发结果非规格化流程。

  2. 非规格化输入的乘法:如果一个输入是非规格化数,其有效位不是1.M1.M而是0.M0.M。乘法器的Booth编码逻辑需要处理最高位为0的输入——这本身不难,但乘积的规格化移位量变得不可预测,可能需要多位左移。

  3. 乘积的非规格化输出:当乘积指数下溢时,需要将乘积尾数右移直到指数变为0。右移量可能很大(最多达尾数位宽),需要一个大型右移位器——这在正常乘法流水线中不存在。

对于FMA运算,非规格化数的问题尤为复杂:乘积和加数cc都可能是非规格化数,而且FMA的中间结果(106位乘积)也可能需要非规格化输出。FMA的内部数据通路已经很宽(160\sim170位),非规格化处理使得移位器和LZD电路的宽度进一步增大。

性能分析 9 — 非规格化数在实际工作负载中的出现频率

虽然非规格化数在大多数计算中极为罕见,但在某些特定场景下会频繁出现:

  • 音频处理:当音频信号的幅度逐渐衰减到接近零时(如混响尾部、淡出效果),浮点采样值进入非规格化范围。一个48 kHz的音频流每秒产生48000个采样,如果其中1%是非规格化数,则每秒触发480次微码辅助。

  • 数值仿真:有限元分析中,当网格单元的刚度矩阵接近奇异时,某些中间计算可能产生极小的非规格化值。

  • 机器学习训练:在低精度(FP16/BF16)训练中,梯度值在深层网络中逐渐缩小,容易进入非规格化范围。BFloat16格式的最小规格化值为21261.18×10382^{-126} \approx 1.18 \times 10^{-38},但其非规格化范围更窄(仅有7位尾数),使得进入非规格化的阈值更高。

  • 物理仿真:粒子物理的蒙特卡罗模拟中,概率权重可能跨越极大的动态范围,偶尔产生非规格化值。

三种处理策略

现代处理器对非规格化数的处理采用三种不同的策略:

策略一:完全硬件支持。 在浮点流水线的每个阶段都集成非规格化数的处理逻辑,确保非规格化数的操作数和结果都能在正常的流水线延迟内完成。这种策略提供了最好的性能一致性,但增加了浮点单元的面积和设计复杂度。Apple的处理器和Intel的Alder Lake之后的设计逐步朝这个方向发展。

策略二:微码辅助。 当浮点流水线检测到非规格化操作数或结果时,触发一个微码辅助(microcode assist)——取消当前操作,转入微码ROM执行一段辅助例程来处理非规格化数。这段微码通常将非规格化数转换为规格化数(通过乘以2的适当幂次),在规格化域中执行运算,然后将结果转换回非规格化形式。微码辅助的延迟极高——通常需要100\sim200个周期,但由于非规格化数在大多数工作负载中极为罕见,这种策略在性能上通常是可接受的。Intel Skylake和AMD Zen系列对非规格化数的输入使用微码辅助。

策略三:Flush-to-Zero / Denormals-Are-Zero。 提供一个模式位(如x86 MXCSR中的FTZ和DAZ位),允许软件选择将非规格化结果冲刷为零(FTZ),或将非规格化输入视为零(DAZ)。这完全绕过了非规格化数的处理,但不符合严格的IEEE 754标准。在性能敏感的应用(如游戏、音视频处理、深度学习推理)中,FTZ/DAZ几乎总是被启用的。

性能分析 10 — 非规格化数的性能影响

在Intel Skylake处理器上,涉及非规格化数的浮点操作的性能惩罚如下:

操作正常延迟非规格化延迟
FADD(非规格化输入)3\sim4周期\sim150周期
FMUL(非规格化结果)4周期\sim120周期
FMA(非规格化输入+结果)4周期\sim170周期

非规格化数导致的性能下降可达30\sim40倍。在音频处理中,当信号幅度趋近于零时,大量的采样值可能进入非规格化范围,导致CPU使用率骤增——这是音频软件中一个著名的性能陷阱。解决方案通常是在混音链的最后注入极小的DC偏移(约102010^{-20})来防止信号进入非规格化范围,或者启用FTZ模式。

设计提示

随着AI和HPC工作负载对浮点性能的要求不断提高,越来越多的处理器正在增加对非规格化数的硬件支持,以消除微码辅助带来的不可预测的性能波动。ARM的Cortex-A系列从A76开始对单精度非规格化数提供了完全的硬件支持,双精度非规格化数则在后续设计中逐步支持。RISC-V的浮点扩展(F/D/Q/Zfh)要求实现支持非规格化数,但允许实现选择是通过硬件还是通过陷阱来处理。对于面向HPC的设计,完全的硬件支持正在成为标准配置。

非规格化数的硬件支持实现细节

对于选择在硬件中完全支持非规格化数的设计,主要的修改集中在浮点流水线的输入、运算和输出三个阶段。

输入阶段的修改:预规格化流程。

非规格化输入的处理需要一个预规格化(pre-normalization)步骤——将非规格化数转换为一种"扩展规格化"的内部表示,使后续的运算逻辑不需要区分规格化数和非规格化数。具体流程为:

  1. 非规格化检测:检查指数字段是否为全零(AND归约的取反)。如果指数为零且尾数不为零,则标记该操作数为非规格化数。

  2. 前导零计数(LZC):对非规格化数的尾数0.M0.M执行前导零计数,得到前导零数量kk。这个LZC电路与31.2.3 节中的LZD结构完全相同,可以共享硬件。对于单精度(23位尾数),LZC需要处理23位输入,输出5位计数值。

  3. 尾数左移:将尾数左移kk位,使最高有效位变为1,形成规格化形式的有效位1.M1.M'

  4. 指数修正:将指数设置为1k1 - k(在偏移编码下)。注意修正后的指数可能是负数——在内部使用扩展位宽的有符号指数来表示。

这个预规格化流程的关键路径为LZC延迟 + 移位器延迟O(logn)+O(logn)=O(logn)\approx O(\log n) + O(\log n) = O(\log n)。这个延迟被添加到流水线的第一级,可能增加该级的关键路径。一种优化是将预规格化拆分为一个独立的半级(half-stage),与正常的输入解包并行执行——对于规格化输入,预规格化路径不在关键路径上;对于非规格化输入,使用预规格化结果。

硬件描述 5 — 非规格化输入的预规格化逻辑

verilog
module denorm_prenormalize (
  input  logic [10:0] biased_exp,     // 11位偏移指数
  input  logic [51:0] mantissa,       // 52位尾数字段
  output logic [52:0] norm_mantissa,  // 53位规格化有效位(含隐含位)
  output logic [10:0] adjusted_exp,   // 修正后的指数
  output logic        is_denorm
);
  logic [5:0] lzc_count;   // 前导零数量(最多52)
  logic [51:0] shifted_mantissa;

  // 检测非规格化数
  assign is_denorm = (biased_exp == 11'b0) && (mantissa != 52'b0);

  // 对尾数执行前导零计数
  lzc_52bit u_lzc (
    .data_in  (mantissa),
    .lzc_out  (lzc_count)
  );

  // 左移尾数消除前导零
  assign shifted_mantissa = mantissa << lzc_count;

  // 输出: 规格化有效位和修正指数
  always_comb begin
    if (is_denorm) begin
      // 非规格化数: 左移后, 有效位变为1.shifted_mantissa[50:0]
      norm_mantissa = {1'b1, shifted_mantissa[51:0]};
      // 指数修正: 1 - lzc_count (以偏移编码)
      // 原始biased_exp=0 -> 实际指数 = 1 - bias
      // 左移lzc_count位 -> 实际指数 = 1 - bias - lzc_count
      // 偏移编码 = 1 - lzc_count (可能为负, 用扩展位宽)
      adjusted_exp = 11'd1 - {5'b0, lzc_count};
    end else begin
      // 规格化数: 正常解包
      norm_mantissa = {1'b1, mantissa};
      adjusted_exp  = biased_exp;
    end
  end
endmodule

输出阶段的修改:结果非规格化流程。

当运算结果的指数下溢(指数0\leq 0时),需要将结果从规格化形式转换为非规格化形式。这个流程是预规格化的"逆过程":

  1. 下溢检测:检查规格化后的结果指数是否0\leq 0。下溢量d=1Eresultd = 1 - E_{\text{result}}给出了需要右移的位数。

  2. 尾数右移:将结果尾数右移dd位。右移量dd的范围为11p1p-1pp为有效位总数),因此需要一个大型右移位器。对于双精度,dd最大为52,需要一个52位的桶形右移位器。

  3. 粘滞位更新:右移过程中被移出的位需要OR归约为新的粘滞位,与原有的G/R/S位合并。

  4. 重新舍入:使用更新后的G/R/S位执行舍入。注意,这次舍入发生在常规舍入之后(结果已经经过了一次规格化和舍入),如果不仔细处理,可能导致双重舍入错误。

  5. 指数固定:将结果指数设置为0(偏移编码),隐含位设置为0。

避免双重舍入的关键技术。 双重舍入的根本问题是:第一次舍入(在规格化结果上)可能丢失了第二次舍入(在非规格化结果上)所需要的精度信息。解决方案是延迟舍入——先检测是否会产生非规格化结果,如果是,则跳过第一次舍入,先执行右移,然后只进行一次舍入。这需要在规格化阶段就预判指数是否会下溢——可以通过比较指数与0的大小来实现。

具体实现中,流水线需要同时计算两个版本的结果:

  • 规格化版本:假设结果是规格化数,执行正常的规格化和舍入。

  • 非规格化版本:假设结果是非规格化数,执行右移和一次性舍入。

然后根据指数是否下溢来选择正确的版本。这种"双版本"策略避免了双重舍入,但增加了面积(需要额外的右移位器和舍入逻辑)。

完全硬件支持非规格化数的面积开销约为浮点单元总面积的5%\sim10%,关键路径的增加通常可以控制在0.5\sim1个FO4门延迟以内(在精心设计的情况下不影响时钟频率)。对于高性能设计来说,这个开销换来的是消除了100\sim200周期的微码辅助惩罚——在非规格化数频繁出现的工作负载中,性能收益是显著的。

浮点异常在乱序处理器中的精确处理

在乱序执行的超标量处理器中,浮点异常的精确处理引入了额外的复杂性。考虑以下场景:一条浮点除法指令DD产生了Division by Zero异常(不被屏蔽),但在DD之前有一条条件分支BB尚未解析。如果BB的预测是错误的,DD是推测执行的结果,不应触发异常。

为了正确处理这种情况,浮点异常标志被存储在ROB的entry中(与指令关联),而非立即写入状态寄存器。异常标志的处理过程如下:

  1. 检测阶段:浮点执行单元在运算过程中检测到异常条件,将异常标志位写入结果的tag字段,随结果一起通过CDB广播。

  2. 暂存阶段:ROB接收到带有异常标志的结果后,将标志存储在对应entry的异常字段中。此时不采取任何操作。

  3. 提交阶段:当指令到达ROB的提交指针时,检查其异常字段。如果有不被屏蔽的异常,触发异常处理:刷新流水线、保存精确的体系结构状态、跳转到异常处理程序。如果所有异常都被屏蔽,则将累积的异常标志OR到浮点状态寄存器中,使用默认结果正常提交。

这种延迟到提交阶段的异常处理方式保证了精确异常(precise exception)的语义——异常发生时,所有在异常指令之前的指令都已完成,所有之后的指令都没有修改体系结构状态。精确异常对于调试和恢复至关重要,但它要求ROB能够存储异常信息,并在提交阶段进行异常决策。

案例研究 5 — 浮点状态寄存器的序列化问题

浮点状态寄存器(如x86的MXCSR)的一个微架构难题是它的序列化效应。由于每条浮点指令都可能修改MXCSR的异常标志位(特别是几乎每条指令都会设置Inexact标志),MXCSR成为了一个"全局共享的可写状态"——所有浮点指令对MXCSR的写入必须按程序序完成。

在朴素的实现中,这意味着MXCSR的异常标志位不能被推测性地修改,必须等到指令提交时才更新。这本身不影响正常的乱序执行(因为很少有指令读取MXCSR的异常标志位),但如果有指令显式读取MXCSR(如x86的STMXCSR指令),它必须等到所有之前的浮点指令都提交后才能执行——这实际上是一个序列化点。

RISC-V的fcsr寄存器面临类似的问题。RISC-V的CSRRS/CSRRW指令对fcsr的读写同样可能引入序列化效应。为了缓解这个问题,某些RISC-V处理器实现(如BOOM)采用了异常标志重命名的技术——将fcsr的异常标志位视为一个可重命名的寄存器,通过OR操作而非覆盖写来更新,从而减少了对前续指令提交的等待。

浮点单元的验证

浮点单元是处理器中最难验证的模块之一。其验证难度来自三个方面:(1)输入空间极大——双精度FMA有三个64位输入,总输入空间为21922^{192},穷举不可能;(2)IEEE 754的正确舍入要求精确到最后一位——不允许"差1个ULP"的近似;(3)特殊值和边界条件众多——NaN、无穷、零、非规格化数的各种组合必须全部正确处理。

验证方法

现代浮点单元的验证通常结合以下方法:

1. 随机仿真测试。 生成大量随机的输入操作数,在RTL仿真中执行运算,然后将结果与高精度软件参考模型(如MPFR库)的结果进行比较。随机测试的覆盖率取决于测试向量的数量和质量——纯随机测试在21922^{192}的输入空间中几乎不可能覆盖边界情况。因此,通常使用约束随机(constrained random)生成——将输入约束在特殊值附近、指数边界附近、尾数全0/全1附近等高风险区域。

2. 定向测试。 手工编写针对特定边界条件的测试用例。关键的定向测试包括:

  • 所有特殊值组合的穷举(5×5=255 \times 5 = 25种对于加法,53=1255^3 = 125种对于FMA)。

  • 最大/最小指数附近的运算(溢出/下溢边界)。

  • 舍入边界——选择使G=1,R=0,S=0G=1, R=0, S=0(tie)的操作数,验证ties-to-even规则。

  • 大规模消去——选择ABA \approx B使得减法后前导零>50> 50位。

  • 非规格化数输入和非规格化结果。

  • 后规格化触发——选择使舍入+1导致尾数溢出的操作数。

3. 形式化验证。 使用数学证明方法来验证浮点运算的正确性——不需要枚举所有输入,而是从逻辑上证明对于所有可能的输入,RTL实现的输出都等于IEEE 754参考模型的输出。形式化验证在浮点领域有成功的应用案例——例如AMD使用ACL2定理证明器验证了其浮点除法和平方根的正确性,Intel在Pentium FDIV Bug之后大幅增加了浮点单元的形式化验证投入。然而,形式化验证的适用范围受限于模型的复杂度——完整的FMA流水线(含数千个寄存器和数万个组合逻辑门)对于当前的形式化工具来说仍然是极大的挑战,通常需要将设计分解为子模块(乘法器、加法器、舍入逻辑等)分别验证。

案例研究 6 — Intel Pentium FDIV Bug的教训

1994年,Intel Pentium处理器中发现了一个浮点除法错误(FDIV Bug):SRT除法器的商位选择查找表中有5个条目的值不正确,导致在某些特定的输入组合下,除法结果的第4位十进制之后出现错误。这个Bug的根本原因是验证不充分——查找表的内容是通过脚本生成的,但脚本中的一个逻辑错误导致5个条目未被正确初始化,而验证过程未能覆盖这些条目对应的输入组合。

Intel为此召回了所有受影响的Pentium处理器,直接成本约4.75亿美元。更深远的影响是推动了整个半导体行业对浮点单元验证的重视——Intel此后建立了专门的浮点验证团队,广泛采用形式化验证方法。这个事件也促进了IBM和AMD等公司加大浮点验证的投入。

从技术角度看,FDIV Bug的教训是:查找表内容的验证至少应该与查找表设计本身同等重要。查找表是浮点除法和平方根实现中最容易出错的部分——因为表的内容是由数学计算生成的,而非由HDL逻辑描述的,传统的RTL仿真可能不会覆盖表中所有条目。对查找表的穷举验证(遍历所有条目并验证其值的正确性)应该是浮点单元验证的必要步骤。

浮点格式转换

浮点格式转换(precision conversion)是处理器浮点单元必须支持的辅助操作,包括不同浮点精度之间的转换(如单精度\leftrightarrow双精度)以及整数与浮点之间的转换(如FCVT.W.SFCVT.D.L等)。虽然格式转换在浮点运算中不如加法和乘法核心,但它的使用频率不低——每次函数调用中的参数传递、混合精度计算中的精度提升/降低、以及数据I/O中的格式适配都涉及格式转换。

浮点精度提升

精度提升(widening conversion),如单精度\to双精度(FCVT.D.S),是一种精确操作——不会引入舍入误差。其硬件实现非常简单:

  1. 符号位:直接复制。

  2. 指数:将8位偏移指数转换为11位偏移指数。设单精度偏移指数为EsE_s(偏移量127),双精度偏移指数为EdE_d(偏移量1023),则Ed=Es127+1023=Es+896E_d = E_s - 127 + 1023 = E_s + 896。这是一个简单的常数加法。

  3. 尾数:将23位尾数字段扩展到52位,低29位补零。

  4. 特殊值处理:零\to零、\infty \to \infty、NaN\toNaN(需要将尾数的NaN有效载荷扩展到更宽的格式)。

  5. 非规格化数:单精度非规格化数在双精度中可能是规格化数(因为双精度的指数范围更宽)。需要对非规格化输入执行规格化:LZC + 左移 + 指数修正。

精度提升的延迟通常为1\sim2个周期,主要由非规格化数的规格化流程决定。对于规格化输入,精度提升可以在1个周期内完成。

浮点精度降低

精度降低(narrowing conversion),如双精度\to单精度(FCVT.S.D),需要舍入——因为目标格式的精度低于源格式。其硬件流程为:

  1. 指数范围检查:双精度指数范围(1022-1022+1023+1023)远大于单精度(126-126+127+127)。如果源指数超出目标范围,产生溢出(±\to \pm\infty)或下溢(\to非规格化或零)。

  2. 尾数截断与舍入:将52位尾数截断到23位,使用G/R/S位进行舍入。G位 = 源尾数第24位,R位 = 第25位,S位 = 第26\sim52位的OR归约。

  3. 非规格化结果:当转换后的指数下溢(0\leq 0),需要将尾数右移并将指数固定为0,同时重新计算G/R/S位——与非规格化数输出的处理完全相同。

精度降低的延迟通常为2\sim3个周期(包含舍入和溢出检测),通常与浮点加法器共享舍入逻辑。

整数与浮点之间的转换

整数到浮点的转换(如FCVT.D.W:32位有符号整数\to双精度)需要以下步骤:

  1. 绝对值化:如果是有符号整数且为负数,取补码得到绝对值,记录符号。

  2. 前导零计数:对绝对值执行LZC,确定最高有效位的位置kk。这决定了指数值E=k+biasE = k + \text{bias}

  3. 尾数对齐:将绝对值左移到规格化位置,使最高有效位成为隐含的"1"位。

  4. 舍入:如果整数的有效位数超过目标浮点格式的精度(如64位整数\to双精度,有效位可能超过53位),需要舍入。32位整数\to双精度不需要舍入(因为32位<<53位精度),但64位整数\to双精度需要舍入。

浮点到整数的转换(如FCVT.W.D:双精度\to32位有符号整数)的流程为:

  1. 指数检查:如果浮点数的幅度231\geq 2^{31}(有符号)或232\geq 2^{32}(无符号),产生溢出。IEEE 754-2019定义溢出结果为目标格式的最大/最小值。

  2. 尾数右移:将尾数右移53Eactual53 - E_{\text{actual}}位(双精度),使小数点对齐到整数位置。

  3. 舍入:根据移出的位执行舍入(如RNE模式下向偶数舍入)。

  4. 符号处理:如果目标是有符号整数且源为负数,取补码。

设计提示

整数\leftrightarrow浮点转换通常在浮点执行单元上执行(而非整数ALU),因为它们需要浮点舍入逻辑和特殊值处理。转换指令的延迟通常为2\sim4个周期。在某些处理器中,转换操作与浮点加法器共享硬件——它们都需要移位器和舍入逻辑。但转换操作需要从整数寄存器文件读取操作数或向整数寄存器文件写入结果,这涉及整数和浮点执行域之间的跨域数据传输——通常需要额外的1\sim2个周期的延迟来完成寄存器文件之间的数据搬移。

低精度格式转换的特殊挑战

随着低精度浮点格式(FP16、BF16、FP8)的广泛使用,格式转换在混合精度计算中的频率大幅增加。深度学习训练的典型模式是:以BF16精度执行前向/反向计算,以FP32精度维护主权重(master weight)——每次参数更新都涉及BF16\toFP32(权重提升)和FP32\toBF16(权重降低)的转换。在一个大型语言模型(如GPT-4级别,\sim1万亿参数)的训练中,每个训练步骤涉及数万亿次格式转换操作。

FP8格式的转换更为复杂。FP8有两种变体:E4M3(4位指数 + 3位尾数)和E5M2(5位指数 + 2位尾数),它们的指数位宽和偏移量都与FP16/FP32不同,转换需要指数重编码(re-biasing)和尾数截断/扩展。此外,FP8的极低精度(1\sim3位十进制)使得转换中的舍入误差相对更大,需要配合量化缩放因子(scaling factor)来使用——转换前先将数据乘以缩放因子以充分利用FP8的动态范围,转换后再除以缩放因子恢复。缩放因子的管理增加了软件和硬件的协同复杂度。

对于高吞吐量的格式转换需求,处理器可以在以下层面进行优化:

  • 向量化转换:SIMD指令一次转换多个值。例如x86的VCVTNE2PS2BF16指令将8个FP32值批量转换为8个BF16值。

  • 转换与运算融合:在FMA的输入或输出端集成格式转换逻辑,避免独立的转换指令。例如"BF16输入\to内部FP32运算\toBF16输出"的FMA可以在一条指令内完成两次转换和一次运算。

  • 隐式转换:某些格式对之间的转换可以通过位拼接或截断隐式完成,不需要显式的转换逻辑。BF16\leftrightarrowFP32是最典型的例子——扩展只需低位补零,截断只需丢弃低16位(加舍入)。

低精度浮点格式与AI加速

近年来,深度学习工作负载的爆发式增长推动了低精度浮点格式在处理器和加速器中的广泛应用。传统的FP32/FP64格式对于深度学习推理和某些训练场景来说精度"过剩",而低精度格式可以在较小的精度损失下获得显著的性能和能效优势。

半精度与BFloat16

两种最重要的低精度浮点格式是:

IEEE 754半精度(FP16/binary16):16位,1位符号 + 5位指数 + 10位尾数。提供约3.3位十进制精度和10510^{-5}6.5×1046.5 \times 10^4的动态范围。FP16的精度和动态范围对于某些深度学习推理任务是充足的,但对训练来说动态范围太窄——梯度值经常超出FP16的可表示范围。

BFloat16(BF16):16位,1位符号 + 8位指数 + 7位尾数。BF16是Google为机器学习训练专门设计的格式——它保持了与FP32相同的8位指数(因此动态范围与FP32相同),但将尾数从23位缩减到7位。这意味着BF16可以直接通过截断FP32的低16位来得到(无需指数转换),而且其动态范围足以容纳训练中的梯度值。

表 31.16对比了各种低精度浮点格式。

格式位宽指数尾数动态范围十进制精度
FP32 (binary32)3282310±3810^{\pm38}\sim7.2位
FP16 (binary16)1651010±510^{\pm5}\sim3.3位
BFloat16168710±3810^{\pm38}\sim2.4位
TF321981010±3810^{\pm38}\sim3.3位
FP8 (E4M3)843$\sim
10^{\pm2}$ | $\sim$1.3位 | | FP8 (E5M2) | 8 | 5 | 2 | $\sim

10^{\pm5}$ | \sim0.9位 |

低精度浮点格式的对比

低精度FMA的硬件设计

低精度浮点运算的硬件设计比高精度运算简单得多,但其核心价值在于并行度——在相同的硅面积和数据通路宽度下,可以同时执行更多的低精度运算。

以256位SIMD通路为例:

  • FP64 FMA:256/64=4256 / 64 = 4个并行运算,每周期8 FLOP

  • FP32 FMA:256/32=8256 / 32 = 8个并行运算,每周期16 FLOP

  • FP16/BF16 FMA:256/16=16256 / 16 = 16个并行运算,每周期32 FLOP

低精度FMA的硬件实现有两种策略:

策略一:分区复用。 将宽数据通路的FMA单元通过分区(partitioning)来支持多个并行的低精度运算。例如,一个53×\times53位的双精度乘法器可以分区为两个24×\times24位的单精度乘法器或四个11×\times11位的半精度乘法器。分区通过在乘法器阵列的适当位置插入"分区屏障"(partition barrier)来实现——屏障阻止进位跨越分区边界传播。

策略二:专用低精度单元。 为低精度运算配置独立的专用执行单元。这种方式可以针对低精度进行更激进的优化(如更少的流水级、更低的供电电压),但增加了总体面积。NVIDIA的GPU中的Tensor Core就是这种策略的代表——Tensor Core是专用的低精度矩阵乘法单元,与通用的FP32/FP64 CUDA Core并行存在。

设计权衡 4 — 低精度浮点运算的精度分区 vs 专用单元

  • 精度分区:面积效率高(复用已有硬件),但分区逻辑增加了关键路径延迟;低精度运算的延迟与高精度相同(因为使用同一流水线),无法单独优化。代表:Intel AVX-512的FP16支持。

  • 专用单元:可以针对低精度进行独立优化(更少流水级、更低延迟),但需要额外的面积和布线资源;需要独立的发射端口和旁路网络。代表:NVIDIA Tensor Core、Intel AMX。

  • 混合方案:在通用FMA单元上通过分区支持FP16/BF16运算,同时配置专用的矩阵单元(如Intel AMX)用于大规模低精度矩阵运算。代表:Intel Sapphire Rapids及后续设计。

BFloat16的硬件简化

BFloat16的一个重要硬件优势是其与FP32的格式兼容性。由于BF16的指数位数与FP32相同(都是8位,偏移量127),BF16\toFP32的转换只需要将7位尾数扩展到23位(低16位补零),不需要任何指数转换——这是一个纯粹的位拼接操作,不需要任何逻辑门。反方向的FP32\toBF16转换也很简单:截断低16位尾数,只需处理舍入(将第16位及以下的位作为G/R/S位)。

这种格式兼容性对硬件设计有以下好处:

  • 混合精度累加:BF16 FMA可以使用BF16乘法器(8×\times8位尾数乘法)但使用FP32的累加器。乘积的16位尾数可以直接扩展到FP32精度后累加,不需要指数转换。这种"BF16乘法 + FP32累加"模式是深度学习训练中最常见的混合精度策略。

  • 数据通路复用:BF16操作数可以直接存储在FP32寄存器的高16位中,低16位为零。这允许BF16运算直接使用FP32的寄存器文件和数据通路,不需要独立的BF16寄存器文件。

性能分析 11 — 低精度FMA的能效优势

低精度浮点运算不仅提供了更高的吞吐量,还具有显著的能效优势。乘法器的能耗大致与位宽的平方成正比(因为乘法器面积n2\propto n^2),加法器的能耗与位宽成线性关系。对于FMA运算(乘法占主导),各精度格式的相对能耗估算如下:

格式尾数乘法位宽相对乘法能耗每焦耳FLOP
FP6453×5353 \times 53\sim100%1×\times
FP3224×2424 \times 24\sim20%5×\times
FP1611×1111 \times 11\sim4.3%23×\times
BF168×88 \times 8\sim2.3%43×\times
FP8 (E4M3)4×44 \times 4\sim0.6%170×\times

这些数字是近似的(实际能耗还取决于寄存器文件读写、控制逻辑、时钟树等),但趋势是明确的:精度每降低一半,FMA的能效提升约4\sim5倍。这是深度学习加速器追求低精度运算的根本驱动力——在给定的功耗预算下,BF16 FMA可以提供约40倍于FP64的算力密度。

案例研究 7 — 从FP32到BF16:深度学习训练的精度演进

深度学习训练中精度格式的演进反映了精度与性能的持续权衡:

  • 2012\sim2017:FP32是训练的标准精度。AlexNet到ResNet都使用FP32训练。

  • 2017\sim2019:混合精度训练(FP16计算 + FP32主权重)成为主流。NVIDIA V100的FP16 Tensor Core是关键推动力。但FP16的窄动态范围需要损失缩放(loss scaling)技术来防止梯度下溢。

  • 2019\sim至今:BFloat16逐渐取代FP16用于训练。BF16的FP32级动态范围消除了损失缩放的需求,简化了训练流程。Google TPU从v2开始原生支持BF16,Intel从Cooper Lake开始在CPU上支持BF16(AVX-512 BF16扩展),AMD从Zen 4开始支持。

  • 2023\sim至今:FP8格式开始用于推理和部分训练场景。NVIDIA H100和AMD MI300支持FP8运算。FP8的极低精度(1\sim3位十进制)要求精心的量化策略和校准。

从处理器设计的角度看,这个趋势意味着浮点单元需要支持越来越多的精度格式,每种格式都需要对应的运算逻辑、舍入逻辑和格式转换逻辑。这对浮点单元的设计复杂度和验证覆盖率提出了巨大挑战——每增加一种精度格式,验证矩阵就扩展一个维度。

浮点流水线的整体集成

前面各节分别讨论了浮点加法器、乘法器、FMA、除法/平方根以及异常处理的设计。本节将这些部件放回处理器微架构的全局视角,讨论浮点执行域的整体组织。

浮点执行端口的配置

现代超标量处理器的浮点执行资源通常组织为多个执行端口(execution port),每个端口连接一条或多条功能单元流水线。图图 31.5展示了一个典型的浮点执行域配置。

典型的浮点执行域配置:4个执行端口分别连接FMA、FADD、FDIV/FSQRT和FP Store等功能单元。
典型的浮点执行域配置:4个执行端口分别连接FMA、FADD、FDIV/FSQRT和FP Store等功能单元。

不同处理器的浮点执行端口配置各有特色:

  • Intel Golden Cove:Port 0(FMA + FADD + FMUL + FDIV/FSQRT),Port 1(FMA + FADD + FMUL)。两个FMA端口共享FDIV/FSQRT的单一实例。

  • AMD Zen 4:Port 0/1各有一个FMA流水线,Port 2/3各有一个FADD流水线。FDIV/FSQRT挂在Port 0上。AMD的配置比Intel多了独立的FADD端口,对加法密集型工作负载更有利。

  • Apple Everest:据推测有4\sim6个FP/SIMD执行端口,包含4个FMA流水线,是目前最激进的浮点执行配置。

  • 香山昆明湖:2个FMA端口 + 1个FADD端口 + 1个FDIV/FSQRT端口。

浮点旁路网络

浮点执行域的旁路网络(bypass/forwarding network)是确保流水线高效运行的关键。当一条浮点指令的结果需要被后续指令立即使用时,旁路网络将结果从产生它的功能单元直接转发到消费它的功能单元的输入端,避免了经过寄存器文件的写入-读取延迟(通常1个周期)。

浮点旁路网络的设计挑战在于:

  • 不同延迟的功能单元:FADD(3周期)、FMA(4周期)、FMUL(4周期)的延迟不同,旁路网络需要在正确的时刻从正确的流水级提取结果。

  • 宽数据通路:浮点旁路的数据宽度等于最宽的浮点数据格式(如512位的AVX-512),旁路多路选择器的面积和延迟都不可忽视。

  • 跨端口旁路:从一个端口的功能单元旁路到另一个端口的功能单元,需要跨越物理布局上的距离,布线延迟可能显著。

一个重要的微架构决策是旁路时刻(bypass timing)。如果FMA的结果在第4个周期(流水线最后一级)可以旁路,那么依赖FMA结果的后续指令可以在第5个周期读取旁路数据——即"1周期旁路延迟"。某些设计选择在第3个周期(FMA的倒数第二级,舍入之前)提供"非精确旁路"(speculative bypass),允许后续指令提前开始执行,但如果舍入改变了结果的低位,需要重新执行——这是一种以正确性复杂度换取性能的激进优化,在工业设计中较少使用。

浮点寄存器文件的组织

浮点物理寄存器文件(FP PRF)的组织与整数PRF有以下不同:

  • 条目宽度:浮点PRF的每个条目宽度通常等于最宽的SIMD寄存器宽度。例如支持AVX-512的处理器中,每个FP PRF条目为512位。对于RISC-V V扩展,条目宽度等于VLEN。

  • 读写端口:浮点PRF需要为FMA指令提供3个读端口(三个源操作数aabbcc)和1个写端口。如果有NN个发射端口,总共需要3N3N个读端口和NN个写端口。这使得FP PRF的读端口数量通常多于INT PRF(整数指令通常只有2个源操作数)。

  • 条目数量:FP PRF的条目数量需要覆盖所有浮点指令的飞行窗口(in-flight window)。典型的FP PRF有128\sim256个条目——少于INT PRF(通常256\sim384个条目),因为浮点指令的飞行窗口通常小于整数指令。

设计提示

浮点PRF的面积开销在SIMD向量处理器中尤为显著。一个支持AVX-512的浮点PRF,如果有256个条目、每个条目512位、10个读端口和4个写端口,其总容量为256×512=128Kbit=16KB256 \times 512 = 128\,\text{Kbit} = 16\,\text{KB}——接近L1数据缓存的大小。在先进工艺下,这样的多端口寄存器文件的面积约为0.51mm20.5\sim 1\,\text{mm}^2,是浮点域中面积最大的单一结构。

减小FP PRF面积的技术包括:(1)PRF分库(banking)——将PRF分为多个库,每个库的端口数量减少,通过仲裁逻辑在发射时分配库端口;(2)物理寄存器共享——当多条指令的目标寄存器值相同(如连续写入同一寄存器),共享同一个物理寄存器条目;(3)窄值压缩——当SIMD寄存器中只有部分通道被使用时(如256位AVX操作在512位寄存器中),只分配所需的存储空间。

整数域与浮点域的跨域数据传输

在现代处理器中,整数执行域和浮点执行域通常拥有独立的物理寄存器文件独立的发射队列。当指令需要在两个域之间传输数据时(如浮点比较结果写入整数寄存器、整数到浮点的格式转换、或浮点位模式的整数重解释),就涉及跨域数据传输(cross-domain transfer)。

跨域传输的微架构实现有几种方式:

方式一:通过加载/存储缓冲区中转。 最简单但最慢的方式:先将数据从一个域的寄存器存储到内存(Store),再从内存加载到另一个域的寄存器(Load)。这种"存储-加载"方式的延迟约为5\sim10个周期(存储延迟 + 存储转发延迟),完全不可接受用于频繁的跨域操作。

方式二:直接寄存器间传输。 处理器提供专用的跨域传输通路,允许数据直接从一个PRF读出并写入另一个PRF。传输延迟通常为1\sim3个周期。x86的MOVD/MOVQ指令和RISC-V的FMV.W.X/FMV.X.W指令使用这种方式。这需要在INT PRF和FP PRF之间增加物理布线——由于两个PRF在芯片上通常不相邻,布线延迟可能较长。

方式三:统一寄存器文件。 某些处理器使用统一的物理寄存器文件(Unified PRF),整数和浮点寄存器共享同一个物理存储结构。在这种设计中,跨域传输只是一个寄存器重命名操作(将浮点逻辑寄存器映射到整数逻辑寄存器已指向的物理寄存器),不需要数据搬移,延迟为零。但统一PRF的缺点是条目宽度必须取整数和浮点中较宽的那个,浪费了窄数据类型的存储空间,且端口数量需要同时满足整数和浮点的需求。

性能分析 12 — 跨域传输延迟的性能影响

跨域数据传输在以下场景中频繁出现:

  • 浮点比较 + 条件分支:RISC-V的FEQ.S rd, fs1, fs2将浮点比较结果写入整数寄存器rd,随后的BNE rd, zero, target需要从整数寄存器读取rd。如果跨域传输延迟为2周期,则分支解析被额外延迟2周期——这对分支预测器的训练和流水线效率有直接影响。

  • 混合整数/浮点计算:在某些算法中(如视频编解码),需要频繁地将浮点计算结果转换为整数格式。每次FCVT.W.S\to整数运算的转换链中,跨域传输引入的延迟可能成为瓶颈。

  • 位操作技巧:如前述的"快速逆平方根"和某些SIMD掩码操作,需要将浮点位模式重解释为整数进行位操作。

在Intel的微架构中(Skylake\simGolden Cove),整数到浮点的传输延迟约为1\sim2周期,浮点到整数约为2\sim3周期(更慢是因为需要等待浮点流水线完成)。AMD Zen系列的跨域延迟类似。Apple的处理器据测量具有更低的跨域延迟(约1周期),可能使用了统一或部分统一的寄存器文件设计。

浮点单元的功耗管理

浮点单元是处理器中功耗最高的执行资源之一。双精度FMA单元(含乘法器阵列、CSA树、CPA、移位器和舍入逻辑)的动态功耗在先进工艺下约为50100mW50\sim 100\,\text{mW}(在典型频率和活动率下),两个FMA单元加上寄存器文件和旁路网络,浮点域的总功耗可达200400mW200\sim 400\,\text{mW}——在一个总功耗约510W5\sim 10\,\text{W}的处理器核心中占到了显著比例。

浮点单元的功耗管理技术包括:

  • 时钟门控(Clock Gating):当没有浮点指令在执行时,关闭浮点单元的时钟信号以消除动态功耗。时钟门控是最基本也最有效的功耗管理技术——在纯整数工作负载中,浮点单元的时钟可以完全关闭。时钟门控的粒度可以从整个浮点域(粗粒度)到单个流水级(细粒度)不等。

  • 操作数门控(Operand Gating):在浮点单元的输入端,当操作数为零或特殊值时,门控后续的数据通路以减少无效的信号翻转。例如,当FMA的c=0c = 0时(退化为纯乘法),加数注入路径可以被门控。

  • 精度自适应功耗优化:当执行单精度运算时,双精度数据通路的高32位可以被门控,减少约一半的动态功耗。类似地,执行半精度运算时可以门控更多的高位。这种精度感知功耗优化在混合精度工作负载中(如深度学习推理中FP16/BF16运算占主导)尤为有效。

  • 向量宽度自适应:当执行128位SIMD运算(而非256位或512位)时,高位通道的功能单元可以被门控。Intel的AVX-512设计中,128位SSE操作只激活最低的128位通道,高384位通道被门控。

设计提示

Intel的AVX-512浮点单元在运行512位向量指令时的功耗极高(两个512位FMA单元同时全速运行),以至于处理器在检测到AVX-512指令时会主动降频(频率降低100\sim200 )以控制总功耗在TDP之内。这种频率节流(frequency throttling)对于延迟敏感的工作负载可能产生负面影响——即使AVX-512指令本身的吞吐量很高,降频后的标量/128位代码性能也会下降。AMD Zen 4采用了不同的策略:将AVX-512指令拆分为两个256位μ\muop,避免了同时激活全宽512位数据通路,因此不需要降频(但吞吐量减半)。这两种策略体现了AVX-512性能与功耗的不同权衡点。

本章小结

本章详细探讨了浮点功能单元的微架构设计,从IEEE 754标准的基础知识到具体的硬件实现。核心要点包括:

  1. IEEE 754标准定义了浮点数的表示(符号+指数+尾数)、四种舍入模式(RNE/RZ/RP/RN)和特殊值(NaN、±\pm\infty±0\pm 0、非规格化数)。偏移量选择2k112^{k-1}-1(而非2k12^{k-1})最大化了动态范围上界,隐含位"1"免费获得了一位额外精度。Guard/Round/Sticky三个精度位是实现正确舍入的关键硬件支撑。这些定义直接决定了浮点硬件的设计约束。

  2. 浮点加法器是最复杂的基本浮点运算单元。大规模对阶与大规模消去的互斥性是双路径设计的数学基础——Near Path(ΔE1|\Delta E| \leq 1,需要LZA和大规模规格化)和Far Path(ΔE>1|\Delta E| > 1,需要大规模对阶)分别处理两种互斥场景,显著缩短了关键路径。LZA通过位级逻辑预测前导零位置(可能有±1\pm 1位误差),与尾数减法并行执行,是Near Path性能的关键。并行舍入技术同时计算截断和增量两个版本,消除了舍入在关键路径上的延迟。现代双精度浮点加法器的典型延迟为3\sim4个周期。

  3. 浮点乘法器的核心是尾数乘法,可以与整数乘法器共享Booth编码器和Wallace树硬件。规格化最多需要1位移位,比加法简单得多。早期截断和混合精度分区乘法器是面积优化的重要手段。

  4. FMA(融合乘加)是现代浮点单元的核心运算单元。单次舍入消除了乘法中间结果的舍入误差项abϵ|ab|\cdot\epsilon,在灾难性消去场景下优势尤为显著——分开计算可能导致100%的相对误差,而FMA始终保持ϵ\epsilon级相对精度。其微架构设计的关键在于将加数cc注入乘法器的CSA压缩树中实现乘加融合,内部数据通路宽度达160\sim170位。典型延迟为4\sim5个周期,吞吐量可达每周期2次。

  5. 浮点除法和平方根使用迭代算法(Newton-Raphson或Goldschmidt),复用乘法器或FMA硬件。Goldschmidt方法的两个独立乘法可以在双FMA端口上并行执行,每次迭代延迟仅为1×TFMA1\times T_{\text{FMA}}(Newton-Raphson为2×TFMA2\times T_{\text{FMA}})。SRT与乘法迭代的混合方法是工业界的主流选择。除法和平方根通常是非全流水线的,延迟为13\sim21个周期。

  6. 浮点异常处理涉及五种IEEE 754异常的检测和报告。异常标志通过流水线逐级传递,在提交阶段按序写入状态寄存器以保证精确异常语义。非规格化数的处理是最大挑战:输入需要额外的前导零计数和预规格化左移,输出需要下溢检测、右移和避免双重舍入的延迟舍入策略。完全硬件支持、微码辅助和FTZ/DAZ模式代表了三种不同的设计权衡,现代处理器正逐步向完全硬件支持演进。

  7. 低精度浮点格式(FP16、BFloat16、FP8)在AI工作负载的驱动下日益重要。BF16的8位指数与FP32兼容,使得格式转换近乎零开销。低精度FMA的实现可以通过分区复用(partition)或专用单元两种策略。浮点执行域的整体集成——包括端口配置、旁路网络和宽SIMD寄存器文件——决定了浮点单元的实际性能上限。

设计权衡 5 — 前向桥接——从功能单元到存储系统

至此,我们完成了乱序执行引擎的核心数据通路:寄存器重命名消除假依赖(第 24.0 章\sim第 26.0 章),发射队列缓冲ILP(第 27.0 章\sim第 29.0 章),功能单元执行运算(本章和第 30.0 章)。然而,一个关键的缺失环节是存储器访问——Load和Store指令是处理器中最复杂的指令类型。Load的延迟不确定(L1命中4周期、L2命中12周期、主存200+周期),这正是发射队列投机唤醒(第 28.0 章)面临的核心挑战。Store的语义要求按程序序可见,但乱序执行引擎的本质就是打破程序序——如何调和这一矛盾?Store Buffer、Store-to-Load转发和内存消歧将在第 32.0 章中展开讨论。

正文与图片:CC BY-NC-SA 4.0 · 本仓库少量站点配置代码:MIT