下载掌阅APP，畅读海量书库
立即打开
尽管大多数现代CPU都有支持浮点运算的FPU硬件，但是为了理解浮点运算背后的原理，花些时间编写浮点算术运算程序是值得的。一般会选择使用更快的汇编语言编写数学函数，因为浮点包的主要设计目标就是速度快。但这里我们实现浮点包只是想看清楚算法流程，因此选择的是更容易编写、阅读和理解的代码。
高级编程语言（如C/C++或Pascal）实现浮点数的加法和减法运算确实很容易，所以我们用高级语言来实现浮点数的加法和减法运算。而浮点数的乘法和除法运算，使用汇编语言实现则更容易，所以我们使用高级汇编语言（HLA）实现乘法和除法程序。
本节我们使用IEEE32位单精度浮点格式（如图4-2所示），带符号的值使用1的补码表示。也就是说，符号位（第31位）为1则数字为负，为0则数字为正。第23～30位是一个8位增码-127格式的阶码，尾数占用余下的24位。由于高位隐含为1，因此这种格式不支持非规约化的值。
加法和减法的实现代码基本上是一样的。毕竟，X-Y就等于X+（-Y）。如果可以把负数加到另一个数上，那么要实现两个数的减法，就可以先取一个数的负值，然后加到另一个数上。而且IEEE浮点格式使用1的补码表示负数值，取负非常简单，只需要把符号位取反即可。
因为使用的是标准的IEEE32位单精度浮点格式，所以理论上我们可以直接使用C/C++的浮点数据类型（假设底层的C/C++编译器也使用这种格式，实际上大多数现代计算机也是这样做的）。但很快我们就会发现，用软件实现浮点计算需要将浮点格式当成位串或者整数值，才能进行各种字段操作。因此，把浮点表示形式的位模式当成32位无符号整型，操作起来更方便。假设C/C++的无符号长整型的实现是32位（本节假设这种长整型就是uint32_t，相当于typedefunsignedlonguint32_t），为了区分实数值和程序中它们的实际整型表示，我们定义了下面这个real数据类型，所有实数变量都用这种类型声明：
采用和C/C++一模一样的浮点格式实现有一个好处：可以将浮点字面值直接赋值给real变量，而且其他浮点操作也可以使用现有的库，比如输入和输出。这种做法也有问题：如果real变量出现在浮点表达式中，C/C++会自动将整型转换成浮点格式（记住，对C/C++来说real只是一个无符号长整型值）。这意味着我们需要让编译器把real变量中保存的位当作float对象来处理。
直接采用 （float）realVariable 强制转换类型没有效果。C/C++编译器会认为 realVariable 中包含的是整数值，编译代码时会将其转换为和这个整数值相等的浮点值。而我们期望C/C++编译器把 realVariable 的位当作浮点表示形式，不做任何转换。下面这个C/C++宏可以巧妙地做到这一点：
这个宏的参数必须是一个real变量，而编译器会认为其结果就是float变量。有了float变量，我们就可以开发两个C/C++函数fpadd（）和fpsub（），分别完成浮点数的加法和减法运算。这两个函数都接受三个参数，分别是运算符的左操作数、右操作数和一个指向结果的指针。两个函数原型如下：
fpsub（）函数对右操作数取负，再用取负的结果调用fpadd（）函数。fpsub（）函数的代码如下：
所有的实际运算都是由fpadd（）函数完成的。我们把fpadd（）函数分解成几个不同的工具函数，这样理解和维护起来更容易。真正的软件浮点库程序是不会这样做的，因为额外的子程序调用会让计算更慢。这里开发fpadd（）函数只是为了学习。再说，实现高性能的浮点数加法运算，使用硬件FPU比使用软件实现更合适。
IEEE浮点格式是典型的打包数据类型。第3章介绍过，打包数据类型能够显著降低一个数据类型要求的存储空间，但是在真正计算时使用打包的字段就不那么方便了。因此，浮点函数要做的第一件事就是把符号、阶码和尾数三个字段从浮点表示形式中提取出来。
第一个函数是extractSign（），负责从打包的浮点表示形式中提取符号位（第31位），返回0（代表正数）或1（代表负数）。
下面这个表达式也可以提取符号位（而且更高效）：
但是，在前面一段代码中，将第31位右移到第0位显然更容易理解。
第二个工具函数是extractExponent（），它提取的是打包实数格式中第23～30位的阶码。具体实现如下：先将实数值向右移动23位，再掩去符号位，最后把得到的增码-127格式的阶码转换为2的补码格式（减去127即可）。
然后是从实数值中提取尾数的extractMantissa（）函数。提取尾数必须屏蔽掉浮点表示形式中的阶码位和符号位，再将隐含为1的高位插入尾数中。只有全部位都为0这一种情况例外，这时必须返回0。
前面讲过，对采用科学记数法表示（IEEE浮点格式就是）的两个值进行加减法运算的时候，必须首先对齐这两个值的阶码。例如，两个十进制数（基数为10）1.2345e3和8.7654e1相加，必须首先调整其中一个数，使两个数的阶码相等。可以右移第一个数的小数点，减少它的阶码。例如，下面这些值都等于1.2345e3：
同样，可以左移小数点让阶码增加。以下这些值都等于8.7654e1：
对于二进制数的浮点加减法运算，尾数左移一位同时递减指数，或者尾数右移一位同时递增指数，都可以把二进制阶码调整到一样。
尾数向右移位会降低数字的精度（因为最终尾数的最低位会被移出去）。为了尽可能地保持计算的准确性，从尾数中移出去的位不应当被直接截断，而是应该将结果舍入到余下的尾数位能够表示的最接近的值。下面按顺序列出了IEEE的舍入规则：
1.如果最后移出的一位是0，则截断结果。
2.如果最后移出的一位是1，且其他移出的位中至少有一位为1，则最终的尾数加 1。
3.如果最后移出的一位是1，且其他移出的位全部为0，尾数的低位为1，则最终的尾数进1。
尾数的移位和舍入是一种相对复杂的操作，而且浮点加法运算的代码会执行多次。因此，可以将此运算实现为一个工具函数。下面列出的代码是在C/C++中实现的shiftAndRound（）函数：
这段代码“巧妙”地利用masks和HOmasks两张查询表就提取到了右移操作会用到的尾数位。masks表中每一项掩码为1的位（置位的位）都是移位时会被移出去的位。HOmasks表中的每一项掩码只有其索引对应的位是1，换句话说，索引0的掩码第0位为1，索引1的掩码第1位为1，依此类推。上面这段代码根据尾数需要右移的位数（索引）来查询两张表中对应的掩码。
将尾数的初始值和masks表中对应的掩码进行逻辑AND运算，如果得到的结果大于HOmasks表中的对应项，那么shiftAndRound（）函数会将移位后的尾数向上舍入为较大的值。如果得到的结果等于HOmasks表中对应的项，那么实现代码会根据移位后尾数值的低位进行舍入（注意，当尾数低位为1时，表达式（*valToShift ＆ 1）的结果为1，反之为0）。最后，如果得到的结果小于HOmasks表中对应的项，则实现代码什么都不用做，因为尾数已经舍入过了。
当通过调整一个操作数使得两个操作数的阶码对齐之后，加法算法下一步要做的就是比较两个值的符号。如果两个操作数的符号相同，则把尾数直接相加（使用标准的整数加法操作）。如果符号不同，则将尾数相减，而不是相加。因为浮点表示形式使用的是1的补码，而标准整数算术运算使用的是2的补码，所以无法直接用正数减去负数。必须用较大的值减去较小的值，再根据初始操作数的符号和量级确定最终结果的符号，具体方法如表4-3所示。
表4-3 两个操作数符号不同时的处理
无论进行加法还是减法运算，两个24位数字的运算结果都可能是25位（实际上这种情况在处理规约化数值时很常见）。浮点运算代码在加减法完成之后必须立即检查结果是否发生了溢出。如果发生了溢出，则尾数需要右移1位并对结果进行舍入，而且阶码也要递增。完成这一步之后，剩下的工作就是将得到的符号、阶码和尾数三个字段打包成32位的IEEE 浮点格式。负责这项工作的packFP（）函数代码如下：
注意，这个函数适用于规约化数值、非规约化数值和零，但不适用于NaN和无穷。
使用这些工具函数，fpadd（）就可以完成两个浮点数的加法运算，得到32位的实数结果：
最后我们用C中的main（）函数对fpadd（）和fsub（）的软件实现进行总结，这个main（）函数演示了如何使用这两个函数：
使用Microsoft Visual C++编译全部代码（并且把uint32_t定义为unsigned long）得到的输出如下：
大多数软件浮点库实际上是用汇编语言而不是高级（编程）语言编写的，而且经过人工优化。从前面章节的介绍我们知道，用高级语言编写浮点程序是可行的，特别是单精度浮点数的加减运算，实现代码很高效。如果采用合适的库程序，还可以使用高级语言编写浮点数的乘法和除法运算程序。但实际上浮点数的乘除法运算使用汇编语言实现更容易，因此本节介绍用HLA实现的单精度浮点乘法和除法算法。
本节使用HLA实现的函数有两个：fpmul（）和fpdiv（），它们的原型如下：
本节的代码和前一节的代码相比，除了编程语言不同（本节使用汇编语言而不是C语言），还有两处主要的差异。第一，本节的代码并没有为实数值创建新的数据类型，而是使用了内置的real32数据类型。因为在汇编语言中，任何32位的内存对象都可以直接转换为real32或dword类型。第二，这些原型只有两个参数，并没有包含结果指针参数。这些函数直接把real32类型的结果放在了EAX寄存器中。
1.浮点数乘法
当两个用科学记数法表示的数值相乘时，最终结果的符号、阶码和尾数的计算方法如下：
●结果的符号是两个操作数的符号的异或。也就是说，如果两个操作数符号相同，则结果为正；如果两个操作数符号不同，则结果为负。
●结果的阶码是两个操作数的阶码之和。
●结果的尾数是两个操作数的尾数的整数（定点数）乘积。
此外，IEEE浮点格式还有一些额外的规则会影响浮点数乘法，由它们直接得到乘法结果：
●两个操作数有一个或者两个都为0，则结果为0（因为0的表示形式特殊，所以需要特殊处理）。
●任意一个操作数为无穷，则结果为无穷。
●任意一个操作数为NaN，则结果也是NaN。
fpmul（）函数首先检查两个操作数是否为0。如果是，则函数立即向调用者返回0.0。接下来fpmul（）函数检查左操作数和右操作数是否为NaN或无穷。只要操作数中有这类值，就把其直接当作结果返回给调用者。
如果fpmul（）的两个操作数都是有效的浮点值，那么fpmul（）需要提取出打包浮点值中的符号、阶码和尾数字段。实际上，这里说提取并不准确，更准确的说法是分离（Isolate）。下面这段代码分离了两个操作数的符号位，计算出结果的符号：
这段代码先对两个操作数进行异或运算，然后屏蔽掉EBX寄存器的第0～30位，只留下第31位作为结果的符号位。fpmul（）并没有像对一般解包数据那样把符号位移到第0位，因为将计算结果重新打包成浮点值的时候还需要把符号位移回第31位。
处理阶码，fpmul（）需要把第23～30位分离出来进行运算。当将使用科学记数法表示的两个数值相乘时，阶码的值必须相加。但是，阶码之和还要减去127，因为增码-127格式的阶码相加相当于加了两次偏移。下面的代码分离出阶码位，调整了多余的偏移，再将阶码相加：
首先，在这段代码中，减去的是126而不是127。原因是稍后尾数相乘的结果需要乘以2。减去126相当于阶码减去127并将尾数乘以2（这省了一条指令）。
如果上面代码中add（eax，ecx）算出的阶码之和太大，超过了8位可以表示的范围，那么ECX会从第30位进位到第31位，这样会设置80x86的溢出标志。如果在乘法运算中发生了溢出，代码将直接返回无穷作为结果。
如果没有发生溢出，那么fpmul（）需要设置两个尾数隐含的高位。下面这段代码完成了这个任务，去掉了尾数中的所有阶码和符号位，并将EAX和EDX中的尾数左移，对齐到第31位。
EAX和EDX中的尾数都移到第31位之后，下面我们就可以使用80x86的mul（）指令完成乘法运算了：
这条指令计算的是EAX和EDX的64位乘积，结果将被保存到EDX:EAX（高位的双字在EDX中，低位的双字在EAX中）中。因为任意两个 n 位整数的乘积可能达到2× n 位，所以mul（）指令计算的结果是EDX:EAX=EAX×EDX。在开始乘法运算之前，EAX和EDX中的尾数需要左对齐，确保乘积的尾数最后保存在EDX中的第7～30位。实际上需要让这些位出现在EDX的第8～31位。这就是为什么前面在调整增码-127格式的值时，减去的是126而不是127的原因（相当于将结果乘以2，和结果左移一位效果一样）。在进行乘法操作前这些数字是规约化的，所以在乘法运算之后除非结果是0，否则EDX中的第30位是1。32位的IEEE实数格式不支持非规范化值，因此在使用32位浮点数时不必担心这种情况。
因为每个尾数都有24位，所以尾数乘积的有效位可能达到48位。而我们的结果尾数最多只有24位，所以需要舍入得到24位的结果（采用IEEE舍入算法，参见4.4 节）。下面这段代码可以将EDX中的值舍入到24个有效位（第8～31位）：
舍入之后的数字可能需要重新规约化。如果尾数所有位全都为1并且还需要向上舍入，则尾数的最高位会发生溢出。如果发生溢出，可上面代码片段最后的rcr（）和inc（）指令会将溢出的位移回到尾数中。
还剩下最后一件事：将最终的符号、阶码和尾数打包到32位EAX寄存器中。实现代码如下：
在这段代码中，唯一有技巧的地方是lea（）（加载有效地址）指令的使用，只用了这一条指令就完成了EDX（尾数）和ECX（阶码）的加法计算，并将结果移动到了EAX中。
2.浮点数除法
浮点数除法运算比乘法运算更复杂一些，因为IEEE浮点标准花了很大篇幅定义了除法运算过程中可能出现的退化情况。我们不打算在这里讨论处理这些情况的全部代码。如果需要，请参考前面fpmul（）关于这些情况的讨论，或者参考本节即将给出的fdiv（）完整代码清单。
假设需要对两个有效的数值进行除法运算，首先需要使用和乘法一样的算法（以及代码）计算结果的符号。当两个使用科学记数法表示的值相除时，阶码必须相减。与乘法的算法相比，除法的阶码处理更方便一些：只需要解包两个除法操作数中的阶码并将它们从增码-127转换为2的补码形式。实现代码如下：
80x86 div（）指令要求商一定不能超过32位。如果不能满足这个条件，CPU可能会中止操作，抛出除法异常。只要除数的高位为1，被除数的高两位是%01，就不会出现除法错误。下面这段代码完成了除法操作之前的操作数准备：
下一步就是真正的除法运算了。如前所述，为了防止除法错误，我们必须将被除数右移1位（将高两位置为%01），实现代码如下所示：
div（）指令执行完之后，商就保存在EAX中的高24位，余数则保持在AL:EDX中。下面要做的是结果的规约化和舍入。舍入相对来说简单一些，因为除法运算的余数就保存在AL:EDX中。当余数小于$80:0000_0000（即80x86的AL寄存器为$80，而EDX为0）时向下舍入，当余数大于$80:0000_0000时向上舍入。如果余数恰好等于$80:0000_0000，则舍入到最接近的值。
实现代码如下：
fpdiv（）最后的操作是将偏移加回到阶码中（并检查溢出），然后将商的符号、阶码和尾数字段打包为32位浮点格式。实现代码如下：
代码确实不少。但是，完整地实现一遍浮点运算有助于我们理解其背后的原理。大家现在能够明白FPU到底都替我们做了哪些工作了吧。