为 QEMU softfloat 添加用于神经网络计算的浮点精度

QEMU 的 softfloat 源码位于 fpu/ 和 include/fpu/ 路径中，代码最初源自 Berkeley SoftFloat IEC/IEEE 浮点运算包的 2a 版本（SoftFloat-2a），后续经过 QEMU 项目贡献者修改。

目前我已经为 softfloat 添加了 tfloat32 和 float8e4m3 和 float8e5m2 的支持，可以从下面的仓库获取：

git clone -b neural-network-softfloat https://gitee.com/gevico/qemu.git

源文件结构如下：

include/fpu
|-- softfloat-helpers.h # standalone helpers，用于初始化，比如设置舍入模式、异常模式等
|-- softfloat-macros.h # QEMU 浮点支持宏
|-- softfloat-types.h # 定义浮点精度和格式
`-- softfloat.h  # 对外使用的浮点接口，比如浮点运算，浮点精度转换

fpu/
|-- meson.build # 构建脚本
|-- softfloat-parts-addsub.c.inc
|-- softfloat-parts.c.inc
|-- softfloat-specialize.c.inc # 浮点运算的一些特殊实现
`-- softfloat.c # 浮点功能的核心源码

为 softfloat 添加一个新的浮点精度，分下面几步完成：

在 softfloat-types.h 中增加新的浮点精度定义
在 softfloat.c 中实现关键浮点运算函数
在 softfloat-specialize.c.inc 中实现特殊的浮点运算函数
在 softfloat.h 中添加浮点精度的接口声明

下面我们以添加 tfloat32（Tensor float32， TF32）浮点精度为例，进行讲解。

一、在 softfloat-types.h 中增加新的浮点精度定义

tfloat32 的浮点格式如下：

tfloat32 内部计算时仍然按照 float32 进行运算，只不过输出的浮点值，仅保留 float32 小数位的前 10 位。

我们来添加 tfloat32 的定义。

// path: include/fpu/softfoloat-types.h

/*
 * Software neural-network floating-point types.
 */
typedef uint16_t bfloat16;
typedef uint32_t tfloat32;  /* 添加 TF32 的定义 */

/* 一些操作宏，参考其他精度定义 */
#define tfloat32_val(x) (x)
#define make_tfloat32(x) (x)

二、在 softfloat.c 中实现关键浮点运算函数

为了提高可移植性和代码复用，softfloat 采用 packet 和 unpacket 的形式，对浮点格式进行解包和打包，具体运算是将浮点值解包以后，对浮点位、指数位、小数位进行分别计算。

因此我们需要先添加 tfloat32 的 float format 的定义，代码如下：

// path: fpu/softfloat.c

/* Expand fields based on the size of exponent and fraction */
#define FLOAT_PARAMS_(E)                                \
    .exp_size       = E,                                \
    .exp_bias       = ((1 << E) - 1) >> 1,              \
    .exp_re_bias    = (1 << (E - 1)) + (1 << (E - 2)),  \
    .exp_max        = (1 << E) - 1

#define FLOAT_PARAMS(E, F)                              \
    FLOAT_PARAMS_(E),                                   \
    .frac_size      = F,                                \
    .frac_shift     = (-F - 1) & 63,                    \
    .round_mask     = (1ull << ((-F - 1) & 63)) - 1

/* 定义 tfloat32 的 floatfmt */
static const FloatFmt tfloat32_params = {
    FLOAT_PARAMS(8, 23)
};

/* Unpack a float to parts, but do not canonicalize.  */
static void unpack_raw64(FloatParts64 *r, const FloatFmt *fmt, uint64_t raw)
{
    const int f_size = fmt->frac_size;
    const int e_size = fmt->exp_size;

    *r = (FloatParts64) {
        .cls = float_class_unclassified,
        .sign = extract64(raw, f_size + e_size, 1),
        .exp = extract64(raw, f_size, e_size),
        .frac = extract64(raw, 0, f_size)
    };
}

