如何計算反向觸發（和sqrt）函數（C語言）的浮點運算中的舍入錯誤？

Question

我有一個相當複雜的函數，它需要幾個double值，它們表示以緯度和經度爲弧度的形式（大小，緯度，經度）的3維空間中的兩個向量以及一個角度。函數的作用是將第一個向量繞第二個向量旋轉指定的角度並返回合成向量。我已經驗證該代碼在邏輯上是正確的並且可行。

函數的預期用途是用於圖形，所以不需要雙精度;然而，在目標平臺上，使用浮點數（sinf，cosf，atan2f，asinf，acosf和sqrtf）的trig（和sqrt）函數在double上的運行速度要快於浮點數（可能是因爲計算這些值的指令實際上可能需要double;如果傳遞一個float值，則該值必須轉換爲double值，這要求將其複製到具有更多內存的區域 - 即開銷）。因此，函數中涉及的所有變量都是雙精度的。

這是問題：我試圖優化我的功能，以便它可以被稱爲每秒更多次。因此，我將這些調用替換爲sin，cos，sqrt等調用這些函數的浮點版本，因爲它們導致總體速度提高3-4倍。這適用於幾乎所有的輸入;然而，如果輸入矢量接近於與標準單位矢量（i，j或k）平行，則各種函數的舍入誤差會足以導致稍後調用sqrtf或反函數trig函數（asinf，acosf， atan2f）在這些函數的域之外傳遞參數僅僅是。因此，我留下了這個困境：要麼我只能調用雙精度函數，並避免這個問題（並且每秒最多可以有大約1,300,000個向量操作），或者我可以試着想出一些東西其他。最終，我想要一種方法來消除對反向trig函數的輸入以處理邊緣情況（對於sqrt，這樣做很簡單：只需使用abs）。分支不是一種選擇，因爲即使是單個條件語句也會增加很多開銷，以至於性能增益都會丟失。

那麼，有什麼想法？

編輯：有人對我使用雙打和浮點操作表示困惑。如果我實際上將所有值存儲在雙倍大小的容器中（I.E. double-type變量），則將函數存儲在浮動大小容器中的速度要快得多。但是，由於顯而易見的原因，浮點精度觸發操作比雙精度觸發操作要快。

Answer 1

基本上，你需要找到一個解決你的問題的numerically stable算法。這種事情沒有通用的解決方案，需要根據您的具體情況完成，如果單個步驟使用的概念如condition number。如果潛在的問題本身是病態的，事實上可能是不可能的。

Answer 2

你看過Graphics Programming Black Book還是將計算交給GPU？

Answer 3

單精度浮點固有地引入錯誤。所以，你需要通過使用epsilon因子來建立你的數學運算，以便所有的比較都有一定程度的「漏洞」，並且你需要清理輸入到有限域的函數。

前者轉移時，如

bool IsAlmostEqual(float a, float b) { return fabs(a-b) < 0.001f; } // or 
bool IsAlmostEqual(float a, float b) { return fabs(a-b) < (a * 0.0001f); } // for relative error 

是很容易的，但是這是混亂的。限制域輸入有點棘手，但更好。關鍵是要使用conditional move operators，這在一般做這樣的事情

float ExampleOfConditionalMoveIntrinsic(float comparand, float a, float b) 
{ return comparand >= 0.0f ? a : b ; } 

在單個運算

，而不會產生一個分支。

這些取決於架構而變化。在x87浮點單元上，你可以用FCMOV conditional-move op來完成，但這很笨拙，因爲它取決於之前設置的條件標誌，所以速度很慢。另外，cmov沒有一致的編譯器。這就是爲什麼我們在可能的情況下避免x87浮點數支持SSE2標量數學的原因之一。

條件此舉是更好的配對與位與一個comparison operator在SSE支持。這是即使是標量運算可取：

// assuming you've already used _mm_load_ss to load your floats onto registers 
__m128 fsel(__m128 comparand, __m128 a, __m128 b) 
{ 
    __m128 zero = {0,0,0,0}; 
    // set low word of mask to all 1s if comparand > 0 
    __m128 mask = _mm_cmpgt_ss(comparand, zero); 
    a = _mm_and_ss(a, mask); // a = a & mask 
    b = _mm_andnot_ss(mask, b); // b = ~mask & b 
    return _mm_or_ss(a, b);  // return a | b 
    } 
} 

編譯器是更好的，但不是很大，約在啓用SSE2標量運算髮射這種用於ternaries格局。您可以使用MSVC上的編譯器標記/arch:sse2或GCC上的-mfpmath=sse來完成此操作。

在PowerPC和許多其他的RISC體系結構，fsel()是硬件操作碼和因此通常編譯器本徵爲好。