[注明] 现有的Intel中的SSE指令 cvtps2dq XMM,XMM/m128 支持把源存储器4个单精度浮点数变成4个双字有符号整数,结果送入目的寄存器,内存变量必须对齐内存16字节,另外还有其他指令支持双精度和整型的转换。

转载本文,作为C算法设计的探讨。

在计算机图形运算中,常常要将浮点数转换为整数,例如在图像的光栅化阶段,就要执行大量的类型转换,以便将浮点数表示的坐标转化为整数表示的屏幕坐标。Ok,it's so easy:
----------------------------------------------------------------------------------------
//
// 强制类型转换
// 小数部分将被裁剪掉
//
int_val = (int)float_val;
----------------------------------------------------------------------------------------
嘿嘿,很高兴你居然和我一样单纯!这个操作实在是太TINY了,以至于我从来没想过它是怎么实现的,直到某天某个家伙跟我说,不要使用标准C类型转换,因为那太慢了!我当时的震惊不下于倒霉的冒险者遇上了龙。

标准C类型转换最大的优点是,它是独立于平台的,无论是在X86上跑,还是在PowerPC上跑,你什么都不用担心,编译器会为你搞定一切。而这也恰恰是它最大的缺点——严重依赖于编译器的实现。而实际测试表明,编译器所生成的代码,其速度实在不尽人意。

一个替代的方法是直接对数据位进行操作。如果你对IEEE浮点数的表示法比较熟悉的话(如果你压根什么都不知道,请先查阅文末附录中的资料),这是显而易见的。它提取指数和尾数,然后对尾数执行移位操作。代码如下:
----------------------------------------------------------------------------------------
//
// 将32位浮点数fval转换为32位整数并存储在ival中
// 小数部分将被裁剪掉
//
void TruncToInt32 (int &ival, float fval)
{
ival = *(int *)&fval;

// 提取尾数
// 注意实际的尾数前面还有一个被省略掉的1
int mantissa = (ival & 0x07fffff) | 0x800000;

// 提取指数
// 以23分界,指数大于23则左移,否则右移
// 由于指数用偏移表示,所以23+127=150
int exponent = 150 - ((ival >> 23) & 0xff);

if (exponent < 0)
ival = (mantissa << -exponent);
else
ival = (mantissa >> exponent);

// 如果小于0,则将结果取反
if ((*(int *)&fval) & 0x80000000)
ival = -ival;
}
----------------------------------------------------------------------------------------
该函数有一个BUG,那就是当fval=0时,返回值是2。原因是对于语句mantissa>>exponent,编译器使用了循环移位指令。解决方法是要么对0作特殊处理,要么直接用汇编来实现。

这个函数比标准的C转换要快,而且由于整个过程只使用了整数运算,可以和FPU并行运行。但缺点是,(1)依赖于硬件平台。例如根据小尾和大尾顺序的不同,要相应地修改函数。(2)对于float和double要使用不同的实现,因为二者的数据位不同。(3)对于float,只能保留24位有效值,尽管int有32位。

更快的方法是使用FPU指令FISTP,它将栈中的浮点数弹出并保存为整数:
----------------------------------------------------------------------------------------
//
// 将64位浮点数fval转换为32位整数并存储在ival中
// 小数部分将四舍五入到偶数
//
inline void RoundToIntFPU (int &ival, double fval)
{
_asm
{
fld fval
mov edx, dword ptr [ival]
fistp dword ptr [edx]
}
}
----------------------------------------------------------------------------------------
很好,无论速度还是精度似乎都相当令人满意。但如果换一个角度来看的话,fistp指令需要6个cycle,而浮点数乘法才仅仅需要3个cycle!更糟的是,当fistp运行的时候,它必须占用FPU,也就是说,其他的浮点运算将不能执行。仅仅为了一次类型转换操作就要付出如此大的代价,光想想就觉得心疼。

当然,它也有很多优点:更快的速度,更精确的数值(四舍五入到偶数),更强的适用性(float和double均可)。要注意的是,FPU默认的四舍五入到偶数(round to even)和我们平常所说的四舍五入(round)是不同的。二者的区别在于对中间值的处理上。考虑十进制的3.5和4.5。四舍五入到偶数是使其趋向于邻近的偶数,所以舍入的结果是3.5->4,4.5->4;而传统的四舍五入则是3.5->4,4.5->5。四舍五入到偶数可以产生更精确,更稳定的数值。

