¿Precisión arbitraria de punto flotante de Python disponible?

Question

Solo por diversión y porque era realmente fácil, he escrito un programa corto para generar Grafting numbers, pero debido a problemas de precisión de punto flotante no está encontrando algunos de los ejemplos más grandes.

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 

Este código pierde al menos un número de injerto conocido. 9999999998 => 99999.99998999999999949999999994999999999374999999912... Parece perder precisión adicional después de multiplicar por 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 

Así que escribió un pequeño programa C++ para ver si era mi CPU truncar el número de coma flotante o pitón de alguna manera.

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

que da salida:

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 

para que parezca que estoy corriendo con fuerza contra los límites de precisión de punto flotante y la CPU está cortando el resto de bits, ya que piensa que la diferencia restante es un error de coma flotante ¿Hay alguna forma de solucionar esto bajo Python? ¿O necesito pasar a C y usar GMP o algo?

Answer 1

En la biblioteca estándar, el módulo decimal puede ser lo que está buscando. Además, he encontrado mpmath para ser bastante útil. El documentation tiene muchos buenos ejemplos también (desafortunadamente mi computadora de oficina no tiene mpmath instalado; de lo contrario, verificaría algunos ejemplos y los publicaría).

Una advertencia sobre el módulo decimal, sin embargo. El módulo contiene varias funciones incorporadas para operaciones matemáticas simples (por ejemplo, sqrt), pero los resultados de estas funciones pueden no coincidir siempre con la función correspondiente en math u otros módulos con mayores precisiones (aunque pueden ser más precisos). Por ejemplo,

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

En Python 3.2.3, estas salidas de las dos primeras líneas

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

la que como se ha dicho, no es exactamente lo que cabría esperar, y se puede ver que a mayor precisión, menos coinciden los resultados. Tenga en cuenta que el módulo decimal tiene más precisión en este ejemplo, ya que se acerca más al actual value.

Answer 2

Puede intentarlo con Decimal en lugar de floatingpoint.

Answer 3

Python no tiene flotantes de precisión arbitraria incorporados, pero hay paquetes de Python de terceros que usan GMP: gmpy y PyGMP.

Answer 4

Para este problema en particular, decimal es una excelente forma de hacerlo, ¡ya que almacena los dígitos decimales como tuplas!

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

Puesto que usted está buscando una propiedad que se expresa de forma natural en la mayoría notación decimal, que es un poco tonto para usar una representación binaria.La página de Wikipedia que vincula a no indicó cuántos pueden aparecer "dígitos no injerto" ante los "dígitos de injerto" comienzan, por lo que este le permite especificar:

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

Creo que hay una buena probabilidad de que el resultado de Decimal.sqrt() será más preciso, al menos para esto, que el resultado de math.sqrt() debido a la conversión entre representación binaria y representación decimal. Considere la siguiente, por ejemplo:

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