浮点数到整数的转换

在计算机图形运算中，常常要将浮点数转换为整数，例如在图像的光栅化阶段，就要执行大量的类型转换，以便将浮点数表示的坐标转化为整数表示的屏幕坐标。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厉害。
--------------------------------------------------------------------------------------
参考文献和相关资源可到鄙人的小店下载：http://greatsorcerer.go2.icpcn.com/info/float2int.html

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 Kotay