软件模拟实现IEEE-754单精度浮点数运算
大多数CPU都有硬件的浮点单元（FPU），但是有一些MCU使用的内核（比如Cortex-M3）没有FPU，或者一些内核只支持单精度，同时大部分CPU都不支持高精度128位的浮点数，如果需要使用这些，则需要软件进行模拟（即使用CPU的定点运算指令），使用软件模拟而不用硬件单元的浮点数通常叫做软浮点数。除此之外，软浮点还有一个作用就是平台无关性，因为对于浮点数，不同的硬件实现多少是有小区别的。
本文就以IEEE-754单精度浮点数为例，实现软浮点数的一些基本运算操作。我最早是为了FPC的微控制器开发实现一个精简的单精度软浮点模块，后来遇到了不少坑，花了不少时间才完全实现，所以就记录一下，考虑到大多数人不熟悉Pascal，所以用C重写了一遍。本文的代码不是最优化的，主要是为了一定的代码可读性，不过尽量去优化效率和代码体积。
其实网上可以找到一个非常好的软浮点实现，http://www.jhauser.us/arithmetic/SoftFloat.html，但是其代码比较复杂，很多非常高水平的操作，我也看不懂很多东西。本文的代码有部分内容（主要是比较运算）参考了其实现。
本文代码实现了两种舍入模式，即向0舍入（截断）和向最近的偶数舍入这两种。默认使用向偶数舍入，通过在include头文件之前定义宏#define ROUND_TRUNC可以改成向0舍入，这样可以降低运算精度减小代码体积（在需要的时候，同时需要注意编译器优化级别，因为可能由于 auto inline 优化导致代码体积膨胀）。
考虑到非规格化数并不常用，且会增加代码复杂度，本文代码忽略了非规格化数的操作，在运算中一律视为0。
完整代码在 https://gitee.com/peazomboss/softfloat32 可以找到，使用MIT许可证。
注意：代码仅使用gcc编译测试，不确定其他C或C++编译器。
浮点数简介
在IEEE-754标准之前，各厂商有自己的浮点格式，但是现在IEEE-754已经一统江湖了，目前最新的标准是IEEE-754 2008，不过这个标准的完整文件是要收钱的，很贵哦，当然咱也没必要看，因为网上可以找到大量资料。
简单来说，目前最常用的是单精度浮点数和双精度浮点数，在AI领域还有NV的半精度浮点数，以及一些科学计算需要的128位浮点数。当然还有80位扩展精度浮点数之类的，但是现在基本不常见了。
以32位单精度浮点数为例，其由1位符号（Sign），8位指数（Exponent）和23位有效数（Significand）组成。其中指数也有叫阶码的（本文叫指数），有效数也有叫尾数的（本文不作区分）。
其中8位指数是用偏移量（Bias）表示的，所以指数127相当于实际指数0，指数取值为1~254，其中0和255有特殊用途；23位有效数是规格化（Normalized，也有叫标准化，规约化的）的结果，并不包含隐含的1，加上之后实际上其有效数有24位。
其他的我也不多说了，网上讲的好的非常多，我不干这些重复的事了。
简单说例如数字52，可以表示为1.101×2^5（即1.625*32），那么其符号为0，指数为5+127=132，有效数为1.101，规格化后的结果为0 10000100 1010000000...，基本就是这样。对于其他精度的浮点数，也就是指数和有效数的不同。
基本代码
头文件softfloat32.h基本内容如下：
#include <stdint.h>
#include <stdbool.h>
typedef uint32_t sfloat32;
// 将浮点数转换到软浮点数
sfloat32 sf32_from_float(float f);
// 将软浮点数转换到浮点数
float sf32_to_float(sfloat32 sf);
// 判断是否是NaN
bool sf32_isnan(sfloat32 sf);
// 判断是否无穷
bool sf32_isinf(sfloat32 sf);
// 获取符号位
bool sf32_sign(sfloat32 sf);
其中在softfloat32.c文件中部分内容如下：
#include "softfloat.h"
#define SF32_INF 0x7f800000 // Inf，指数为0xff，有效位为0
#define SF32_NAN 0xffc00000 // NaN，指数为0xff，有效位不为0
#define SF32_MAX 0x7f7fffff // 单精度浮点表示的最大值
#define SF32_EXP(sf) (((sf) >> 23) & 0xff) // 取8位指数
#define SF32_SIG(sf) ((sf) & 0x7fffff) // 取23位有效位
#define SF32_SIGN(sf) ((sf) >> 31) // 取符号位
/* 构造一个32位浮点数 */
static inline sfloat32 sf32_pack(bool sign, uint16_t exp, uint32_t sig)
{
return ((uint32_t)sign << 31) | ((uint32_t)exp << 23) | (sig & 0x7fffff);
}
inline sfloat32 sf32_from_float(float f)
{
union
{
float f;
uint32_t u;
} usf;
usf.f = f;
return usf.u;
}
inline float sf32_to_float(sfloat32 sf)
{
union
{
float f;
uint32_t u;
} usf;
usf.u = sf;
return usf.f;
}
bool sf32_isnan(sfloat32 sf)
{
return (sf & 0x7fffffff) > SF32_INF;
}
bool sf32_isinf(sfloat32 sf)
{
return (sf & 0x7fffffff) == SF32_INF;
}
bool sf32_sign(sfloat32 sf)
{
return SF32_SIGN(sf);
}
这些是基本的操作，使用union进行变量类型的转换可以一定程度上避免ub。
比较运算
把比较运算放到前面，是因为这个实现起来非常简单。由于IEEE-754的指数是按照正常值加一个偏置值形成的，因为8位的指数在高位，而尾数在低位，所以除了一些特殊情况，按照定点数的规则即可比大小。
特殊情况是指NaN，任何数和NaN比较，只有不等于返回true，其他全是false，这个规则比较重要。
具体的函数声明如下：
// 判断相等
bool sf32_eq(sfloat32 sf1, sfloat32 sf2);
// 判断不相等
bool sf32_ne(sfloat32 sf1, sfloat32 sf2);
// 判断小于
bool sf32_lt(sfloat32 sf1, sfloat32 sf2);
// 判断小于等于
bool sf32_le(sfloat32 sf1, sfloat32 sf2);
// 判断大于
bool sf32_gt(sfloat32 sf1, sfloat32 sf2);
// 判断大于等于
bool sf32_ge(sfloat32 sf1, sfloat32 sf2);
以下是相等和不相等的实现，对于前面NaN的规则来说，不相等就是相等的简单取反：
bool sf32_eq(sfloat32 sf1, sfloat32 sf2)
{
if (sf32_isnan(sf1) || sf32_isnan(sf2))
return false;
return
(sf1 == sf2) || // 内容完全一样
(((sf1 | sf2) << 1) == 0); // 对于0则忽略符号位
}
bool sf32_ne(sfloat32 sf1, sfloat32 sf2)
{
return !sf32_eq(sf1, sf2);
}
但是对于其他操作数来说，则不是这样，比如大于不是小于等于的简单取反，而是要单独判断NaN，所以下面的代码也就好理解了：
static bool sf32_real_lt(sfloat32 sf1, sfloat32 sf2)
{
bool sign1 = SF32_SIGN(sf1);
bool sign2 = SF32_SIGN(sf2);
if (sign1 != sign2) // 符号不同，也可以用异或判断
// 若左边为负数且两数不为0，则说明左边小于右边
return sign1 && ((sf1 | sf2) << 1);
// 若符号相同，则左边小于右边的条件是两数不相等
// 同时，若两数为正数，绝对值小的小，否则绝对值大的小
return (sf1 != sf2) && (sign1 ^ (sf1 < sf2));
}
static bool sf32_real_le(sfloat32 sf1, sfloat32 sf2)
{
bool sign1 = SF32_SIGN(sf1);
bool sign2 = SF32_SIGN(sf2);
if (sign1 != sign2)
// 与之前不同，两数可以为0
return sign1 || (((sf1 | sf2) << 1) == 0);
// 与之前不同，相等或小于
return (sf1 == sf2) || (sign1 ^ (sf1 < sf2));
}
bool sf32_lt(sfloat32 sf1, sfloat32 sf2)
{
if (sf32_isnan(sf1) || sf32_isnan(sf2))
return false;
return sf32_real_lt(sf1, sf2);
}
bool sf32_le(sfloat32 sf1, sfloat32 sf2)
{
if (sf32_isnan(sf1) || sf32_isnan(sf2))
return false;
return sf32_real_le(sf1, sf2);
}
bool sf32_gt(sfloat32 sf1, sfloat32 sf2)
{
if (sf32_isnan(sf1) || sf32_isnan(sf2))
return false;
return !sf32_real_le(sf1, sf2);
}
bool sf32_ge(sfloat32 sf1, sfloat32 sf2)
{
if (sf32_isnan(sf1) || sf32_isnan(sf2))
return false;
return !sf32_real_lt(sf1, sf2);
}
在sf32_xx函数中，执行NaN的判断，而真正执行的sf32_real_xx函数里，就是纯粹的大小比较了。由于符号位的存在，实际大小比较需要针对同号和异号进行单独的判断，同时规定了正负0是相等的。注意到这些细节这块代码也就不难理解了。
舍入说明
在正式讲四则运算的代码之前，先讨论一下舍入的问题。因为想要获得较为精确的结果，必须妥善处理舍入的问题，以免出现精度损失导致误差变大。IEEE规定二进制浮点运算有4种舍入模式，其中默认的也是最常用的是向最近偶数舍入，同时由于向0舍入较为容易实现，所以本文也会涉及。至于上舍入和下舍入，本文不作讨论。
本文将向最近偶数舍入模式简称为舍入模式，将向0舍入模式简称为截断模式。
比如在加减法的时候，需要做对齐指令的操作，此时小的那个数的尾数会右移，这样必然导致部分尾数的丢弃，此时我们需要保留一部分被移出的位，以确保运算的精度，同时在最后的包装阶段统一进行舍入操作。
按照一些书籍上说的，需要引入保护位（Guard bit），舍入位（Round bit）和粘贴位（Sticky bit），合称为GRS位，以确保舍入过程可以保证一定的精度。
保护位其实就是右移后剩下的最后一位，舍入位则是被移出去的最左边一位，粘贴位则取决于其余被移出位的逻辑或关系（也可以理解为是否不为0）。
举个例子，假如有一个6位的数11 0100（其实就是52的二进制），我们将其右移5位，那么剩下的就是00 0001，因此其保护位就是最后的1；被移出去的是1 0100，其粘贴位就是被移出去的最左边位，也是1；其余4位0100的逻辑或是1，或者说移出去的0100不为0，所以粘贴位就是1。如果我们将这个6位的数右边补0扩展到8位，即1101 00 00，那么将其右移5位后的结果就是0000 01 11，多出来的RS位可以参与运算，从而保留了一定的精度。
实现右移保留粘贴位的具体代码如下：
// 对于32位数，右移保留粘贴位
static uint32_t shift_right_sticky32(uint32_t sig, uint16_t n)
{
uint32_t res = 0;
if (n < 31) { // 范围判断
res = sig >> n;
if (sig << (32 - n)) // 判断粘贴位
res = res | 1;
}
else {
if (sig) // 如果右移超过31位，只需判断是否存在粘贴位
res = 1;
}
return res;
}
// 对于64位数，右移保留粘贴位
static uint64_t shift_right_sticky64(uint64_t sig, uint16_t n)
{
uint64_t res = 0;
if (n < 63) {
res = sig >> n;
if (sig << (64 - n))
res = res | 1;
}
else {
if (sig)
res = 1;
}
return res;
}
上述两个函数分别针对32位和64位的数，用来实现右移过程中粘贴位的保留。
向最近偶数舍入
这个方法相对比较复杂，因为要减少舍入误差，所以所有运算必须要使用GRS位来进行精度控制。
理论上来说，只需额外保持舍入位和粘贴位就可以实现舍入，但是实际上为了保证精度，我们通常也会在最后的舍入阶段额外保留完整的GRS位。
以十进制为例，通常我们的四舍五入是遇到了0.5就会+1，但是实际上这样会导致一个问题：1-4舍，5-9入，导致舍入发生的概率并不相同。因此有这样一种舍入，叫四舍六入五成双。即1-4舍，6-9入，而5看左边一位是不是偶数，是偶数则舍，否则入。比如说1.5舍入后就是2，但是2.5舍入后也是2，也就是遇到5的情况下使得左边一位是偶数。
那么对于二进制来说，也就是要保证最低位是0，确保结果是偶数。比如1.1应该舍入为10.0，10.1也是10.0。
对于保留了GRS位的情况，则GRS=100需要判断奇偶，其他情况就简单舍入。
所以舍入代码如下：
// 对于32位数，进行向偶数舍入操作，低3位分别为GRS位
static uint32_t round_even_grs_32(uint32_t sig)
{
uint8_t grs = sig & 7;
sig = sig >> 3;
if ((grs > 4) || ((grs == 4) && (sig & 1)))
sig++;
return sig;
}
由于舍入通常是最后一步，所以在舍入后就可以顺便将GRS位丢弃了。
不过可能会存在舍入后进位的情况，也就是所有位都是1的时候，+1就会导致溢出。
向0舍入（截断）
这个是最简单的了，直接放弃多余的位，右移就行了。但是在减法对齐指数的时候，不能直接丢弃，而加法对齐却是可以丢弃，所以看似简单其实也有一些玄机。
还有一个，就是加法和乘法结果如果出现上溢，在此模式下结果应该设置为浮点数可以表示的最大值，而不是正负无穷，这点需要注意。
加减法运算
说实话，加减法应该是最难的，二者相对来说加法稍简单一些。
与定点数加减法不同，由于浮点数不是按照补码存储的，所以不能简单将加减法结合在一起，而是要分开计算。
对于加法来说，有以下4种情况：
- 正 + 正（1）
- 正 + 负（2）
- 负 + 正（3）
- 负 + 负（4）
对于情况（1）和情况（4）来说，它们的结果的符号是确定的，所以可以直接加法。而情况（2）和情况（3）则实际上是减法。
对于减法来说，则是如下4种情况：
- 正 - 正（5）
- 正 - 负（6）
- 负 - 正（7）
- 负 - 负（8）
对于情况（6）和情况（7），实际上结果的符号是确定的，所以也是加法就可以搞定。情况（5）和情况（8）才是真正的减法。
因此我们可以把上述的1、4、6、7这4种情况按照加法处理，符号则是左操作数的符号；把2、3、5、8这4种情况按照减法处理，符号则需要进一步比大小来判断。
可以用下面的代码来进一步包装实际的加减法：
sfloat32 sf32_add(sfloat32 sf1, sfloat32 sf2)
{
if (SF32_SIGN(sf1) == SF32_SIGN(sf2))
return real_sf32_add(sf1, sf2);
return real_sf32_sub(sf1, sf2);
}
sfloat32 sf32_sub(sfloat32 sf1, sfloat32 sf2)
{
if (SF32_SIGN(sf1) == SF32_SIGN(sf2))
return real_sf32_sub(sf1, sf2);
return real_sf32_add(sf1, sf2);
}
即对于加法操作，同号才是真正的加法，否则实际是减法；对于减法操作，同号才是真正的减法，否则实际是加法。
多嘴一句，上面的代码判断符号是否相等还可以继续优化，其在Mingw的64位环境下用gcc -O3生成的代码如下：
.seh_proc sf32_add
sf32_add:
.seh_endprologue
movl %ecx, %r8d
movl %edx, %eax
shrl $31, %r8d
shrl $31, %eax
cmpl %eax, %r8d
je .L242
jmp real_sf32_sub
.p2align 4,,10
.L242:
jmp real_sf32_add
.seh_endproc
实际上更好的写法应该是if ((sf1 ^ sf2) >> 31 == 0)或者if ((int32_t)(sf1 ^ sf2) >= 0)，后者更好一点（至少对于x86来说），不过对于gcc，都会生成如下代码：
.seh_proc sf32_add
sf32_add:
.seh_endprologue
movl %ecx, %eax
xorl %edx, %eax
jns .L242
jmp real_sf32_sub
.p2align 4,,10
.L242:
jmp real_sf32_add
.seh_endproc
这样编译器就能生成极致的代码了，直接少了3条指令，不过对于理解原理倒也没必要抠这么多细节，但是追求极致确是令人振奋的。
好吧扯得有点远了。
加法
加法的一般流程是：
- 拆开两个数，获取符号、指数、有效数
- 判断并处理特殊情况（0、无穷、NaN）
- 根据指数之差对齐尾数
- 对有效数做加法，并进行舍入操作
- 判断上溢
- 包装符号、指数、有效数
先把完整的代码贴出来，后面详细说明具体的部分：
static sfloat32 real_sf32_add(sfloat32 sf1, sfloat32 sf2)
{
int exp1 = SF32_EXP(sf1);
int exp2 = SF32_EXP(sf2);
uint32_t sig1 = SF32_SIG(sf1);
uint32_t sig2 = SF32_SIG(sf2);
bool sign = SF32_SIGN(sf1);
int expdiff = exp1 - exp2;
if (expdiff == 0) { // 两数指数相同
if (exp1 == 0) // 忽略非规格化数
return ((uint32_t)sign << 31);
if (exp1 == 0xff) // 判断无穷或NaN
return sf32_nan_inf(sign, sig1 | sig2);
}
if (exp1 == 0xff)
return sf32_nan_inf(sign, sig1);
if (exp2 == 0xff)
return sf32_nan_inf(sign, sig2);
if (exp1 == 0)
return (sf2 & 0x7fffffff) | ((uint32_t)sign << 31);
if (exp2 == 0)
return (sf1 & 0x7fffffff) | ((uint32_t)sign << 31);
sig1 |= 0x800000;
sig2 |= 0x800000;
#ifndef ROUND_TRUNC
sig1 = sig1 << 2;
sig2 = sig2 << 2;
#endif
if (expdiff < 0) {
expdiff = -expdiff;
#ifdef ROUND_TRUNC
if (expdiff > 24)
sig1 = 0;
else
sig1 = sig1 >> expdiff;
#else
sig1 = shift_right_sticky32(sig1, expdiff);
#endif
exp1 = exp2;
}
else if (expdiff > 0) {
#ifdef ROUND_TRUNC
if (expdiff > 24)
sig2 = 0;
else
sig2 = sig2 >> expdiff;
#else
sig2 = shift_right_sticky32(sig2, expdiff);
#endif
}
uint32_t sig = sig1 + sig2;
#ifdef ROUND_TRUNC
if (sig > 0xffffff) {
sig = sig >> 1;
exp1++;
}
#else
if (sig < 0x4000000)
sig = sig << 1;
else
exp1++;
sig = round_even_grs_32(sig);
if (sig > 0xffffff)
exp1++;
#endif
if (exp1 > 254)
#ifdef ROUND_TRUNC
return (SF32_MAX | ((uint32_t)sign << 31));
#else
return (SF32_INF | ((uint32_t)sign << 31));
#endif
return sf32_pack(sign, exp1, sig);
}
这里引入了函数sf32_nan_inf，其具体如下，主要是用来区分无穷和NaN：
static inline sfloat32 sf32_nan_inf(bool sign, uint32_t sig)
{
if (sig)
return SF32_NAN;
return (SF32_INF | ((uint32_t)sign << 31));
}
根据上述流程，我分块讲解一下real_sf32_add()的代码：
- 拆包，不多解释，符号就是第一个数的符号。
int exp1 = SF32_EXP(sf1);
int exp2 = SF32_EXP(sf2);
uint32_t sig1 = SF32_SIG(sf1);
uint32_t sig2 = SF32_SIG(sf2);
bool sign = SF32_SIGN(sf1);
- 特殊情况判断，对于非规格化数（即指数为0，不论尾数如何），将其视为0。
int expdiff = exp1 - exp2;
if (expdiff == 0) { // 两数指数相同
if (exp1 == 0) // 忽略非规格化数
return ((uint32_t)sign << 31);
if (exp1 == 0xff) // 判断无穷或NaN
return sf32_nan_inf(sign, sig1 | sig2);
}
if (exp1 == 0xff)
return sf32_nan_inf(sign, sig1);
if (exp2 == 0xff)
return sf32_nan_inf(sign, sig2);
if (exp1 == 0)
return (sf2 & 0x7fffffff) | ((uint32_t)sign << 31);
if (exp2 == 0)
return (sf1 & 0x7fffffff) | ((uint32_t)sign << 31);
上述的return sf32_nan_inf(sign, sig1 | sig2);这句对两数的有效数求算术或，这样可以判断两数是否存在NaN。
- 首先给有效数第24位补上规格化舍去的1，然后进行对齐操作。对于截断模式，不需要保留GRS位，而对于舍入模式，则需要保留GRS位，但是一开始仅保留2位，剩下1位在之后设置。
sig1 |= 0x800000;
sig2 |= 0x800000;
#ifndef ROUND_TRUNC
sig1 = sig1 << 2;
sig2 = sig2 << 2;
#endif
if (expdiff < 0) {
expdiff = -expdiff;
#ifdef ROUND_TRUNC
if (expdiff > 24)
sig1 = 0;
else
sig1 = sig1 >> expdiff;
#else
sig1 = shift_right_sticky32(sig1, expdiff);
#endif
exp1 = exp2;
}
else if (expdiff > 0) {
#ifdef ROUND_TRUNC
if (expdiff > 24)
sig2 = 0;
else
sig2 = sig2 >> expdiff;
#else
sig2 = shift_right_sticky32(sig2, expdiff);
#endif
}
这里需要注意一个细节，就是截断模式对齐如果指数差大于了24，实际上并不需要右移，直接置0即可。而对于舍入模式则需要保留粘贴位。
- 有效数相加，判断溢出。
uint32_t sig = sig1 + sig2;
#ifdef ROUND_TRUNC
if (sig > 0xffffff) {
sig = sig >> 1;
exp1++;
}
#else
if (sig < 0x4000000)
sig = sig << 1;
else
exp1++;
sig = round_even_grs_32(sig);
if (sig > 0xffffff)
exp1++;
#endif
对于截断模式，只需判断结果是否大于了24位，然后移掉多余的位；而舍入模式，需要判断结果是否大于27位，在舍入之后可能存在溢出，此时只需指数+1即可，甚至不需要右移，因为已经都是0了（只有全为1的情况下才会溢出，那么结果必然全为0，而规格化是不需要保留1的）。
由于之前左移了2位，因此实际上应该是2个26位的数相加，所以结果有可能是26位或27位。对于26位的情况，需要将其扩展到27位，而在舍入时会舍弃低3位，所以实际上结果是24位，这个过程中没有出现实际的溢出；对于27位的情况，由于出现了溢出，所以需要将指数+1。
在之前仅左移2位的目的就是这里，可以通过一句简单的if else将两种情况结合且不出现任何右移操作，这样在保证精度的同时还可以减少生成的代码体积。
可优化点：之前保留GRS位共3位用来舍入，这是理论上最少需要保留的位数，实际上我们完全可以先左移7位，使得尾数有31位，这样当相加产生溢出的时候就变成了32位，如果此时把这个数当作有符号数判断，那么可以发现它一定是小于0的，于是我们可以使用类似于if ((int32_t)a < 0)这样的写法，比if (a < 0x7ffffff)这样的代码生成的汇编更短。当然如果这么做需要同时修改round_even_grs_32函数。
- 判断上溢，打包返回结果
if (exp1 > 254)
#ifdef ROUND_TRUNC
return (SF32_MAX | ((uint32_t)sign << 31));
#else
return (SF32_INF | ((uint32_t)sign << 31));
#endif
return sf32_pack(sign, exp1, sig);
这个没什么好说的，关于上溢结果前面已经提到了。
减法
减法的流程和加法挺不一样的，具体体现在符号判断，大小判断等。
基本流程如下：
- 拆开两个数，获取符号、指数、有效数
- 判断指数差
- 指数相同，判断特殊值，判断尾数确定符号
- 指数不同，直接确定符号，判断特殊值
- 按照具体情况做减法并舍入
- 判断下溢
- 包装符号、指数、有效数
不过流程是这样说的，实际代码还是有点不一样，直接放出完整代码，具体看一下代码的注释：
static sfloat32 real_sf32_sub(sfloat32 sf1, sfloat32 sf2)
{
int exp1 = SF32_EXP(sf1);
int exp2 = SF32_EXP(sf2);
uint32_t sig1 = SF32_SIG(sf1);
uint32_t sig2 = SF32_SIG(sf2);
uint32_t sig;
bool sign = SF32_SIGN(sf1); // 先取被减数的符号
int expdiff = exp1 - exp2;
if (expdiff == 0) { // 指数相同
if (exp1 == 0xff) // 不需要判断是不是无穷，因为无穷减无穷也是NaN
return SF32_NAN;
if (sig1 == sig2) // 尾数相同
return 0;
if (sig1 > sig2) { // 大减小
sig = sig1 - sig2;
}
else {
sig = sig2 - sig1;
sign = !sign; // 由于被减数小于减数，需要翻转符号
}
sig = sig << 3; // 留出GRS位，虽然对这里没什么用
}
else if (expdiff < 0) { // 被减数小于减数
sign = !sign; // 翻转符号
if (exp2 == 0xff)
return sf32_nan_inf(sign, sig2);
if (exp1 == 0)
return (sf2 & 0x7fffffff) | ((uint32_t)sign << 31);
exp1 = exp2;
expdiff = -expdiff;
sig1 = (sig1 | 0x800000) << 3; // 补上24位的1，留出GRS位
sig1 = shift_right_sticky32(sig1, expdiff); // 保留粘贴位，很重要
sig = ((sig2 | 0x800000) << 3) - sig1; // 减数减去被减数
}
else {
if (exp1 == 0xff)
return sf32_nan_inf(sign, sig1);
if (exp2 == 0)
return (sf1 & 0x7fffffff) | ((uint32_t)sign << 31);
sig2 = (sig2 | 0x800000) << 3;
sig2 = shift_right_sticky32(sig2, expdiff);
sig = ((sig1 | 0x800000) << 3) - sig2; // 这里就正常减法
}
// 为了舍入和规格化，需要把结果补齐到27位
#if 1 // 使用前导0的方法
int shift = count_leading_zero(sig) - 5;
exp1 -= shift;
sig <<= shift;
#else // 不使用前导0
while (sig < 0x4000000) {
sig = sig << 1;
exp1--;
}
#endif
#ifdef ROUND_TRUNC
sig = sig >> 3;
#else
sig = round_even_grs_32(sig);
if (sig > 0xffffff) // 和之前一样
exp1++;
#endif
if (exp1 < 1) // 下溢判断
return ((uint32_t)sign << 31);
return sf32_pack(sign, exp1, sig);
}
代码中的count_leading_zero的一个通用实现如下（具体可以看 Hacker's Delight 的第 5.3 节）：
static int count_leading_zero(uint32_t x)
{
if (x == 0)
return 32;
int n = 1;
if ((x >> 16) == 0) { n += 16; x <<= 16; }
if ((x >> 24) == 0) { n += 8; x <<= 8; }
if ((x >> 28) == 0) { n += 4; x <<= 4; }
if ((x >> 30) == 0) { n += 2; x <<= 2; }
n -= (x >> 31);
return n;
}
对于ARM的一些CPU来说，存在汇编指令clz，可以直接获取一个寄存器的前导0个数。
对于Intel的CPU来说，使用bsr以及一系列操作可以达到同样的效果：
test eax,eax ; 假设这里eax是被求数
jnz @@BSR ; 如果eax非0，则求前导0个数
mov eax,32 ; 否则结果就是32
jmp @@END
@@BSR:
bsr eax,eax ; 获取eax最高位1的下标（从0开始）
xor eax,31 ; 与31异或，相当于用31去减
@@END:
这段汇编只是一个示例，仅供参考。
乘除法运算
相比于加减法，乘除法就简单了不少，其中乘法最容易实现。
和加减法不同的是，乘除法的符号位非常容易判断，对两个数的符号求异或即可。
乘法
乘法流程较为简单：
- 拆开两个数，获取符号、指数、有效数
- 判断特殊值（NaN，无穷，0）
- 指数相加，尾数相乘，结果舍入
- 判断上下溢
- 打包结果
特殊值情况相对复杂，这个涉及到NaN和无穷以及0的一系列操作。任何数乘NaN都是NaN，任何数乘无穷都是无穷（除了0，结果是NaN），0乘任何数都是0，代码实现中判断顺序很重要。
尾数都是24位的，取值范围是[0x800000~0xFFFFFF]，相乘会产生47~48位的结果，对于大部分32位CPU来说，都有64位数乘法的指令，所以可以直接计算，并没有非常大开销，而64位结果右移时的开销也很小。
完整代码如下，具体看注释，理解了加减法很多东西就没难度了：
sfloat32 sf32_mul(sfloat32 sf1, sfloat32 sf2)
{
int exp1 = SF32_EXP(sf1);
int exp2 = SF32_EXP(sf2);
uint32_t sig1 = SF32_SIG(sf1);
uint32_t sig2 = SF32_SIG(sf2);
bool sign = SF32_SIGN(sf1) ^ SF32_SIGN(sf2);
if (exp1 == 0xff) { // NaN或无穷
if (sig1 || ((exp2 == 0xff) && sig2) || (exp2 == 0)) // NaN*a或a*NaN或Inf*0
return SF32_NAN;
return (SF32_INF | ((uint32_t)sign << 31));
}
if (exp2 == 0xff) {
if (sig2 || (exp1 == 0)) // a*NaN或0*Inf
return SF32_NAN;
return (SF32_INF | ((uint32_t)sign << 31));
}
if ((exp1 == 0) || (exp2 == 0)) // 存在0
return ((uint32_t)sign << 31);
uint32_t sig;
int exp = exp1 + exp2 - 127; // 要减去偏移量
sig1 |= 0x800000;
sig2 |= 0x800000;
#ifdef ROUND_TRUNC
sig = ((uint64_t)sig1 * sig2) >> 23; // 相乘移位
if (sig > 0xffffff) { // 结果必然是24位或25位
sig = sig >> 1;
exp++;
}
#else
sig = shift_right_sticky64((uint64_t)sig1 * sig2, 21); // 少移2位，保留RS
if (sig < 0x4000000) // 这块代码和之前加法的操作如出一辙啊
sig = sig << 1;
else
exp++;
sig = round_even_grs_32(sig);
if (sig > 0xffffff)
exp++;
#endif
if (exp < 1)
return ((uint32_t)sign << 31);
if (exp > 254)
#ifdef ROUND_TRUNC // 和加法一样
return (SF32_MAX | ((uint32_t)sign << 31));
#else
return (SF32_INF | ((uint32_t)sign << 31));
#endif
return sf32_pack(sign, exp, sig);
}
和之前提到的差不多，bool sign = SF32_SIGN(sf1) ^ SF32_SIGN(sf2);可以优化为bool sign = (sf1 ^ sf2) >> 31;，显然指令至少会少一条。
除法
除法比乘法稍微复杂一些。
其基本流程和乘法很像：
- 拆开两个数，获取符号、指数、有效数
- 判断特殊值（NaN，无穷，0）
- 指数相减，尾数相除，结果舍入
- 判断上下溢
- 打包结果
特殊值只要搞清楚0作为被除数和除数的区别，以及0与无穷的一些关系，就不赘述了。
唯一稍复杂的是尾数的相除，由于尾数是定点小数，不能直接相除，否则会引入极大的误差，所以需要将被除数扩大再使用定点数的除法，因为24位的定点数想要不损失精度，至少需要扩展到48位，同时由于舍入问题，还需要获取余数。
和乘法不同的是，不是所有的32位CPU都支持64位被除数去除以一个32位的除数（也就是长除法）。对于英特尔x86来说，大家都很熟悉div指令，如果操作数是32位寄存器，那么会将eax作为低位，edx作为高位去除以这个寄存器的值，将商存储在eax，余数存储在edx。而对于其他一些32位CPU来说，不一定支持这样的指令。
所以需要自行实现一个长除法，而且为了使用这个除法操作，需要确保两数相除的商要在32位。这实际上不难实现，请看如下代码：
static inline uint32_t ulong_div_mod(uint64_t n, uint32_t base, uint32_t *r)
{
uint64_t rem = n;
uint32_t quo = 0;
for (int i = 31; i >= 0; i--) {
uint64_t b = (uint64_t)base << i;
if (rem > b) {
rem -= b;
quo |= (uint32_t)(1 << i);
}
}
*r = rem; // 需要自行确保传入的参数非空，可以少一个判断
return quo;
}
其实原理很简单啊，很容易想到的，就是移位，依次比大小，减得下就减，然后给结果置位。
下面是一个更好的代码（对于32位CPU来说），从文末推荐的 Hacker's Delight 这本书学到的，有兴趣可以去看看，在书上 9.4.1 章，其实用的就是非常经典的恢复余数法：
static inline uint32_t ulong_div_mod(uint64_t n, uint32_t base, uint32_t *r)
{
uint32_t t, x, y; // x高32位，y低32位
x = n >> 32;
y = n;
for (int i = 1; i <= 32; i++) {
t = x >> 31;
x = (x << 1) | (y >> 31);
y = y << 1;
if ((x | t) >= base) {
x -= base;
y++;
}
}
*r = x;
return y;
}
其实书上第 9.4.2 章所讲到的短除法实现长除法对于有32位硬件除法的CPU来说是更高效的，不过这里就不深入了。
虽然这样实现了长除法，但是效率仍然很低，即使使用了两次短除法实现一次长除法，除了两次除法指令的开销，还有多次乘法的开销，而且代码量还大了不少。
对于不支持长除法的CPU来说，还有一种方法可以快速求的定点小数的除法，而且不需要依赖除法指令，只需要数次乘法和一系列基础的指令即可。这便是使用牛顿迭代法求近似倒数并计算尾数的除法。有关牛顿法求倒数的介绍，在 https://www.52pojie.cn/thread-1839972-1-1.html 可以找到。
这里放一下这里需要用到的求倒数的代码：
const uint8_t initreciptable8x8[8] = {
0xf1, 0xd8, 0xc3, 0xb2, 0xa4, 0x98, 0x8d, 0x84
};
static inline uint32_t approx_recip(uint32_t n)
{
uint32_t x, t;
x = initreciptable8x8[(n >> 28) & 7] << 24;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
return x;
}
使用长除法或者牛顿迭代求倒数法，就可以实现浮点数的除法操作了：
sfloat32 sf32_div(sfloat32 sf1, sfloat32 sf2)
{
int exp1 = SF32_EXP(sf1);
int exp2 = SF32_EXP(sf2);
uint32_t sig1 = SF32_SIG(sf1);
uint32_t sig2 = SF32_SIG(sf2);
bool sign = SF32_SIGN(sf1) ^ SF32_SIGN(sf2);
if (exp1 == 0xff) {
if (sig1 || (exp2 == 0xff)) // NaN/a或Inf/NaN或Inf/Inf
return SF32_NAN;
return (SF32_INF | ((uint32_t)sign << 31));
}
if (exp2 == 0xff) { // a/NaN
if (sig2)
return SF32_NAN;
return ((uint32_t)sign << 31);
}
if (exp1 == 0) {
if (exp2 == 0) // 0/0
return SF32_NAN;
return ((uint32_t)sign << 31); // 0/a
}
if (exp2 == 0) // a/0
return (SF32_INF | ((uint32_t)sign << 31));
int exp = exp1 - exp2 + 127; // 指数相减
sig1 = sig1 | 0x800000;
sig2 = sig2 | 0x800000;
uint32_t sig;
#ifdef LONG_DIV // 使用长除法实现的代码
uint64_t sig64;
if (sig1 < sig2) { // 比较大小，这是为了确保商一定是32位的
exp--;
sig64 = (uint64_t)sig1 << 32;
}
else {
sig64 = (uint64_t)sig1 << 31;
}
#ifdef ROUND_TRUNC
sig = ulong_div(sig64, sig2) >> 8; // 直接右移
#else
uint32_t rem;
sig = ulong_div_mod(sig64, sig2, &rem) >> 5; // 右移5位，留出GRS
if (rem) // 余数存在最低位置1，相当于粘贴位
sig |= 1;
sig = round_even_grs_32(sig); // 舍入处理
#endif
#else // !LONG_DIV，即不使用长除法，而是倒数法
if (sig1 < sig2) { // 和之前思路一样，确保不溢出
exp--;
sig1 <<= 8;
}
else {
sig1 <<= 7;
}
sig2 <<= 8; // 从Q1.23扩展到Q1.31
sig = (uint64_t)sig1 * approx_recip(sig2) >> 31; // 求近似的商
uint64_t rem = (uint64_t)sig1 << 32;
rem -= (uint64_t)sig * sig2; // 计算余数
while (rem >= sig2) { // 修正商
sig++;
rem -= sig2;
}
#ifdef ROUND_TRUNC
sig >>= 8;
#else
sig >>= 5;
if (rem)
sig |= 1;
sig = round_even_grs_32(sig);
#endif
#endif // LONG_DIV
if (exp < 1)
return ((uint32_t)sign << 31);
if (exp > 254)
return (SF32_INF | ((uint32_t)sign << 31));
return sf32_pack(sign, exp, sig);
}
在长除法的代码中，对于截断模式，直接使用ulong_div，该函数和ulong_div_mod完全一样，但是少了一个余数的参数，生成的代码更精简一点（不过不这样写编译器大概率会优化，所以这样做聊胜于无吧）。在倒数法中，都需要求余数来修正估算的商，所以区别就在舍入部分。
和整数的转换
这个比较常用，也比较简单。
int32_t sf32_to_i32(sfloat32 sf)
{
int exp = SF32_EXP(sf) - 127;
if (exp < 0) // 指数过小返回0
return 0;
bool sign = SF32_SIGN(sf);
if (exp > 30) // 指数过大就返回负数最大值，也是主流做法
return (int32_t)0x80000000;
uint32_t sig = SF32_SIG(sf) | 0x800000;
int shift = 23 - exp;
int32_t res;
if (shift > 0)
res = sig >> shift;
else
res = sig << (-shift);
if (sign)
res = -res;
return res;
}
sfloat32 sf32_from_i32(int32_t i)
{
if (i == 0)
return 0;
bool sign = i >> 31;
if (sign)
i = -i;
int exp = 127 + 23;
#if 1
int shift = count_leading_zero(i) - 8; // 计算移位数
if (shift > 0) { // 小于0x800000左移
i <<= shift;
exp -= shift;
}
else if (shift < 0) { // 大于0x1000000右移
shift = -shift;
i >>= shift;
exp += shift;
}
#else
while (i >= 0x1000000) { // 大了就右移，指数+1
i = i >> 1;
exp++;
}
while (i < 0x800000) { // 小了就左移，指数-1
i = i << 1;
exp--;
}
#endif
return sf32_pack(sign, exp, i);
}
对于64位整数同理，只需要适当修改一下代码即可。
使用前导0可以获得更快的速度，即使会造成一定程度上的指令数量膨胀，而且实际上可以针对不同CPU使用专门的方法，效果是明显好于循环的。
文末感言
有些浮点数的细节是很难完全说清楚的，这个需要自己实践一下才知道，比如我就是和CPU自带的浮点单元的结果进行比较判断是不是没有问题，我甚至专门写了个随机数的压力测试，基本上能扛住较大的次数基本就说明正常的运算没有问题了。当然如果出现问题，就需要大量的调试操作了，因为这些运算细节很少有人说，所以只能一点一点摸索出来。
网上好多资料都是和计组有关的，极少数和硬件设计Verilog有关系，基本上就没有软件浮点数的。所以具体的参照非常少，也导致写起来非常缓慢。
而且IEEE-754 2008规定的基本运算平方根我也暂时没有高效率的实现，所以就不放出来献丑了，目前可以选择广为流传的卡马克平方根倒数的方法，直接二进制操作的方法目前尚未研究，毕竟不是学硬件的，就算学硬件的现在有那么多IP可以调谁还写这些是吧。
照着我这个思路，半精度浮点，双精度浮点应该都不难实现，实际上难的是怎么优化得效率更高一些，比如双精度浮点必然涉及更大的数的乘除法，这些如何优化是个问题。
结尾推荐一些好书，比如 Computer Arithmetic Algorithms and Hardware Designs 2nd edition，作者 Behrooz Parhami，这本书基本是属于硬件设计那块的了，不过也可以看看了解硬件的实现原理。还有 Hacker's Delight，这本也非常不错，我在论坛电子书屋板块也发布了，登录了论坛的朋友在论坛搜索“算法心得：高效算法的奥秘”就可以找到了。
尤其这本 Hacker's Delight，我最近在看，以后有什么好的思路会更新一下。
由于书签存在乱码问题，所以 Computer Arithmetic Algorithms and Hardware Designs 一书目前没有发在论坛，等我有空重新制作书签后再发。
[复制链接]
发表于 2023-9-6 12:20
|
发表于 2023-10-7 18:27
