浮点数算法的内部实现

科学计算当中会用到不少浮点数的操作，这些浮点数可能是16位，32位，64位，80位甚至是128位。开源项目SoftFloat提供了一个高效的浮点运算实现，可以在没有硬件支持的情况下，高效模拟浮点数的各种操作。

那么，浮点数之间的比较，基本运算这些究竟是怎么实现的呢，可以拿32位浮点数作为例子。

这是32位浮点数加法的实现，首先声明了一个结构体float32_t。

typedef struct { uint32_t v; } float32_t;

这提供了32位浮点数的底层位表示，同时还声明了一个union。

union ui32_f32 { uint32_t ui; float32_t f; };

一方面保存了浮点数的位表示，另一方面也可以转换为32位无符号整型直接进行比较，这在后面的算法当中会直接涉及。先看看加法。

float32_t f32_add( float32_t a, float32_t b )
{
    union ui32_f32 uA;
    uint_fast32_t uiA;
    union ui32_f32 uB;
    uint_fast32_t uiB;
#if ! defined INLINE_LEVEL || (INLINE_LEVEL < 1)
    float32_t (*magsFuncPtr)( uint_fast32_t, uint_fast32_t );
#endif

    uA.f = a;
    uiA = uA.ui;
    uB.f = b;
    uiB = uB.ui;
#if defined INLINE_LEVEL && (1 <= INLINE_LEVEL)
    if ( signF32UI( uiA ^ uiB ) ) {
        return softfloat_subMagsF32( uiA, uiB );
    } else {
        return softfloat_addMagsF32( uiA, uiB );
    }
#else
    magsFuncPtr =
        signF32UI( uiA ^ uiB ) ? softfloat_subMagsF32 : softfloat_addMagsF32;
    return (*magsFuncPtr)( uiA, uiB );
#endif

}

这里uiA和uiB是存储无符号整型的，signF32UI是提取符号位的。signF32UI(uiA ^ uiB)判断符号位是否相同，如果相同则调用加法，如果符号位不相同则调用减法，因为没有浮点数，所以只能通过整型去模拟，另外，union存储浮点和整型有一个名词，似乎叫类型双关技术？不过这里union存储的只是位表示，并不是真的浮点数。

float32_t f32_sub( float32_t a, float32_t b )
{
    union ui32_f32 uA;
    uint_fast32_t uiA;
    union ui32_f32 uB;
    uint_fast32_t uiB;
#if ! defined INLINE_LEVEL || (INLINE_LEVEL < 1)
    float32_t (*magsFuncPtr)( uint_fast32_t, uint_fast32_t );
#endif

    uA.f = a;
    uiA = uA.ui;
    uB.f = b;
    uiB = uB.ui;
#if defined INLINE_LEVEL && (1 <= INLINE_LEVEL)
    if ( signF32UI( uiA ^ uiB ) ) {
        return softfloat_addMagsF32( uiA, uiB );
    } else {
        return softfloat_subMagsF32( uiA, uiB );
    }
#else
    magsFuncPtr =
        signF32UI( uiA ^ uiB ) ? softfloat_addMagsF32 : softfloat_subMagsF32;
    return (*magsFuncPtr)( uiA, uiB );
#endif

}

减法则是在判断符号那里反过来，其它一样。这时候可以看看比较运算怎么做。

bool f32_le( float32_t a, float32_t b )
{
    union ui32_f32 uA;
    uint_fast32_t uiA;
    union ui32_f32 uB;
    uint_fast32_t uiB;
    bool signA, signB;

    uA.f = a;
    uiA = uA.ui;
    uB.f = b;
    uiB = uB.ui;
    if ( isNaNF32UI( uiA ) || isNaNF32UI( uiB ) ) {
        softfloat_raiseFlags( softfloat_flag_invalid );
        return false;
    }
    signA = signF32UI( uiA );
    signB = signF32UI( uiB );
    return
        (signA != signB) ? signA || ! (uint32_t) ((uiA | uiB)<<1)
            : (uiA == uiB) || (signA ^ (uiA < uiB));

}

最后的表达式有点绕，一步一步拆分。首先符号不相等（一正一负）的话，如果A的符号是1，也就是负数，肯定比B小，否则走 || 后的分支。把A和B的最高位（符号位）剔除，判断是否相同，也就是+0和-0的情况，这里记得别漏了前面的!符号，因为判断两者是否都为0；如果A和B同号的话，如果都是正数则直接比较，如果都是负数，则前面的signA会对结果取反。