Python float無聲溢出精度錯誤

Question

StackOverflow在浮點表示方面有很多主題，關於異常，截斷和精度問題。我一直在試圖解釋這一點，但仍然沒有弄清楚。

from operator import add, sub, mul, div 

fun = 'add(add(sub(sub(safeDiv(xL1, add(D, mul(sub(add(sub(A, 0), mul(D, rv)), mul(xL1, add(D, 3))), xL1))), add(0, rv)), safeDiv(mul(sub(D, sigma2), safeDiv(sub(safeDiv(xL1, A), 1), add(safeDiv(safeDiv(B1, xL1), sub(4, xL2)), add(sigma1, xL1)))), sigma1)), add(4, B1)), add(add(A, A), sub(add(xL1, xL1), mul(xL2, safeDiv(xL1, add(sub(add(mul(D, -4), add(add(safeDiv(mul(sigma2, sigma2), safeDiv(B1, sigma1)), sub(add(D, safeDiv(xL2, B1)), D)), sub(4, B1))), A), add(mul(sigma2, xL1), mul(xL1, mul(rv, xL2)))))))))' 

d = [(
51.696521954140991, 
31.156112806482234, 
54.629863633907163, 
27.491618858013698, 
26.223584534107289, 
77.10005191617563, 
2708.4145268939428, 
0.20952943771134946, 
15.558278150405643, 
102.0, 
225.0)] 

arglabels = ['xL1', 'sigma1', 'xL2', 'sigma2', 'A', 'B1', 'D', 'rv'] 
other = {'add': add, 'sub':sub, 'mul':mul,'safeDiv':div} 
inputs = dict(zip(arglabels, d[0][: -4] + (d[0][-3]*d[0][-4],))) 
inputs.update(other) 
print eval(fun, inputs) 

此代碼應產生一個介於225和240之間的結果，而是返回一個負數。就是這樣，也不例外，沒有警告，什麼也沒有。所以它必定是一個精確的錯誤，將結果完全關閉。

通過四捨五入最大我可以用來得到一個合理的結果是1位小數（這使我接近207 ...），在某些情況下從numpy幫助longdoubles，但還不夠。我已經完成了它（如此相當精確的損失，並獲得240）。

另一個細節，在筆記本與主腳本運行起來我有this behaviour：

當我添加當地人字典首次返回一個非常合理的結果，但接下來的時間它回到了負值。必須有一些導入影響這一點，但我也找不到它。

我該怎麼做才能避免這種情況？我怎樣才能得到某種警告？我怎樣才能跟蹤哪裏出錯？

編輯：接受的答案正確地確定了問題，看看答案下面的評論瞭解更多詳情。但是，它沒有討論如何避免它或糾正函數。也許這應該是MathOverflow討論...

Answer 1

如果預期的結果是在範圍225...240下列問題是可能的：

• 在d
• 功能add，sub，mul和safeDiv代fun
• 不正確的值的過程中的誤差應該做一些比浮點加法，減法，乘法和除法更復雜的事情。

在問題中提供的輸入不能給出除-2786.17215265之外的任何內容，因爲它是一個完美的數字結果。沒有浮點舍入錯誤或溢出。在測試腳本的輸出結果下面顯示了算術函數的詳細版本，所有浮點運算都已定義好。沒有什麼可以給出明顯的舍入誤差。當扣除接近值時存在一些有風險的操作：

add: -10626.8589858 + 10627.794501 = 0.935515251547 

但是，這遠不是舍入誤差。

通過C程序或數學工具（MATLAB/Octave）可以獲得相同的結果。

---

在屏幕截圖的不同的輸出是由於未顯示有本地變量的值。由於Out[108]與Out[110]相同，我假定ds與data[17]相同。局部變量用於輸出Out[109]和Out[110]，因此差異值爲rv，因爲只有該變量在In[110]中更改。如果所有其他變量的值是固定的，則可以看出，如果rv等於以下值之一，則可以獲得Out[109]（230.62977145178198）：4.075164147,4.485709922，51.72476610。以下測試腳本中的最後一個輸出行也說明了這一點。

注意，如果fun被分析爲rv的函數，它有兩個極點（約3.3和51.7）。所以，從技術上講，函數可以給出從負無窮到正無窮的結果。

---

測試

from operator import add, sub, mul, div 

fun = 'add(add(sub(sub(safeDiv(xL1, add(D, mul(sub(add(sub(A, 0), mul(D, rv)), mul(xL1, add(D, 3))), xL1))), add(0, rv)), safeDiv(mul(sub(D, sigma2), safeDiv(sub(safeDiv(xL1, A), 1), add(safeDiv(safeDiv(B1, xL1), sub(4, xL2)), add(sigma1, xL1)))), sigma1)), add(4, B1)), add(add(A, A), sub(add(xL1, xL1), mul(xL2, safeDiv(xL1, add(sub(add(mul(D, -4), add(add(safeDiv(mul(sigma2, sigma2), safeDiv(B1, sigma1)), sub(add(D, safeDiv(xL2, B1)), D)), sub(4, B1))), A), add(mul(sigma2, xL1), mul(xL1, mul(rv, xL2)))))))))' 

