¿Está disponible la precisión arbitraria de punto flotante?
Solo por diversión y porque fue realmente fácil, escribí un programa corto para generar números de injerto , pero debido a problemas de precisión de punto flotante no encuentro 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
A este código le falta 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 escribí un programa corto en C++ para ver si era mi CPU el que truncaba el número de punto flotante o Python 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;
}
Qué salidas:
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
Entonces parece que estoy corriendo contra los límites de la precisión de punto flotante y la CPU está cortando los bits restantes porque piensa que la diferencia restante es un error de punto flotante. ¿Hay alguna manera de solucionar este problema en Python? ¿O necesito pasar a C y usar GMP o algo así?
En la biblioteca estándar, el decimalmódulo puede ser lo que estás buscando. Además, mpmath me ha resultado bastante útil. La documentación también tiene muchos ejemplos excelentes (desafortunadamente, la computadora de mi oficina no la tiene mpmathinstalada; de lo contrario, verificaría algunos ejemplos y los publicaría).
Sin embargo, una advertencia sobre el decimalmódulo. El módulo contiene varias funciones integradas para operaciones matemáticas simples (por ejemplo sqrt), pero es posible que los resultados de estas funciones no siempre coincidan con la función correspondiente en mathu otros módulos con precisiones más altas (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, esto genera las dos primeras líneas
math.sqrt: 0.37796447300922719758631274089566431939601898193359375
decimal.sqrt: 0.377964473009227227214516536234
actual value: 0.3779644730092272272145165362341800608157513118689214
lo cual, como se indicó, no es exactamente lo que cabría esperar y puede ver que cuanto mayor es la precisión, menos coinciden los resultados. Tenga en cuenta que el decimalmódulo tiene más precisión en este ejemplo, ya que coincide más con el valor real .
Para este problema en particular, decimales una excelente manera de hacerlo, porque 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)
Dado que estás buscando una propiedad que se exprese de forma más natural en notación decimal, es un poco tonto usar una representación binaria. La página de Wikipedia a la que se vinculó no indicó cuántos "dígitos no injertados" pueden aparecer antes de que comiencen los "dígitos injertados", por lo que esto 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 muchas posibilidades de que el resultado de Decimal.sqrt()sea 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 lo 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')
Puedes probar con decimal en lugar de punto flotante.
