2012-07-17 88 views
30

只是爲了好玩,因爲它非常簡單,我寫了一個簡短的程序來生成Grafting numbers,但由於浮點精度問題,它沒有找到一些較大的例子。Python浮點任意精度可用?

def isGrafting(a): 
    for i in xrange(1, int(ceil(log10(a))) + 2): 
    if a == floor((sqrt(a) * 10**(i-1)) % 10**int(ceil(log10(a)))): 
     return 1 

a = 0 
while(1): 
    if (isGrafting(a)): 
    print "%d %.15f" % (a, sqrt(a)) 
    a += 1 

該代碼遺漏了至少一個已知嫁接編號。 9999999998 => 99999.99998999999999949999999994999999999374999999912...它似乎降低了乘以10**5後額外的精度。

>>> a = 9999999998 
>>> sqrt(a) 
99999.99999 
>>> a == floor((sqrt(a) * 10**(5)) % 10**int(ceil(log10(a)))) 
False 
>>> floor((sqrt(a) * 10**(5)) % 10**int(ceil(log10(a)))) 
9999999999.0 
>>> print "%.15f" % sqrt(a) 
99999.999989999996615 
>>> print "%.15f" % (sqrt(a) * 10**5) 
9999999999.000000000000000 

因此,我寫了一個簡短的C++程序,以查看它是否是我的CPU以某種方式截斷了浮點數或python。

#include <cstdio> 
#include <cmath> 
#include <stdint.h> 

int main() 
{ 
    uint64_t a = 9999999998; 
    printf("%ld %.15f %.15f %.15f %.15f\n", a, sqrt((double)a), sqrt((double)a)*1e4, sqrt((double)a)*1e5, sqrt((double)a)*1e6); 
    a = 999999999998; 
    printf("%ld %.15f %.15f %.15f %.15f\n", a, sqrt((double)a), sqrt((double)a)*1e5, sqrt((double)a)*1e6, sqrt((double)a)*1e7); 
    a = 99999999999998; 
    printf("%ld %.15f %.15f %.15f %.15f\n", a, sqrt((double)a), sqrt((double)a)*1e6, sqrt((double)a)*1e7, sqrt((double)a)*1e8); 
    return 0; 
} 

,輸出:

9999999998 99999.999989999996615 999999999.899999976158142 9999999999.000000000000000 99999999990.000000000000000 
999999999998 999999.999998999992386 99999999999.899993896484375 999999999999.000000000000000 9999999999990.000000000000000 
99999999999998 9999999.999999899417162 9999999999999.900390625000000 99999999999999.000000000000000 999999999999990.000000000000000 

因此,它看起來像我跑起來難防的浮點精度的極限,因爲它認爲該CPU被斬去其餘位,其餘的差異是浮點錯誤。有沒有辦法在Python下解決這個問題?或者我需要遷移到C並使用GMP或其他?

+1

可以使用['fractions'模塊](https://docs.python.org/3/library/fractions.html)執行有理數的精確算術。 – jfs 2015-09-30 14:55:47

回答

31

在標準庫中,decimal模塊可能是您要查找的內容。另外,我發現mpmath是相當有幫助的。 documentation也有很多很好的例子(不幸的是我的辦公室電腦沒有安裝mpmath;否則我會驗證一些例子併發布它們)。但是,有關decimal模塊的一個警告。該模塊包含用於簡單數學運算的幾個內置函數(例如sqrt),但這些函數的結果可能並不總是與更高精度的math或其他模塊中的相應函數相匹配(儘管它們可能更準確)。例如,

from decimal import * 
import math 

getcontext().prec = 30 
num = Decimal(1)/Decimal(7) 

print(" math.sqrt: {0}".format(Decimal(math.sqrt(num)))) 
print("decimal.sqrt: {0}".format(num.sqrt())) 

在Python 3.2.3,這個輸出前兩行

math.sqrt: 0.37796447300922719758631274089566431939601898193359375 
decimal.sqrt: 0.377964473009227227214516536234 
actual value: 0.3779644730092272272145165362341800608157513118689214 

其作爲說明,不正是你所期望的,你可以看到,在高精度越高,結果越少。請注意,decimal模塊在此示例中具有更高的準確性,因爲它更接近actual value

+2

+1 mpmath。使用十進制數字的問題在於,您無法在十進制對象的數學函數中做很多事情,所以如果您只是在玩耍,那麼它非常有限。 – DSM 2012-07-17 13:08:13

+0

@DSM我同意。我只用'mpmath'來解決相當簡單的問題,但是我發現它是一個很好的,可靠的軟件包。 – 2012-07-17 13:10:40

+1

只是要清楚 - 我敢肯定,在你的'math.sqrt'與'Decimal.sqrt()'的測試中,'math.sqrt'產生的結果是_less_正確的,因爲binary-to - 轉換。考慮'decimal.Decimal(math.sqrt(num)** 2)* 7'的輸出與'decimal.Decimal(num.sqrt()** 2)* 7'的輸出。 – senderle 2012-07-17 13:30:58

7

你可以用Decimal而不是浮點數來嘗試。

5

Python沒有內置的任意精度浮點數,但有第三方Python包使用GMP:gmpyPyGMP

8

對於這個特定的問題,decimal是一個很好的方法,因爲它將十進制數字存儲爲元組!

>>> a = decimal.Decimal(9999999998) 
>>> a.as_tuple() 
DecimalTuple(sign=0, digits=(9, 9, 9, 9, 9, 9, 9, 9, 9, 8), exponent=0) 

既然你要找的是十進制的最自然的表達的屬性,這是一個有點傻使用二進制表示。維基百科頁面鏈接到您沒有說明有多少「嫁接數字」開始之前「非嫁接數字」可能出現,所以這讓你指定:

>>> def isGrafting(dec, max_offset=5): 
...  dec_digits = dec.as_tuple().digits 
...  sqrt_digits = dec.sqrt().as_tuple().digits 
...  windows = [sqrt_digits[o:o + len(dec_digits)] for o in range(max_offset)] 
...  return dec_digits in windows 
... 
>>> isGrafting(decimal.Decimal(9999999998)) 
True 
>>> isGrafting(decimal.Decimal(77)) 
True 

我覺得有一個很好的機會的結果由於二進制表示和十進制表示之間的轉換,因此至少對於此,Decimal.sqrt()將比math.sqrt()的結果更準確。考慮以下內容,例如:

>>> num = decimal.Decimal(1)/decimal.Decimal(7) 
>>> decimal.Decimal(math.sqrt(num) ** 2) * 7 
Decimal('0.9999999999999997501998194593') 
>>> decimal.Decimal(num.sqrt() ** 2) * 7 
Decimal('1.000000000000000000000000000')