Precisione arbitraria in virgola mobile Python disponibile?

Question

Solo per divertimento e perché è stato davvero facile, ho scritto un breve programma per generare Grafting numbers, ma a causa di problemi di precisione in virgola mobile non è trovare alcuni degli esempi più grandi.

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 

Questo codice manca almeno un numero di innesto conosciuto. 9999999998 => 99999.99998999999999949999999994999999999374999999912... Sembra diminuire la precisione extra dopo averlo moltiplicato per 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 

così ho scritto un breve programma C++ per vedere se era la mia CPU troncare il numero in virgola mobile o pitone in qualche modo.

#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; 
} 

quali uscite:

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 

così sembra che io sto correndo duro contro i limiti della precisione in virgola mobile e la CPU è tagliando i restanti bit perché pensa che la differenza residua è un errore in virgola mobile. C'è un modo per aggirare questo sotto Python? O devo passare a C e usare GMP o qualcosa del genere?

Answer 1

Nella libreria standard, il modulo decimal potrebbe essere quello che stai cercando. Inoltre, ho trovato mpmath essere abbastanza utile. Il documentation ha anche molti ottimi esempi (purtroppo il mio computer dell'ufficio non ha installato mpmath, altrimenti verificherei alcuni esempi e li pubblicheremo).

Un avvertimento sul modulo decimal, tuttavia. Il modulo contiene diverse funzioni integrate per semplici operazioni matematiche (ad esempio sqrt), ma i risultati di queste funzioni potrebbero non corrispondere sempre alla funzione corrispondente in math o altri moduli con precisioni più elevate (sebbene possano essere più precisi). Ad esempio,

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())) 

In Python 3.2.3, questo emette le prime due righe

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

che, come detto, non è esattamente quello che ci si aspetta, e si può vedere che maggiore è la precisione, meno risultati corrispondono. Notare che il modulo decimal ha maggiore precisione in questo esempio, poiché corrisponde più strettamente allo actual value.

Answer 2

Puoi provare con Decimal invece di virgola mobile.

Answer 3

Python non ha galleggianti di precisione arbitraria incorporati, ma ci sono pacchetti Python di terze parti che utilizzano GMP: gmpy e PyGMP.

Answer 4

Per questo problema particolare, decimal è un ottimo modo per andare, perché memorizza le cifre decimali come tuple!

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

Dal momento che siete alla ricerca di una proprietà che è più naturalmente espresso in notazione decimale, è un po 'stupido per usare una rappresentazione binaria.La pagina di Wikipedia si è collegato al non ha indicato quanti "cifre non-innesto" possono apparire prima dei "innesto cifre" iniziano, quindi questo permette di specificare:

>>> 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 

Penso che ci sia una buona probabilità il risultato di Decimal.sqrt() sarà più preciso, almeno per questo, del risultato di math.sqrt() a causa della conversione tra rappresentazione binaria e rappresentazione decimale. Si consideri il seguente, ad esempio:

>>> 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')