吾爱破解 - LCG - LSG |安卓破解|病毒分析|www.52pojie.cn

 找回密码
 注册[Register]

QQ登录

只需一步,快速开始

查看: 4111|回复: 16
上一主题 下一主题
收起左侧

[C&C++ 原创] 软件模拟实现IEEE-754单精度浮点数运算

  [复制链接]
跳转到指定楼层
楼主
DEATHTOUCH 发表于 2023-9-6 12:20 回帖奖励
本帖最后由 DEATHTOUCH 于 2024-2-28 12:55 编辑

软件模拟实现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条指令,不过对于理解原理倒也没必要抠这么多细节,但是追求极致确是令人振奋的。

好吧扯得有点远了。

加法

加法的一般流程是:

  1. 拆开两个数,获取符号、指数、有效数
  2. 判断并处理特殊情况(0、无穷、NaN)
  3. 根据指数之差对齐尾数
  4. 对有效数做加法,并进行舍入操作
  5. 判断上溢
  6. 包装符号、指数、有效数

先把完整的代码贴出来,后面详细说明具体的部分:

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()的代码:

  1. 拆包,不多解释,符号就是第一个数的符号。
    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);
  1. 特殊情况判断,对于非规格化数(即指数为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。

  1. 首先给有效数第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即可。而对于舍入模式则需要保留粘贴位。

  1. 有效数相加,判断溢出。
    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函数。

  1. 判断上溢,打包返回结果
    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);

这个没什么好说的,关于上溢结果前面已经提到了。

减法

减法的流程和加法挺不一样的,具体体现在符号判断,大小判断等。

基本流程如下:

  1. 拆开两个数,获取符号、指数、有效数
  2. 判断指数差
  3. 指数相同,判断特殊值,判断尾数确定符号
  4. 指数不同,直接确定符号,判断特殊值
  5. 按照具体情况做减法并舍入
  6. 判断下溢
  7. 包装符号、指数、有效数

不过流程是这样说的,实际代码还是有点不一样,直接放出完整代码,具体看一下代码的注释:

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:

这段汇编只是一个示例,仅供参考。

乘除法运算

相比于加减法,乘除法就简单了不少,其中乘法最容易实现。

和加减法不同的是,乘除法的符号位非常容易判断,对两个数的符号求异或即可。

乘法

乘法流程较为简单:

  1. 拆开两个数,获取符号、指数、有效数
  2. 判断特殊值(NaN,无穷,0)
  3. 指数相加,尾数相乘,结果舍入
  4. 判断上下溢
  5. 打包结果

特殊值情况相对复杂,这个涉及到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;,显然指令至少会少一条。

除法

除法比乘法稍微复杂一些。

其基本流程和乘法很像:

  1. 拆开两个数,获取符号、指数、有效数
  2. 判断特殊值(NaN,无穷,0)
  3. 指数相减,尾数相除,结果舍入
  4. 判断上下溢
  5. 打包结果

特殊值只要搞清楚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 一书目前没有发在论坛,等我有空重新制作书签后再发。

免费评分

参与人数 15威望 +3 吾爱币 +135 热心值 +15 收起 理由
chuan9 + 1 + 1 谢谢@Thanks!
sam喵喵 + 1 推荐的第一本书希望楼主也可以在论坛发布一下谢谢
schtg + 1 + 1 用心讨论,共获提升!
sanmylc + 1 + 1 用心讨论,共获提升!
allspark + 1 + 1 用心讨论,共获提升!
fengbolee + 2 + 1 欢迎分析讨论交流,吾爱破解论坛有你更精彩!
Zz4794zZ + 1 + 1 谢谢@Thanks!
Bluesky10 + 1 + 1 热心回复!
canty胖胖 + 1 + 1 用心讨论,共获提升!
ilil + 1 + 1 谢谢@Thanks!
笙若 + 1 + 1 谢谢@Thanks!
月夜克星 + 2 + 1 厉害 而且相当有耐心了 之前用到点相关的东西 研究了好几天才凑合着实现了.
侃遍天下无二人 + 2 + 100 + 1 感谢发布原创作品,吾爱破解论坛因你更精彩!
wushaominkk + 1 + 20 + 1 感谢发布原创作品,吾爱破解论坛因你更精彩!
lyl610abc + 2 + 1 有点底层啊,666

查看全部评分

本帖被以下淘专辑推荐:

发帖前要善用论坛搜索功能,那里可能会有你要找的答案或者已经有人发布过相同内容了,请勿重复发帖。

来自 15#
 楼主| DEATHTOUCH 发表于 2023-10-7 18:27 |楼主
双精度double的加减乘除已经基本搞定了,到时候研究一下平方根,再整一个精简方便的软浮点库出来。
推荐
 楼主| DEATHTOUCH 发表于 2023-9-6 14:51 |楼主
本帖最后由 DEATHTOUCH 于 2023-9-6 15:13 编辑
weiya909 发表于 2023-9-6 13:55
楼总,你的代码写得很六,就是有点费眼。

论坛的markdown支持似乎不是那么好,也有可能我不太会用

我在代码仓库也上传了一份,可以看看 tutorial.md · Peacoor/softfloat32 - Gitee.com
推荐
月夜克星 发表于 2023-9-7 19:58
厉害 而且相当有耐心了 之前用到点相关的东西 研究了好几天才凑合着实现了一部分 感谢作者分享出来
沙发
zhu639 发表于 2023-9-6 13:25
专业支持一下
3#
weiya909 发表于 2023-9-6 13:55
楼总,你的代码写得很六,就是有点费眼。
5#
KMT-jiang 发表于 2023-9-6 15:15
技术大佬,牛批
6#
bjxiaoyao 发表于 2023-9-6 16:01
大佬厉害,很专业
7#
亡屿歌 发表于 2023-9-7 09:46
写得不错,gitee上看起来效果很好&#128077;。
9#
ttspop 发表于 2023-9-7 21:09
单精度是不是可以用IQMATH
10#
 楼主| DEATHTOUCH 发表于 2023-9-7 22:41 |楼主
ttspop 发表于 2023-9-7 21:09
单精度是不是可以用IQMATH

查了一下好像挺不错,不过不知道有没有源码,网上东西有点少
您需要登录后才可以回帖 登录 | 注册[Register]

本版积分规则

快速回复 收藏帖子 返回列表 搜索

RSS订阅|小黑屋|处罚记录|联系我们|吾爱破解 - LCG - LSG ( 京ICP备16042023号 | 京公网安备 11010502030087号 )

GMT+8, 2024-7-23 17:18

Powered by Discuz!

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表