除此之外,还有更好,更快的方法吗?有的,那就是华丽的 Magic Number !

请看下面的函数:
----------------------------------------------------------------------------------------
//
// 将64位浮点数转换为32位整数
// 小数部分将四舍五入到偶数
//
inline void RoundToInt64 (int &val, double dval)
{
static double magic = 6755399441055744.0;
dval += magic;
val = *(int*)&dval;
}
----------------------------------------------------------------------------------------
快,这就是它最伟大的地方!你所需要的仅仅是一次浮点数加法,你还能再奢望些什么呢?

当然,绝对不要使用莫名其妙的代码,现在马上就让我们来看看它是怎么变的戏法。

首先,我们要搞清楚FPU是怎样执行浮点数加法的。考虑一下8.75加2^23。8.75的二进制表示是1000.11,转化为IEEE标准浮点数格式是1.00011*2^3。假设二者均是32位的float,它们在内存中的存储为:
----------------------------------------------------------------------------------------
符号位(31), 指数(30-23), 尾数(22-0)

8.75: 0, 10000010, 0001 1000 0000 0000 0000 000

2^23: 0, 10010110,0000 0000 0000 0000 0000 000
----------------------------------------------------------------------------------------
FPU所执行的操作是:(1)提升指数较小者的指数,使二者指数相同;(2)将二者的尾数相加;(3)将结果规整为标准格式。让我们具体来看看整个过程:
----------------------------------------------------------------------------------------
1,提升8.75的指数,尾数要相应地右移23-3=20位:

8.75 = 1.00011*2^3 = .0000000000000000000100011*2^23

2,将二者相加。必须注意的是,
   由于float只有23位尾数,所以8.75的低位被移除掉了:

8.75: 0, 10010110, 0000 0000 0000 0000 0001 000
+
2^23: 0, 10010110,0000 0000 0000 0000 0000 000
=
0, 10010110, 0000 0000 0000 0000 0001 000

3,将规整为标准格式:

别忘了2^23还有个前导1,所以结果是规整的,无需执行任何操作
----------------------------------------------------------------------------------------
你发现什么了吗?不错,将结果的尾数部分提取出来,正是int(8.75)!是的,magic number的奥妙就在这里,通过迫使FPU将尾数移位,从而获得我们需要的结果。

但是别高兴得太早,让我们来看看负数的情况。以-8.75为例,-8.75加2^23相当于2^23减去8.75,过程如下:
----------------------------------------------------------------------------------------
2^23: 0, 10010110,0000 0000 0000 0000 0000 000
-
8.75: 0, 10010110, 0000 0000 0000 0000 0001 000
=
0, 10010110, 1111 1111 1111 1111 1110 000
----------------------------------------------------------------------------------------
很好,尾数部分正是int(-8.75)=-8的补码。然而,2^23的前导1在减法过程中被借位了,所以结果的尾数前面并没有1,FPU将执行规整操作,最后我们得到的是:
----------------------------------------------------------------------------------------
0, 10010110, 1111 1111 1111 1111 1100 000
----------------------------------------------------------------------------------------
功亏一篑!等等,如果将2^23的尾数的最高位置1,那么在减法过程中不就可以保全前导1了吗?完全正确,这就是我们需要的。所以最后的magic number是0,10010110,1000 0000 0000 0000 0000 000,也就是1.5*2^23。

然而,我们要清楚这个方法的一些局限性。首先,在将结果的float值保存为整数时,我们需要使用某些mask值将22-31位屏蔽掉。其次,float只能保有最多23位的有效值,在将尾数最高位置1后,有效值更是降到了22位,这意味着我们对大于2^23-1的数值无能为力。

相比之下,如果我们只处理double,那么所有的问题就都迎刃而解了。首先,double的指数位,符号位和尾数的最高位全部都在高32位,无需屏蔽操作。其次,double有多达52位的尾数,对于32位的int来说,实在是太奢侈了!

用于double的magic number是1.5*2^52=6755399441055744.0,推导过程是一样的。

