浮点数算法的内部实现
科学计算当中会用到不少浮点数的操作,这些浮点数可能是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会对结果取反。
结语
最近处于校招阶段,正在准备,有时间会分享自己的心得和体会,希望尽早上岸。