解決浮點舍入問題C++

我開發了一個科學應用程序（模擬在細胞核中移動的染色體）。染色體被分成小片段，使用4x4旋轉矩陣圍繞隨機軸旋轉。解決浮點舍入問題C++

問題是模擬執行了數千億次的旋轉，因此浮點舍入誤差會堆積並呈指數級增長，所以碎片會隨着時間流逝而「飄離」並與染色體的其餘部分分離。

我使用C++的雙精度。目前在CPU上軟運行，但將被移植到CUDA，並且仿真最多可持續1個月。因爲所有的片段都被鏈接在一起（你可以看到它是一個雙向鏈表），但我認爲如果可能的話，這將是最好的想法。

你有什麼建議嗎？我感覺有點失落。

非常感謝你，

編輯：增加了一個簡單的示例代碼。你可以假設所有的矩陣數學都是經典的實現。

// Rotate 1000000 times 
for (int i = 0; i < 1000000; ++i) 
{ 
    // Pick a random section start 
    int istart = rand() % chromosome->length; 

    // Pick the end 20 segments further (cyclic) 
    int iend = (istart + 20) % chromosome->length; 

    // Build rotation axis 
    Vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position; 
    axis.normalize(); 

    // Build rotation matrix and translation vector 
    Matrix4 rotm(axis, rand()/float(RAND_MAX)); 
    Vector4 oldpos = chromosome->segments[istart].position; 

    // Rotate each segment between istart and iend using rotm 
    for (int j = (istart + 1) % chromosome->length; j != iend; ++j, j %= chromosome->length) 
    { 
     chromosome->segments[j].position -= oldpos; 
     chromosome->segments[j].position.transform(rotm); 
     chromosome->segments[j].position += oldpos; 
    } 
}

來源

2011-04-11 Robert Troipartrois

數值分析和穩定性是一個巨大的領域。沒有一個正確的答案。沒有看到一些示例代碼，很難給出任何具體的建議。 – 2011-04-11 22:04:09

你說得對，我添加了一些代碼，如果這可能有幫助。 – 2011-04-11 22:17:26

順便說一句，這聽起來像一個很酷的項目。 – 2011-04-11 22:26:26

您需要爲您的系統找到一些約束，並努力將其保持在合理範圍內。我做了大量的分子碰撞模擬，並且在這些系統中總能量是守恆的，所以每一步我都要仔細檢查系統的總能量，如果它的變化達到某個閾值，那麼我知道我的時間步長選擇的很差（太大或太小），我選擇一個新的時間步並重新運行。這樣我可以實時跟蹤系統發生的情況。

對於這個模擬，我不知道你有多少守恆量，但是如果你有一個守恆量，你可以試着保持不變。請記住，縮短時間步長並不總是會提高準確度，您需要根據精確度來優化步長。我已經進行了數週的CPU時間數值模擬，並且守恆數量總是在10^8的1個部分內，所以有可能，你只需要玩一些。另外，正如Tomalak所說，也許試着總是引用你的系統到開始時間，而不是上一步。因此，不是總是移動你的染色體，而是將染色體保存在其起始位置，並與它們一起存儲一個轉換矩陣，從而使你到達當前位置。當你計算你的新旋轉時，只需修改變換矩陣。它可能看起來很愚蠢，但有時這很好，因爲錯誤的平均值爲0.

例如，假設我有一個位於（x，y）處的粒子和我計算的每個步驟（dx，dy）和移動粒子。分步的方式將做到這一點

t0 (x0,y0) 
t1 (x0,y0) + (dx,dy) -> (x1, y1) 
t2 (x1,y1) + (dx,dy) -> (x2, y2) 
t3 (x2,y2) + (dx,dy) -> (x3, y3) 
t4 (x3,30) + (dx,dy) -> (x4, y4) 
...

如果你總是參考T0，你可以在任何時候，TN做到這一點

t0 (x0, y0) (0, 0) 
t1 (x0, y0) (0, 0) + (dx, dy) -> (x0, y0) (dx1, dy1) 
t2 (x0, y0) (dx1, dy1) + (dx, dy) -> (x0, y0) (dx2, dy2) 
t3 (x0, y0) (dx2, dy2) + (dx, dy) -> (x0, y0) (dx3, dy3)

所以，要想讓你的，你需要做的實際位置（ x0，y0）+（dxn，dyn）

現在爲我的例子簡單的翻譯，你可能不會贏得很多。但是對於輪換而言，這可以成爲一種生活救星。只需保留與每個染色體相關聯的歐拉角的矩陣並更新該矩陣，而不是染色體的實際位置。至少這樣他們不會飄走。

來源

2011-04-11 23:05:55 miked

+1用於檢查能量。從確保數值收斂到特定解決方案的角度出發，我會建議。（但是，我會阻止參考長時間模擬的開始時間，因爲浮點時間值將失去精度並可能失速。） – Potatoswatter 2011-04-11 23:23:19

節約能源是一個很好的方法，至少要注意您的時間系統出問題了。如何糾正能量的收益/損失當然可能具有挑戰性。這也不是完美的，因爲系統的一部分可能會獲得，而另一部分會損失，等於0. – 2011-04-12 03:59:05