根据同样的原理,我们还可以计算出将float和double转化为定点数的magic number。例如对于16-16的定点数,尾数右移的位数比原先转化为整数时要少16,所以对于double来说,相应的magic number就是1.5*2^36。如果要转化为8-24的定点数,则移位次数要少24,magic number是1.5*2^28。对于其他格式的定点数,以此类推。比起int(float_val*65536)这样的表达式,无论是速度还是精度都要优胜得多。

另外,不得不再次重申的是,对于在最后的结果中被移除掉的数值,FPU会将其四舍五入到偶数。然而,有时我们确实需要像floor和ceil那样的操作,那该怎么办呢?很简单,将原数加上或减去一个修正值就可以了,如下所示:
----------------------------------------------------------------------------------------
// 修正值
static double magic_delta=0.499999999999;

// 截取到整数
inline void Floor64 (int &val, double dval)
{
RoundToInt64(val,dval-delta);
}

// 进位到整数
inline void Ceil64 (int &val, double dval)
{
RoundToInt64(val,dval+delta);
}
----------------------------------------------------------------------------------------
这个世界上没有免费的午餐。我们获得了速度,相对应地,也必须付出一些代价,那就是可移植性。上述方法全都基于以下几个假设:(1)在x86上跑;(2) 符合IEEE的浮点数标准;(3)int为32位,float为32位,double为64位。局限性其实是蛮大的,相比之下,int_val= (int)float_val就要省事多了。当然,你也可以使用条件编译。

最后,让我们来看两组实际的测试数值。这份报告来自于Sree Kotay和他的朋友Herf:
----------------------------------------------------------------------------------------
平台:Pentium IV/1.2

64位浮点数到32位整数的转换:

int(f):        2772.65 ms
fistp:         679.269 ms
magic number: 622.519 ms

64位浮点数到16-16定点数的转换:

int(f*65536): 2974.57 ms
fistp:         3100.84 ms
magic number: 606.80 ms

floor函数:

ANSI floor:    7719.400 ms
magic number: 687.177 ms
----------------------------------------------------------------------------------------
平台:Pentium D/3.2

64位浮点数到32位整数的转换:

int(f):        1948.470 ms
fistp:         523.792 ms
magic number: 333.033 ms

64位浮点数到16-16定点数的转换:

int(f*65536): 2163.56 ms
fistp:         7103.66 ms
magic number: 335.03 ms

floor函数:

ANSI floor:    3771.55 ms
magic number: 380.25 ms
----------------------------------------------------------------------------------------
性能的差距实在令人惊讶,不是吗?我想说的是,写编译器的家伙绝对不是傻瓜(恐怕比你我都要聪明得多),他们当然知道所有这些优秀的算法。但为什么编译器的表现会如此糟糕呢?其中一个理由是,为了使浮点数运算尽可能精确和迅速,IEEE在算法的选择上有很多的限制。而另一方面,IEEE的舍入规则(四舍五入到偶数)尽管从维持浮点数的连贯性上来看非常棒,但却不符合ANSI C在浮点数到整数的类型转换上的说明(截尾)。于是,编译器不得不做一大堆的工作来保证行为的正确性(符合标准)。这在大部分情况下都不是什么问题——除非你在写图形/声音/多媒体之类的代码。

刚刚在我的赛扬1.1G上实际测试了一下。编译器是VC2003,使用PerformenceCounter来计算时间,在DEBUG模式下,C标准转换/FISTP/MagicNumber三种方法所耗费时间的比约为5/3/2。但在优化全开的RELEASE模式下,标准C类型转换的速度却是遥遥领先于所有其他的方法!也不知道是我的测试方法有问题,还是该赞VS厉害。
--------------------------------------------------------------------------------------
参考文献:
1,What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg(PDF):
这篇论文囊括了关于浮点数的所有东西,正如其名字所示,任何想要了解浮点数的人都必读的文献。(其中关于精度的讨论实在令我受益非浅。)

2,Let's Get to The Floating Point by Chris Hecker(PDF):
早期的关于浮点数的文章,写得非常棒,值得一看。

3,IEEE Standard 754 for Binary Floating Point Arithmetic by Prof.W.Kahan(PDF):
关于IEEE754标准的说明。

4,IA64 Floating Point Operations and the IEEE Standard for Binary Floating Point Arithmetic(PDF):
关于IA64的浮点数实现。

5,Know Your FPU: Fixing Floating Fast by Sree Kot