Raíz cuadrada entera en python

¿Hay una raíz cuadrada entera en algún lugar de python, o en bibliotecas estándar? Quiero que sea exacto (es decir, devuelva un número entero), y ladrar si no hay solución.

En este momento yo lancé mi propia ingenua:

def isqrt(n): i = int(math.sqrt(n) + 0.5) if i**2 == n: return i raise ValueError('input was not a perfect square') 

Pero es feo y realmente no confío en él para los enteros grandes. Podría recorrer los cuadrados y rendirme si hubiera excedido el valor, pero asumo que sería un poco lento hacer algo así. También creo que probablemente estaría reinventando la rueda, algo como esto seguramente ya debe existir en Python

El método de Newton funciona perfectamente bien en enteros:

 def isqrt(n): x = n y = (x + 1) // 2 while y < x: x = y y = (x + n // x) // 2 return x 

Esto devuelve el entero más grande x para el cual x * x no excede n . Si desea verificar si el resultado es exactamente la raíz cuadrada, simplemente realice la multiplicación para verificar si n es un cuadrado perfecto.

Analizo este algoritmo y otros tres algoritmos para calcular raíces cuadradas en mi blog .

Lo siento por la respuesta muy tardía; Me topé con esta página. En caso de que alguien visite esta página en el futuro, el módulo de python gmpy2 está diseñado para trabajar con entradas muy grandes e incluye, entre otras cosas, una función de raíz cuadrada entera.

Ejemplo:

 >>> import gmpy2 >>> gmpy2.isqrt((10**100+1)**2) mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001L) >>> gmpy2.isqrt((10**100+1)**2 - 1) mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L) 

Por supuesto, todo tendrá la etiqueta “mpz”, pero los mpz son compatibles con int:

 >>> gmpy2.mpz(3)*4 mpz(12) >>> int(gmpy2.mpz(12)) 12 

Vea mi otra respuesta para una discusión del rendimiento de este método en relación con algunas otras respuestas a esta pregunta.

Descargar: https://code.google.com/p/gmpy/

Aquí hay una implementación muy sencilla:

 def i_sqrt(n): i = n.bit_length() >> 1 # i = floor( (1 + floor(log_2(n))) / 2 ) m = 1 << i # m = 2^i # # Fact: (2^(i + 1))^2 > n, so m has at least as many bits # as the floor of the square root of n. # # Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2) # >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED. # while m*m > n: m >>= 1 i -= 1 for k in xrange(i-1, -1, -1): x = m | (1 << k) if x*x <= n: m = x return m 

Esto es sólo una búsqueda binaria. Inicialice el valor m para que sea la potencia más grande de 2 que no exceda la raíz cuadrada, luego verifique si se puede establecer cada bit más pequeño manteniendo el resultado no más grande que la raíz cuadrada. (Verifique los bits uno a la vez, en orden descendente.)

Para valores razonablemente grandes de n (por ejemplo, alrededor de 10**6000 , o alrededor de 20000 bits), esto parece ser:

  • Más rápido que la implementación del método de Newton descrita por user448810 .
  • Mucho, mucho más lento que el gmpy2 incorporado gmpy2 en mi otra respuesta .
  • Comparable a, pero algo más lento que la raíz cuadrada de Longhand descrita por nibot .

