浮点数到整数的快速转换

浮点数到整数的快速转换

 

     在计算机图形运算中，常常要将浮点数转换为整数，例如在图像的光栅化阶段，就要执行大量的类型转换，以便将浮点数表示的坐标转化为整数表示的屏幕坐标。

 

 ----------------------------------------------------------------------------------------

 //

 // 强制类型转换

 // 小数部分将被裁剪掉

 //

 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在浮点数到整数的类型转换上的说明（截尾）。于是，编译器不得不做一大堆的工作来保证行为的正确性（符合标准）。这在大部分情况下都不是什么问题——除非你在写图形/声音/多媒体之类的代码。