在計算機圖形運算中,常常要將浮點數轉換為整數,例如在圖像的光柵化階段,就要執行大量的類型轉換,以便將浮點數表示的坐標轉化為整數表示的屏幕坐標。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
posted on 2008-01-31 22:23
小果子 閱讀(7185)
評論(6) 編輯 收藏 引用 所屬分類:
學習筆記