Todos estos enfoques tienen éxito en las entradas de este tamaño, pero en mi máquina, esta función toma alrededor de 1.5 segundos, mientras que @ Nibot toma alrededor de 0.9 segundos, @ user448810's toma alrededor de 19 segundos, y el método incorporado gmpy2 toma menos de un milisegundo. (!). Ejemplo:

 >>> import random >>> import timeit >>> import gmpy2 >>> r = random.getrandbits >>> t = timeit.timeit >>> t('i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # This function 1.5102493192883117 >>> t('exact_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # Nibot 0.8952787937686366 >>> t('isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # user448810 19.326695976676184 >>> t('gmpy2.isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # gmpy2 0.0003599147067689046 >>> all(i_sqrt(n)==isqrt(n)==exact_sqrt(n)[0]==int(gmpy2.isqrt(n)) for n in (r(1500) for i in xrange(1500))) True 

Esta función se puede generalizar fácilmente, aunque no es tan buena porque no tengo una estimación inicial tan precisa para m :

 def i_root(num, root, report_exactness = True): i = num.bit_length() / root m = 1 << i while m ** root < num: m <<= 1 i += 1 while m ** root > num: m >>= 1 i -= 1 for k in xrange(i-1, -1, -1): x = m | (1 << k) if x ** root <= num: m = x if report_exactness: return m, m ** root == num return m 

Sin embargo, tenga en cuenta que gmpy2 también tiene un método i_root .

De hecho, este método podría adaptarse y aplicarse a cualquier función f (no negativa, creciente) para determinar un "entero inverso de f ". Sin embargo, para elegir un valor inicial eficiente de m , aún querría saber algo sobre f .

Edición: Gracias a @Greggo por señalar que la función i_sqrt se puede reescribir para evitar el uso de multiplicaciones. ¡Esto produce un aumento de rendimiento impresionante!

 def improved_i_sqrt(n): assert n >= 0 if n == 0: return 0 i = n.bit_length() >> 1 # i = floor( (1 + floor(log_2(n))) / 2 ) m = 1 << i # m = 2^i # # Fact: (2^(i + 1))^2 > n, so m has at least as many bits # as the floor of the square root of n. # # Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2) # >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED. # while (m << i) > n: # (m<>= 1 i -= 1 d = n - (m << i) # d = nm^2 for k in xrange(i-1, -1, -1): j = 1 << k new_diff = d - (((m<<1) | j) << k) # n-(m+2^k)^2 = nm^2-2*m*2^k-2^(2k) if new_diff >= 0: d = new_diff m |= j return m 

Tenga en cuenta que, por construcción, el k ½ bit de m << 1 no está establecido, por lo tanto en modo bit a bit o se puede usar para implementar la adición de (m<<1) + (1< . En última instancia, tengo (2*m*(2**k) + 2**(2*k)) escrito como (((m<<1) | (1< , por lo que son tres turnos y uno en modo bit-o (seguido de una resta para obtener new_diff ). Tal vez todavía hay una manera más eficiente de conseguir esto? En cualquier caso, es mucho mejor que multiplicar m*m ! Comparar con lo anterior:

 >>> t('improved_i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. 0.10908999762373242 >>> all(improved_i_sqrt(n) == i_sqrt(n) for n in xrange(10**6)) True 

Algoritmo de raíz cuadrada de mano larga

Resulta que hay un algoritmo para calcular las raíces cuadradas que puedes calcular a mano, algo así como la división larga. Cada iteración del algoritmo produce exactamente un dígito de la raíz cuadrada resultante y consume dos dígitos del número cuya raíz cuadrada busca. Si bien la versión de “mano larga” del algoritmo se especifica en decimal, funciona en cualquier base, siendo el binario más simple de implementar y quizás el más rápido de ejecutar (según la representación bignum subyacente).

Debido a que este algoritmo opera con números dígito a dígito, produce resultados exactos para cuadrados perfectos arbitrariamente grandes, y para cuadrados no perfectos, puede producir tantos dígitos de precisión (a la derecha del lugar decimal) como se desee.

Hay dos reseñas interesantes en el sitio de “Dr. Math” que explican el algoritmo:

  • Raíces cuadradas en binario
  • Raíces Cuadradas A Mano

Y aquí hay una implementación en Python:

 def exact_sqrt(x): """Calculate the square root of an arbitrarily large integer. The result of exact_sqrt(x) is a tuple (a, r) such that a**2 + r = x, where a is the largest integer such that a**2 <= x, and r is the "remainder". If x is a perfect square, then r will be zero. The algorithm used is the "long-hand square root" algorithm, as described at http://mathforum.org/library/drmath/view/52656.html Tobin Fricke 2014-04-23 Max Planck Institute for Gravitational Physics Hannover, Germany """ N = 0 # Problem so far a = 0 # Solution so far # We'll process the number two bits at a time, starting at the MSB L = x.bit_length() L += (L % 2) # Round up to the next even number for i in xrange(L, -1, -1): # Get the next group of two bits n = (x >> (2*i)) & 0b11 # Check whether we can reduce the remainder if ((N - a*a) << 2) + n >= (a<<2) + 1: b = 1 else: b = 0 a = (a << 1) | b # Concatenate the next bit of the solution N = (N << 2) | n # Concatenate the next bit of the problem return (a, Na*a) 

Podría modificar fácilmente esta función para realizar iteraciones adicionales para calcular la parte fraccionaria de la raíz cuadrada. Estaba más interesado en calcular raíces de cuadrados grandes y perfectos.

No estoy seguro de cómo esto se compara con el algoritmo del "método entero de Newton". Sospecho que el método de Newton es más rápido, ya que en principio puede generar múltiples bits de la solución en una iteración, mientras que el algoritmo de "mano larga" genera exactamente un bit de la solución por iteración.

Fuente de repo: https://gist.github.com/tobin/11233492

Una opción sería usar el módulo decimal y hacerlo en flotadores suficientemente precisos:

 import decimal def isqrt(n): nd = decimal.Decimal(n) with decimal.localcontext() as ctx: ctx.prec = n.bit_length() i = int(nd.sqrt()) if i**2 != n: raise ValueError('input was not a perfect square') return i 

que creo que debería funcionar:

 >>> isqrt(1) 1 >>> isqrt(7**14) == 7**7 True >>> isqrt(11**1000) == 11**500 True >>> isqrt(11**1000+1) Traceback (most recent call last): File "", line 1, in  isqrt(11**1000+1) File "", line 10, in isqrt raise ValueError('input was not a perfect square') ValueError: input was not a perfect square 

Aquí comparé todas las funciones (correctas) tanto en entradas pequeñas (0… 2 22 ) como grandes (2 50001 ). Los ganadores claros en ambos casos son gmpy2.isqrt sugerido por mathmandan en primer lugar, seguido por la receta ActiveState vinculada por NPE en segundo lugar. La receta de ActiveState tiene un montón de divisiones que se pueden reemplazar por turnos, lo que la hace un poco más rápida (pero aún detrás de gmpy2.isqrt ):

 def isqrt(n): if n > 0: x = 1 << (n.bit_length() + 1 >> 1) while True: y = (x + n // x) >> 1 if y >= x: return x x = y elif n == 0: return 0 else: raise ValueError("square root not defined for negative numbers") 

Resultados de referencia:

  • gmpy2.isqrt() (mathmandan) : 0.08 µs pequeño, 0.07 ms grande
  • int(gmpy2.isqrt()) *: 0.3 µs pequeño, 0.07 ms grande
  • ActiveState (optimizado como se indica arriba) : 0.6 µs pequeño, 17.0 ms grande
  • ActiveState (npe) : 1.0 µs pequeño, 17.3 ms grande
  • Castlebravo de mano larga : 4 µs pequeña, 80 ms grande
  • Mathmandan mejorado : 2.7 µs pequeño, 120 ms grande
  • martineau (con esta corrección ): 2.3 µs pequeño, 140 ms grande
  • nibot : 8 µs pequeño, 1000 ms grande
  • Mathmandan : 1.8 µs pequeño, 2200 ms grande
  • Método de Castlebravo Newton : 1.5 µs pequeño, 19000 ms grande
  • user448810 : 1.4 µs pequeño, 20000 ms grande

(* Dado que gmpy2.isqrt devuelve un objeto gmpy2.mpz , que se comporta principalmente pero no exactamente como un int , es posible que deba volver a convertirlo en un int para algunos usos).

Parece que podrías comprobar así:

 if int(math.sqrt(n))**2 == n: print n, 'is a perfect square' 

Actualizar:

Como señalaste, lo anterior falla para valores grandes de n . Para ellos, el siguiente parece prometedor, que es una adaptación del código C de ejemplo, por Martin Guy @ UKC, junio de 1985, para el método de cálculo de números binarios número por dígito de aspecto relativamente simple mencionado en el artículo de Wikipedia Métodos de cálculo de raíces cuadradas :

 from math import ceil, log def isqrt(n): res = 0 bit = 4**int(ceil(log(n, 4))) if n else 0 # smallest power of 4 >= the argument while bit: if n >= res + bit: n -= res + bit res = (res >> 1) + bit else: res >>= 1 bit >>= 2 return res if __name__ == '__main__': from math import sqrt # for comparison purposes for i in range(17)+[2**53, (10**100+1)**2]: is_perfect_sq = isqrt(i)**2 == i print '{:21,d}: math.sqrt={:12,.7G}, isqrt={:10,d} {}'.format( i, sqrt(i), isqrt(i), '(perfect square)' if is_perfect_sq else '') 

Salida:

  0: math.sqrt= 0, isqrt= 0 (perfect square) 1: math.sqrt= 1, isqrt= 1 (perfect square) 2: math.sqrt= 1.414214, isqrt= 1 3: math.sqrt= 1.732051, isqrt= 1 4: math.sqrt= 2, isqrt= 2 (perfect square) 5: math.sqrt= 2.236068, isqrt= 2 6: math.sqrt= 2.44949, isqrt= 2 7: math.sqrt= 2.645751, isqrt= 2 8: math.sqrt= 2.828427, isqrt= 2 9: math.sqrt= 3, isqrt= 3 (perfect square) 10: math.sqrt= 3.162278, isqrt= 3 11: math.sqrt= 3.316625, isqrt= 3 12: math.sqrt= 3.464102, isqrt= 3 13: math.sqrt= 3.605551, isqrt= 3 14: math.sqrt= 3.741657, isqrt= 3 15: math.sqrt= 3.872983, isqrt= 3 16: math.sqrt= 4, isqrt= 4 (perfect square) 9,007,199,254,740,992: math.sqrt=9.490627E+07, isqrt=94,906,265 100,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,020,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001: math.sqrt= 1E+100, isqrt=10,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001 (perfect square) 

Su función falla para entradas grandes:

 In [26]: isqrt((10**100+1)**2) ValueError: input was not a perfect square 

Hay una receta en el sitio de ActiveState que, con suerte, debería ser más confiable, ya que utiliza solo números enteros de matemáticas. Se basa en una pregunta anterior de StackOverflow: escribir su propia función de raíz cuadrada

Los flotadores no se pueden representar con precisión en las computadoras. Puede probar la configuración de proximidad deseada de épsilon en un valor pequeño dentro de la precisión de los flotadores de python.

 def isqrt(n): epsilon = .00000000001 i = int(n**.5 + 0.5) if abs(i**2 - n) < epsilon: return i raise ValueError('input was not a perfect square') 

He comparado los diferentes métodos dados aquí con un bucle:

 for i in range (1000000): # 700 msec r=int(123456781234567**0.5+0.5) if r**2==123456781234567:rr=r else:rr=-1 

encontrando que este es el más rápido y no necesita importación matemática. Mucho tiempo puede fallar, pero mira esto

 15241576832799734552675677489**0.5 = 123456781234567.0 

Pruebe esta condición (sin computación adicional):

 def isqrt(n): i = math.sqrt(n) if i != int(i): raise ValueError('input was not a perfect square') return i 

Si necesita devolver un int (no un float con un cero final), asigne una segunda variable o compute int(i) dos veces.