只是爲了澄清，在這樣的系統中，能量並不守恆，因爲它不是閉合的。相反，熵是最大化的：每個模擬步驟應該具有相對較高的降低總能量的概率，然後隨機波動將溫度恢復到正常。 – Potatoswatter 2011-04-12 07:59:21

寫您的公式，這樣的時間步T數據不僅僅從浮點數據的時間步長T-1派生。儘量確保浮點錯誤的產生僅限於一個時間步。

在沒有更具體的問題解決的情況下，很難說出更具體的內容。

來源

2011-04-11 22:04:16

謝謝，不幸的是，模擬的本質不是*連續*，而是完全依賴於兩個時間步長。 – 2011-04-11 22:15:44

@Heandel：我從你的代碼中看到了這個問題（並且記得幾次讓它感到沮喪）。除了使用一些'bignum'庫來最大限度地提高浮點類型的精度之外，我不太清楚你可以做些什麼。雖然我可能會錯過一些明顯的東西;我很久沒有做圖形數學了。 – 2011-04-11 22:27:32

我想過使用'bignum'或'apfloat'，但是我還沒有找到任何CUDA的實現！ – 2011-04-12 08:35:36

我認爲這取決於你使用的編譯器。

的Visual Studio編譯器支持/ fp的開關，它告訴浮點運算

可以read more about it的行爲。基本上，/ fp：嚴格是最苛刻的模式

來源

2011-04-11 22:06:12 Eric

謝謝，我使用的是gcc，並嘗試了所有可能的想法，但沒有任何區別。 MSVC使用/ fp：strict和/ fp：precise也給了我相同的結果。 – 2011-04-11 22:16:32

我想這取決於所需的精度，但您可以使用基於'整數'的浮點數。通過這種方法，您可以使用整數併爲小數位數提供您自己的偏移量。

例如，與4個小數點精確，你將有

浮點值 - > int值 1.0000 - > 10000 1.0001 - > 10001 0.9999 - > 09999

你必須當你進行乘法和除法時要小心，當你應用精度偏移時要小心。其他方面，你可以很快得到溢出錯誤。

1.0001 * 1.0001變爲10001 *一萬分之一萬〇一

來源

2011-04-11 22:09:31 yamspog

問題描述相當含糊，所以這裏有一些比較模糊的建議。

選項1：

尋找一些約束集，從而（1）應始終持有，（2），如果他們失敗了，但只是剛剛，可以很容易地調整系統，使他們這樣做，（3 ）如果他們全都成立，那麼你的模擬不會太瘋狂，以及（4）當系統開始瘋狂時，約束開始失敗，但只是輕微。例如，或許染色體相鄰位之間的距離應該至多爲d，對於某些d，並且如果少數距離略大於d，則可以（例如）從一端沿着染色體行走，修復通過將下一個片段及其所有後繼片段移動到其前一個片段的距離太大。或者其他的東西。

然後經常檢查約束條件以確保任何違規行爲在被捕獲時仍然很小;當你發現違規行爲時，請糾正。（你應該安排，當你修理東西時，你「不僅僅是滿足」約束條件）。

如果檢查限制的時間很長，那麼當然你可以做到這一點。（這樣做可以使你做的修正更便宜，例如，如果這意味着任何違反總是很小的。）

選項2：

查找描述，使系統的狀態的一個新途徑這個問題不可能出現。例如，也許（我懷疑這一點），你可以爲每個相鄰的碎片對存儲一個旋轉矩陣，並強制它總是一個正交矩陣，然後讓這些碎片的位置由這些旋轉矩陣隱式確定。

方案3：

而不是你的約束思維的限制，提供一些小的「復原力」，這樣，當某樣東西脫節也容易被拉向後朝事情應該是這樣。請注意，沒有任何問題時，恢復力爲零或至少可以忽略不計，以免它們比原始數字錯誤更嚴重地干擾您的結果。

來源

2011-04-11 22:14:18

你比較模糊的選項3看起來非常有前途。謝謝！ – 2011-04-12 08:15:21

如果我正確讀取此代碼，任何時候任何兩個相鄰染色體片段之間的距離都應該改變。在這種情況下，在主循環計算每對相鄰點之間和主循環之後的距離之前，必要時移動每個點，以便與前一點具有適當的距離。

根據具體情況，您可能需要在主循環中多次強制執行此約束。

來源

2011-04-11 23:44:54 bmcnett

基本上，您需要避免來自這些（不精確的）矩陣運算符的錯誤累積，並且在大多數應用程序中有兩種主要方法。

而不是寫的位置上多次操作的一些初始位置，你可以寫出來，究竟該運營商將個運算後明確。例如，假設你有一個位置x，並且你正在增加一個值e（你不能完全表示）。比計算x + = e好得多;大量的時間將是計算x + EN;其中EN是一些更準確的方式來表示N次後的操作。你應該想一下，你是否有更精確地表示許多旋轉動作的方法。
稍微有些人爲的是把你的新發現的點和項目從你的旋轉中心預期的半徑差異。這將保證它不會漂移（但不一定能保證旋轉角度準確。）

來源

2011-04-12 03:33:46

解決浮點舍入問題C++

回答

相關問題