浮點數的美麗與哀愁

這幾年個人在影像處理程式優化的領域打滾, 如果問到感到棘手的工作, floating point 的處理應該可以排上很前面的名次

在許多演算來說由於同時對於 precision 與 dynamic range 的需求, 因此在計算過程中對於浮點數的使用是非常常見的 (若要避免使用會有很高的專業與困難度), floating 主要優點在於可以表示極大與極小值, 相較整數能大幅避免 overflow 與 underflow, 缺點是有效位數的減少, 而且現今多數的計算單元都俱備 floating 的支援, 已經讓一些人疏於了使用浮點的問題, (包含486與之前的時代 FPU是高檔貨, ARM 也自 ARMv7 才列標配)

然而若橫跨了 PC 與 手機, CPU 與 GPU, CPU 與 DSP, 甚至於三者 ( PC, 手機, GPU), floating point 就變成非常難以考量與處理的負擔, 而為了區分是程式錯誤或是誤差就必須耗費相當的心力

為了簡化問題, 因此文中談到 floating point 若無指名, 一律是指 32bit single precision, 但相同的問題 64bit double precision 中一樣存在

IEEE 754

首先必須要談的是問題的核心 - IEEE 754

對於計算機而言, 浮點數是以上圖的格式存放
fraction 一共 23 bits 存放一個介於 1~2 之間的數目, 這 b22 ~ b0 存放的是二進位小數以下的部分, 也就是說 fraction 所表示的數值為:

1 + b22*(2^-1) + b21*(2^-2) + ... + b0*(2^-23)

而 exponent 代表著指數, 一樣以 2 的只數次方表示, 具有 8bit, 因此值域為 0~255, 但是預設會減去 127, 所以即為 -127 ~ 128, 因此對於 exponent 本身所表示的數值為:

2^(exponent - 127)

而 sign 就不用多談了, 這是用以表示正負

然而 IEEE 754 定義的不僅僅只是 format 而已, 還有著 rounding mode, required operations 以及 exception handling, 符合了 IEEE 754 的規範下, 才可能有相同的輸出結果 (這當然只是一個最低門檻)

Format 本身的問題

扣除浮點數因格式問題不可能表示全部的數目外
格式本身最大的問題是因為 dynamic range 的移動, 像是 (A+B)+C 與 A+(B+C) 單以代數考量這無疑是相等的, 但是若以 floating 格式去思考, 你就會意會到輸出結果很有可能會不同, 而原始演算實作所採用的累加或相乘的順序, 必須在優化實作上努力維持才能產生一致的輸出結果, 這樣的問題對於程式優化影響最大的部分是平行計算, 無論 TLP or ILP, 因為平行度優化考量, 都會有分割與不同面向個別累計的需求, 如此勢必都會產生一定的誤差

CPU 間的問題

或許有些人會認為只使用 CPU 那麼就不會碰到浮點數問題了, 這樣說只算對了一半, 而且你還是必須要只使用一種平台以及指令集, 對於多數演算法設計工程師而言, 他們很慣於使用 PC 平台, 甚至會使用 MATLAB, x86 PC 上的程式預設會使用 x87 浮點數協同單元 指令集, 而這是許多問題的開始 - x87 內部使用 80bit 浮點數表示

通常 x86 CPU 是在 x87 FPU 的內部以 80b 計算得到結果後, 再 truncate 為 IEEE-754 的 float(32b) or double(64b), 這就表示這與 IEEE 754 FPU 的輸出結果會有微小的差異, 目前常見的手機的 ARM 架構, 即為 IEEE 754 compatible FPU, 所以光是 ARM 與 x86 PC 相同的程式碼其輸出結果基本上就會有所不同, 而對於 ARM 與 x86 的部分, 就必須仰賴以 IEEE 754 設計的 SSE2 指令集, 若是 gcc 與 clang 很早就可透過 -mfpmath=sse2 的編譯參數來達到, 但是 Visual Studio 必須是 2013 版後才有正確的 code generation 實作, (也就是說 Windows 的使用者要安裝 VS 2012以後的版本 才有可能透過 /arch:sse2 有一致的輸出)

對於 ARM 與 x86 平台的一致性方案反而又揭露了另一個層面問題:

就算使用了單一CPU架構, 在 ISA 指令集間的支援還是會有所不一致!

類似於 x86 平台上 x87 指令與 SSE2 指令有著不同輸出, 同樣地 ARM VFP 指令(IEEE 754 相容)與 NEON 指令(非完全 IEEE 754 相容) 也可能會有輸出結果不同, 而這樣的問題還會再帶到 libmath 的實作方式, 讓要處理一致性的問題再度的變得更嚴重

GPU

GPU 本身有著龐大的浮點數計算能力, 但是通常為了能達到更高的吞吐量以及加速上的考量, 在計算結果與 IEEE-754 可能存在差異, 不同代的 GPU 或是不同架構都有可能有所不同. 像是 CUDA 是在 compute compatibilty v2.0 之後才完備了 IEEE 754 的支援, 除此之外許多硬體加速的數學函數的輸出上也不保證與 CPU 一致, 這點 Nvidia 在 2011 GTC 中給的 Floating Point and IEEE-754 Compliance for Nvidia GPUs 簡報中有很詳盡的說明. 對於其他GPU 以及各CPU/GPU 平台上的 OpenCL 中的 built-in functions 的實作與支援也有著相同的道理.

DSP

對於 DSP 而言這樣的痛苦並不在於 IEEE 754 本身, 而是多數的多媒體面向的 DSP 為了考量計算能力與面積, 結果多半是直接不俱備 floating 能力的, 像是 Qualcomm S82x 中的 Hexagon 680 HVX 就不俱備 floating 運算的 SIMD 指令, 而通常的處理作法是採用 fixed-point (quantization) 的浮點模擬, 然而若採用靜態位數的方式容易失真, 而動態的方式有著實作上的複雜度以及多餘計算的負擔. 而數學函數上的實作若難以避免則通常必須透過相當紆迴的方式.

Lookup-Table or Frame-based Parameters

對於跨裝置的正確性驗證, 由單一裝置輸出的 Lookup Table(單一的 math function 像是 sin, cos, log, exp 等等) 或是一整張預先透過單一裝置計算的 Frame-based Parameters(複雜的並結合多個 math function 的運算) , 是常用來確認誤差單純是由 floating 計算造成的技巧. 以此來確保實作上的流程與邏輯無誤.

延伸閱讀: The pitfalls of verifying floating-point computations