/* Pack a float from parts, but do not canonicalize.  */
static uint64_t pack_raw64(const FloatParts64 *p, const FloatFmt *fmt)
{
    const int f_size = fmt->frac_size;
    const int e_size = fmt->exp_size;
    uint64_t ret;

    ret = (uint64_t)p->sign << (f_size + e_size);
    ret = deposit64(ret, f_size, e_size, p->exp);
    ret = deposit64(ret, 0, f_size, p->frac);
    return ret;
}

/* 定义解包函数 */
static void QEMU_FLATTEN tfloat32_unpack_raw(FloatParts64 *p, tfloat32 f)
{
    unpack_raw64(p, &tfloat32_params, f);
}

static void tfloat32_unpack_canonical(FloatParts64 *p, tfloat32 f,
                                                             float_status *s)
{
    tfloat32_unpack_raw(p, f);
    parts_canonicalize(p, s, &tfloat32_params);
}

/* 定义组包函数，组包时只保留小数位前10位  */
static tfloat32 QEMU_FLATTEN tfloat32_pack_raw(const FloatParts64 *p)
{
    /* The fraction is kept to only its first 10 bits. */
    return pack_raw64(p, &tfloat32_params) & 0xFFFC0000;
}

/* 增加组包的舍入方式的封装 */
static tfloat32 tfloat32_round_pack_canonical(FloatParts64 *p,
                                  float_status *s)
{
    parts_uncanon(p, s, &tfloat32_params);
    return tfloat32_pack_raw(p);
}

然后我们添加具体的计算函数，数量比较多，我们以 addsub 和精度转换为例：

static tfloat32 QEMU_FLATTEN
tfloat32_addsub(tfloat32 a, tfloat32 b, float_status *status, bool subtract)
{
    FloatParts64 pa, pb, *pr;

    tfloat32_unpack_canonical(&pa, a, status);
    tfloat32_unpack_canonical(&pb, b, status);
    pr = parts_addsub(&pa, &pb, status, subtract);

    return tfloat32_round_pack_canonical(pr, status);
}

tfloat32 tfloat32_add(tfloat32 a, tfloat32 b, float_status *status)
{
    return tfloat32_addsub(a, b, status, false);
}

tfloat32 tfloat32_sub(tfloat32 a, tfloat32 b, float_status *status)
{
    return tfloat32_addsub(a, b, status, true);
}

/* 以 tfloat32 和 float32 的精度相互转换为例 */

float32 tfloat32_to_float32(tfloat32 a, float_status *s)
{
    FloatParts64 p;

    tfloat32_unpack_canonical(&p, a, s);
    parts_float_to_float(&p, s);
    return float32_round_pack_canonical(&p, s);
}

tfloat32 float32_to_tfloat32(float32 a, float_status *s)
{
    FloatParts64 p;

    float32_unpack_canonical(&p, a, s);
    parts_float_to_float(&p, s);
    return tfloat32_round_pack_canonical(&p, s);
}

三、在 softfloat-specialize.c.inc 中实现特殊的浮点运算函数

需要实现的函数只有两个，tfloat32_is_quiet_nan() 和 tfloat32_is_signaling_nan()。

/*----------------------------------------------------------------------------
| Returns 1 if the tfloat32 value `a' is a quiet
| NaN; otherwise returns 0.
*----------------------------------------------------------------------------*/

bool tfloat32_is_quiet_nan(tfloat32 a_, float_status *status)
{
    if (no_signaling_nans(status)) {
        return tfloat32_is_any_nan(a_);
    } else {
        uint32_t a = tfloat32_val(a_);
        if (snan_bit_is_one(status)) {
            return (((a >> 22) & 0x1FF) == 0x1FE) && (a & 0x003FFFFF);
        } else {
            return ((uint32_t)(a << 1) >= 0xFF800000);
        }
    }
}

/*----------------------------------------------------------------------------
| Returns 1 if the tfloat32 value `a' is a signaling
| NaN; otherwise returns 0.
*----------------------------------------------------------------------------*/