d = [(
51.696521954140991, 
31.156112806482234, 
54.629863633907163, 
27.491618858013698, 
26.223584534107289, 
77.10005191617563, 
2708.4145268939428, 
0.20952943771134946, 
15.558278150405643, 
102.0, 
225.0)] 

def add_verbose(a,b): 
    res = a + b; 
    print "add: {0} + {1} = {2}".format(a,b,res); 
    return res; 

def sub_verbose(a,b): 
    res = a - b; 
    print "sub: {0} - {1} = {2}".format(a,b,res); 
    return res; 

def div_verbose(a,b): 
    res = a/b; 
    print "div: {0}/{1} = {2}".format(a,b,res); 
    return res; 

def mul_verbose(a,b): 
    res = a * b; 
    print "mul: {0} * {1} = {2}".format(a,b,res); 
    return res; 

arglabels = ['xL1', 'sigma1', 'xL2', 'sigma2', 'A', 'B1', 'D', 'rv'] 
other = {'add': add_verbose, 'sub':sub_verbose, 'mul':mul_verbose,'safeDiv':div_verbose} 
inputs = dict(zip(arglabels, d[0][: -4] + (d[0][-3]*d[0][-4],))) 
inputs.update(other) 

# out 108 
print eval(fun, inputs) 

# set locals that can give out 109 
safeDiv = div; 
rv = 4.0751641470166795256; 
inputs.update(locals()); 

# out 109 
print eval(fun, inputs) 

輸出

sub: 26.2235845341 - 0 = 26.2235845341 
mul: 2708.41452689 * 3.25991727261 = 8829.20729761 
add: 26.2235845341 + 8829.20729761 = 8855.43088215 
add: 2708.41452689 + 3 = 2711.41452689 
mul: 51.6965219541 * 2711.41452689 = 140170.700616 
sub: 8855.43088215 - 140170.700616 = -131315.269734 
mul: -131315.269734 * 51.6965219541 = -6788542.72473 
add: 2708.41452689 + -6788542.72473 = -6785834.3102 
div: 51.6965219541/-6785834.3102 = -7.61830006318e-06 
add: 0 + 3.25991727261 = 3.25991727261 
sub: -7.61830006318e-06 - 3.25991727261 = -3.25992489091 
sub: 2708.41452689 - 27.491618858 = 2680.92290804 
div: 51.6965219541/26.2235845341 = 1.97137511414 
sub: 1.97137511414 - 1 = 0.971375114142 
div: 77.1000519162/51.6965219541 = 1.49139727397 
sub: 4 - 54.6298636339 = -50.6298636339 
div: 1.49139727397/-50.6298636339 = -0.029456869265 
add: 31.1561128065 + 51.6965219541 = 82.8526347606 
add: -0.029456869265 + 82.8526347606 = 82.8231778914 
div: 0.971375114142/82.8231778914 = 0.0117283004453 
mul: 2680.92290804 * 0.0117283004453 = 31.4426693361 
div: 31.4426693361/31.1561128065 = 1.00919744165 
sub: -3.25992489091 - 1.00919744165 = -4.26912233256 
add: 4 + 77.1000519162 = 81.1000519162 
add: -4.26912233256 + 81.1000519162 = 76.8309295836 
add: 26.2235845341 + 26.2235845341 = 52.4471690682 
add: 51.6965219541 + 51.6965219541 = 103.393043908 
mul: 2708.41452689 * -4 = -10833.6581076 
mul: 27.491618858 * 27.491618858 = 755.789107434 
div: 77.1000519162/31.1561128065 = 2.47463643475 
div: 755.789107434/2.47463643475 = 305.414200171 
div: 54.6298636339/77.1000519162 = 0.708558065477 
add: 2708.41452689 + 0.708558065477 = 2709.12308496 
sub: 2709.12308496 - 2708.41452689 = 0.708558065477 
add: 305.414200171 + 0.708558065477 = 306.122758237 
sub: 4 - 77.1000519162 = -73.1000519162 
add: 306.122758237 + -73.1000519162 = 233.02270632 
add: -10833.6581076 + 233.02270632 = -10600.6354013 
sub: -10600.6354013 - 26.2235845341 = -10626.8589858 
mul: 27.491618858 * 51.6965219541 = 1421.22107785 
mul: 3.25991727261 * 54.6298636339 = 178.088836061 
mul: 51.6965219541 * 178.088836061 = 9206.57342319 
add: 1421.22107785 + 9206.57342319 = 10627.794501 
add: -10626.8589858 + 10627.794501 = 0.935515251547 
div: 51.6965219541/0.935515251547 = 55.2599456488 
mul: 54.6298636339 * 55.2599456488 = 3018.84329521 
sub: 103.393043908 - 3018.84329521 = -2915.4502513 
add: 52.4471690682 + -2915.4502513 = -2863.00308224 
add: 76.8309295836 + -2863.00308224 = -2786.17215265 
-2786.17215265 
230.629771452 

Answer 2

當我嘗試你的榜樣，我得到的答案：

>>> print eval(fun, inputs) 
-2786.17215265 

如果我使用gmpy2和設置精度到200位，指數範圍〜1E9，我得到了一個答案：

>>> print eval(fun,inputs) 
-2786.1721526580839894614784542009831125135156833413835128962432 

它看起來像函數返回一個穩定的結果。所以這個函數可能有問題。

我會遵循@ Prune的建議，並將複雜的功能分成更小的步驟。