bool tfloat32_is_signaling_nan(tfloat32 a_, float_status *status)
{
    if (no_signaling_nans(status)) {
        return 0;
    } else {
        uint32_t a = tfloat32_val(a_);
        if (snan_bit_is_one(status)) {
            return ((uint32_t)(a << 1) >= 0xFF800000);
        } else {
            return (((a >> 22) & 0x1FF) == 0x1FE) && (a & 0x003FFFFF);
        }
    }
}

四、在 softfloat.h 中添加浮点精度的接口声明

除了要声明在 softfloat.c 和 softfloat-specialize.c.inc 中定义的 tfloat32 相关的函数，在 softfloat.h 头文件里还需要添加一些常用的内联函数。

代码实现如下, 摘要重要部分：

static inline tfloat32 tfloat32_abs(tfloat32 a)
{
    /* Note that abs does *not* handle NaN specially, nor does
     * it flush denormal inputs to zero.
     */
    return make_tfloat32(tfloat32_val(a) & 0x7fffffff);
}

static inline tfloat32 tfloat32_chs(tfloat32 a)
{
    /* Note that chs does *not* handle NaN specially, nor does
     * it flush denormal inputs to zero.
     */
    return make_tfloat32(tfloat32_val(a) ^ 0x80000000);
}

static inline bool tfloat32_is_infinity(tfloat32 a)
{
    return (tfloat32_val(a) & 0x7fffffff) == 0x7f800000;
}

static inline bool tfloat32_is_neg(tfloat32 a)
{
    return tfloat32_val(a) >> 31;
}

static inline bool tfloat32_is_zero(tfloat32 a)
{
    return (tfloat32_val(a) & 0x7fffffff) == 0;
}

static inline bool tfloat32_is_any_nan(tfloat32 a)
{
    return ((tfloat32_val(a) & ~(1 << 31)) > 0x7f800000UL);
}

static inline bool tfloat32_is_zero_or_denormal(tfloat32 a)
{
    return (tfloat32_val(a) & 0x7f800000) == 0;
}

static inline bool tfloat32_is_normal(tfloat32 a)
{
    return (((tfloat32_val(a) >> 23) + 1) & 0xff) >= 2;
}

static inline bool tfloat32_is_denormal(tfloat32 a)
{
    return tfloat32_is_zero_or_denormal(a) && !tfloat32_is_zero(a);
}

static inline bool tfloat32_is_zero_or_normal(tfloat32 a)
{
    return tfloat32_is_normal(a) || tfloat32_is_zero(a);
}

static inline tfloat32 tfloat32_set_sign(tfloat32 a, int sign)
{
    return make_tfloat32((tfloat32_val(a) & 0x7fffffff) | (sign << 31));
}

static inline bool tfloat32_eq(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare(a, b, s) == float_relation_equal;
}

static inline bool tfloat32_le(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare(a, b, s) <= float_relation_equal;
}

static inline bool tfloat32_lt(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare(a, b, s) < float_relation_equal;
}

static inline bool tfloat32_unordered(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare(a, b, s) == float_relation_unordered;
}

static inline bool tfloat32_eq_quiet(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare_quiet(a, b, s) == float_relation_equal;
}

static inline bool tfloat32_le_quiet(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare_quiet(a, b, s) <= float_relation_equal;
}

static inline bool tfloat32_lt_quiet(tfloat32 a, tfloat32 b, float_status *s)
{
    return tfloat32_compare_quiet(a, b, s) < float_relation_equal;
}

static inline bool tfloat32_unordered_quiet(tfloat32 a, tfloat32 b,
                                           float_status *s)
{
    return tfloat32_compare_quiet(a, b, s) == float_relation_unordered;
}

#define tfloat32_zero 0
#define tfloat32_half 0x3f000000
#define tfloat32_one 0x3f800000
#define tfloat32_one_point_five 0x3fc00000
#define tfloat32_two 0x40000000
#define tfloat32_three 0x40400000
#define tfloat32_infinity 0x7f800000

后续文章，我们将讲解如何基于 qemu 的 fpu 测试框架，验证我们新增浮点精度的正确性。