Raíz cuadrada entera

Los siguientes programas en Python son implementaciones sencillas.

def isqrt ( y : int ) -> int : """  Raíz cuadrada entera  (búsqueda lineal, ascendente)  """ # subestimación inicial, L <= isqrt(y) L = 0 while ( L + 1 ) * ( L + 1 ) <= y : L += 1regresar L

def isqrt ( y : int ) -> int : """  Raíz cuadrada entera  (búsqueda lineal, descendente)  """ # sobreestimación inicial, isqrt(y) <= R R = y while ( R * R > y ): R -= 1devolver R

En el programa anterior (búsqueda lineal, ascendente) se puede reemplazar la multiplicación por la suma, utilizando la equivalencia $${\displaystyle (L+1)^{2}=L^{2}+2L+1=L^{2}+1+\sum _ {i=1}^{L}2.}$$

def isqrt ( y : int ) -> int : """  Raíz cuadrada entera  (búsqueda lineal, ascendente) usando suma  """ L = 0 a = 1 d = 3mientras a <= y : a = a + d d = d + 2 L = L + 1regresar L

La búsqueda lineal comprueba secuencialmente cada valor hasta que encuentra el más pequeño.${\displaystyle x}$dónde${\displaystyle x^{2}>y}$.

Se consigue una mayor velocidad utilizando la búsqueda binaria .

def isqrt ( y : int ) -> int : """ Raíz cuadrada entera (búsqueda binaria) """ L = 0 # límite inferior de la raíz cuadrada R = y + 1 # límite superior de la raíz cuadradamientras ( L != R - 1 ): M = ( L + R ) // 2 # punto medio para probar si ( M * M <= y ): L = M sino : R = Mregresar L

Ejemplos numéricos

• ${\displaystyle \operatorname {isqrt} (0)=0}$,${\displaystyle \operatorname {isqrt} (1)=1}$.
• Utilizando la búsqueda binaria, el cálculo de${\displaystyle \operatorname {isqrt} (131072)}$converge a${\displaystyle 362}$en${\displaystyle 17}$pasos de iteración a través del${\displaystyle [L,R]}$secuencia

${\displaystyle {\begin{aligned}&[0,131073]\to [0,65536]\to [0,32768]\to [0,16384]\to [0,8192]\to [0,4096]\rightarrow [0,2048]\to [0,1024]\to [0,512]\\&\to [256,512]\to [256,384]\to [320,384]\to [352,384]\to [352,368]\to [360,368]\to [360,364]\to [362,364]\to [362,363]\end{aligned}}}$

• El cálculo de${\displaystyle \operatorname {isqrt} (2000000)}$converge a${\displaystyle 1414}$en${\displaystyle 21}$pasos a través del${\displaystyle [L,R]}$secuencia

${\displaystyle {\begin{aligned}&[0,2000001]\to [0,1000000]\to [0,500000]\to [0,250000]\to [0,125000]\to [0,62500]\to [0,31250]\to [0,15625]\\&\to [0,7812]\to [0,3906]\to [0,1953]\to [976,1953]\to [976,1464]\to [1220,1464]\to [1342,1464]\to [1403,1464]\\&\to [1403,1433]\to [1403,1418]\to [1410,1418]\to [1414,1418]\to [1414,1416]\to [1414,1415]\end{aligned}}}$

Búsqueda lineal (ascendente, comenzando desde${\displaystyle 0}$) necesidades1414 escalones.

Para computación${\displaystyle \lfloor {\sqrt {n}}\rfloor }$Se puede utilizar el cociente de la división euclidiana para ambas operaciones de división. Esto tiene la ventaja de usar solo números enteros para cada valor intermedio, lo que hace innecesario el uso de representaciones de punto flotante .

Dejar${\displaystyle n>0}$y suposición inicial${\displaystyle x_{0}>0}$. Defina la secuencia de enteros:

${\displaystyle x_{k+1}=\left\lfloor {\frac {x_{k}+\left\lfloor n/x_{k}\right\rfloor }{2}}\right\rfloor ,\quad k=0,1,2,\dots }$

1. Positividad: Todos los términos son números enteros positivos:${\displaystyle x_{k}>0}$a pesar de${\displaystyle k}$.

2. Monotonía:

• Si${\displaystyle x_{k}>{\sqrt {n}}}$, entonces${\displaystyle \lfloor n/x_{k}\rfloor \leq n/x_{k}}$;

entonces${\displaystyle x_{k+1}=\left\lfloor {\frac {x_{k}+\lfloor n/x_{k}\rfloor }{2}}\right\rfloor <{\frac {x_{k}+n/x_{k}}{2}}<x_{k}}$.

Por lo tanto, la secuencia disminuye.

• Si${\displaystyle x_{k}<{\sqrt {n}}}$, entonces${\displaystyle \lfloor n/x_{k}\rfloor \geq n/x_{k}-1}$;

entonces${\displaystyle x_{k+1}\geq {\frac {x_{k}+n/x_{k}-1}{2}}>x_{k}-1}$.

Por lo tanto, la secuencia aumenta o permanece igual.

3. Acotación: La sucesión está acotada inferiormente por 1 y superiormente por${\displaystyle x_{0}}$, por lo tanto, está delimitado.

4. Estabilización / Oscilación: Una secuencia de enteros monótona y acotada se estabiliza u oscila entre dos enteros consecutivos:

${\displaystyle x_{k+1}=x_{k}}$o${\displaystyle x_{k+1}=x_{k}\pm 1}$.

5. Condición de "punto fijo" entero: En la estabilización o la oscilación:

${\displaystyle x_{k+1}=\lfloor (x_{k}+\lfloor n/x_{k}\rfloor )/2\rfloor }$.

Esto garantiza que la secuencia esté en${\displaystyle \lfloor {\sqrt {n}}\rfloor }$o oscilando entre los dos enteros más cercanos alrededor${\displaystyle {\sqrt {n}}}$.

6. Conclusión: La secuencia finalmente se estabiliza en${\displaystyle \lfloor {\sqrt {n}}\rfloor }$o oscila entre${\displaystyle \lfloor {\sqrt {n}}\rfloor }$y${\displaystyle \lceil {\sqrt {n}}\rceil }$.

Observación:

• ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$es un punto fijo estricto a menos que${\displaystyle n+1}$es un cuadrado perfecto.
• Si${\displaystyle n+1}$es un cuadrado perfecto, la secuencia oscila entre${\displaystyle \lfloor {\sqrt {n}}\rfloor }$y${\displaystyle \lceil {\sqrt {n}}\rceil }$.

def isqrt ( n : int , x0 : int = 1 ) -> int : """  isqrt mediante iteración de Newton-Heron con estimación inicial especificada.  Utiliza detección de oscilación de 2 ciclos. Precondiciones:  n >= 0 # isqrt(0) = 0  x0 > 0, por defecto es 1 # estimación inicial Salida:  isqrt(n)  """ assert n >= 0 and x0 > 0 , "Entrada no válida"# isqrt(0) = 0; isqrt(1) = 1 si n < 2 : devuelve nprev2 = - 1 # x_{i-2} prev1 = x0 # x_{i-1}mientras sea verdadero : x1 = ( prev1 + n // prev1 ) // 2# Caso 1: convergencia (valor estable) si x1 == prev1 : devolver x1# Caso 2: oscilación (ciclo 2) si x1 == prev2 y x1 != prev1 : # Estamos alternando entre prev1 y prev2 # Elegimos el menor (la raíz cuadrada del entero verdadero) devolver min ( prev1 , x1 )# Avanzar prev2 , prev1 = prev1 , x1

La llamada isqrt(2000000)converge a${\displaystyle 1414}$en 14 pasos a través de while:

${\displaystyle {\begin{aligned}&1000000\to 500001\to 250002\to 125004\to 62509\to 31270\to 15666\to 7896\to 4074\to 2282\\&\to 1579\to 1422\to 1414\rightarrow 1414\end{aligned}}}$.

Se obtiene una iteración al establecer x0en${\displaystyle \lfloor n/2\rfloor }$con la llamada isqrt(2000000, 1000000). Aunque el método de Heron converge cuadráticamente cerca de la solución, se gana menos de un bit de precisión por iteración al principio. Esto significa que la elección de la estimación inicial es crítica para el rendimiento del algoritmo. ^{[ nota 5 ]} Cuando se dispone de un cálculo rápido para la parte entera del logaritmo binario o para la longitud en bitsn.bit_length() (como por ejemplo en Python ), es mejor comenzar en$${\displaystyle x_{0}=2^{\lfloor (\log _{2}n)/2\rfloor +1},}$$que es la menor potencia de dos mayor que${\displaystyle {\sqrt {n}}}$. En el ejemplo de la raíz cuadrada entera de 2000000 ,${\displaystyle \lfloor \log _{2}n\rfloor =20}$,${\displaystyle x_{0}=2^{11}=2048}$y la secuencia resultante es $${\displaystyle 2048\rightarrow 1512\rightarrow 1417\rightarrow 1414\rightarrow 1414.}$$En este caso solo se necesitan cuatro pasos de iteración. Esto corresponde a la llamada isqrt(2000000, 2048).

Si se trabaja en base 2 , la elección del dígito se simplifica a uno entre 0 (el "candidato pequeño") y 1 (el "candidato grande"), y las manipulaciones de dígitos se pueden expresar en términos de operaciones de desplazamiento binario. Siendo *la multiplicación, el desplazamiento a la izquierda y el desplazamiento lógico a la derecha, un algoritmo recursivo para hallar la raíz cuadrada entera de cualquier número natural es:<<>>

def isqrt_recursive ( n : int ) -> int : assert n >= 0 , "n debe ser un entero no negativo" if n < 2 : return n# Llamada recursiva: small_cand = isqrt_recursive ( n >> 2 ) << 1 # equivalente a 2 * isqrt(n // 4) large_cand = small_cand + 1 if large_cand * large_cand > n : return small_cand else : return large_cand

Programa equivalente no recursivo: ^{[ 2 ] [ nota 6 ]}

def isqrt_iterative ( x : int ) -> int : """ Guy, Martin (1985). "Raíz cuadrada de enteros rápida mediante el algoritmo del ábaco del Sr. Woo" """ assert x >= 0 , "x debe ser un entero no negativo"op = x ; res = 0# nota: i << 1 es 2 * i, i << 2 es 4 * i, # i >> 1 es i // 2, e i >> 2 es i // 4...# "uno" comienza en la potencia más alta de cuatro <= x uno = 1 mientras uno <= op : uno <<= 2 uno >>= 2mientras uno != 0 : dltasqr = res + uno si op >= dltasqr : op -= dltasqr res += uno << 1 res >>= 1 uno >>= 2retorno res

Véase Métodos para calcular raíces cuadradas §  Sistema de numeración binario (base 2) para un ejemplo. ^{[ 2 ]}

Paul Zimmermann (1999) proporciona el siguiente algoritmo. ^{[ 3 ]}

${\displaystyle {\text{Algorithm }}{\text{SqrtRem}}(n=a_{3}b^{3}+a_{2}b^{2}+a_{1}b+a_{0})}$

${\displaystyle {\text{Input: }}0\leq a_{i}<b,{\text{ with }}a_{3}\geq b/4}$

${\displaystyle {\text{Output: }}(s,r){\text{ such that }}s^{2}\leq n=s^{2}+r<(s+1)^{2}}$

${\displaystyle \qquad (s',r')\gets {\text{SqrtRem}}(a_{3}b+a_{2})}$

${\displaystyle \qquad (q,u)\gets {\text{DivRem}}(r'b+a_{1},2s')}$

${\displaystyle \qquad s\gets s'b+q}$

${\displaystyle \qquad r\gets ub+a_{0}-q^{2}}$

${\displaystyle \qquad {\text{if }}r<0{\text{ then }}}$

${\displaystyle \qquad \qquad r\gets r+2s-1}$

${\displaystyle \qquad \qquad s\gets s-1}$

${\displaystyle \qquad {\text{return }}(s,r)}$

Debido a que solo se realiza una llamada recursiva por nivel, la complejidad total permanece${\displaystyle O(n)}$en número de bits. Cada nivel realiza únicamente aritmética de tiempo lineal sobre números de tamaño reducido.

La raíz cuadrada al estilo Karatsuba se utiliza principalmente para aritmética de precisión arbitraria en enteros muy grandes, donde se combina eficientemente con la división de Burnikel-Ziegler y la multiplicación de Karatsuba . Fue analizada formalmente por primera vez por Paul Zimmermann (1999). ^{[ 3 ]} Trabajos prácticos anteriores incluyen a Martin Guy (1985), ^{[ 2 ]} y versiones recursivas aparecen en Donald Knuth (1998). ^{[ 4 ]} Las bibliotecas modernas GMP y MPIR implementan técnicas recursivas similares.

El siguiente programa en Python implementa el algoritmo de Zimmermann. Dado un número entero${\displaystyle n\geq 0}$, SqrtRemcalcula simultáneamente su raíz cuadrada entera${\displaystyle s=\lfloor {\sqrt {n}}\rfloor }$y el resto correspondiente${\displaystyle r=n-s^{2}}$La elección isqrt()es a voluntad . ^{[ nota 5 ]}

def SqrtRem ( n : int , word_bits : int = 32 ) -> tuple [ int , int ]: """  Implementación basada en el algoritmo de raíz cuadrada entera estilo Karatsuba de Zimmermann  [Zimmermann, 1999]. Divide recursivamente la entrada n en "ramas"  de tamaño `word_bits` y combina resultados parciales para calcular la raíz cuadrada entera. Argumentos:  n (int): Entero no negativo del que calcular la raíz cuadrada.  word_bits (int, opcional): Número de bits por "rama" o fragmento  utilizado al dividir recursivamente n. El valor predeterminado es 32. Cada  rama representa una parte de tamaño fijo de n para el algoritmo. Devuelve:  tupla[int, int]: s = raíz cuadrada entera de n, r = resto (n - s*s). Notas:  El tamaño de la rama controla la granularidad de la recursión. Un valor mayor de word_bits  reduce la profundidad de la recursión, pero aumenta el tamaño de los subproblemas;  un valor menor de word_bits aumenta la profundidad de la recursión, pero trabaja con fragmentos más pequeños. Referencia:  Zimmermann, P. (1999). "Karatsuba Square Root", Informe de investigación n.° 3805, Inria.  Archivado en: https://inria.hal.science/inria-00072854v1/file/RR-3805.pdf  """ if n < 0 : raise ValueError ( "n debe ser no negativo" ) if n == 0 : return 0 , 0 # caso trivial# Determinar el número de ramas del tamaño de una palabra (imita la división de "ramas" en Zimmermann) limblen = ( n . bit_length () + word_bits - 1 ) // word_bits# Caso base: extremidad única — calcular directamente si limblen <= 1 : s = isqrt ( n ) # cualquier isqrt, por ejemplo, math.isqrt o personalizado r = n - s * s return s , r# --- Paso 1: Dividir n en partes alta y baja --- half_limbs = limblen // 2 shift = half_limbs * word_bits hi = n >> shift # mitad alta, corresponde a a3*b + a2 lo = n & (( 1 << shift ) - 1 ) # mitad baja, corresponde a a1*b + a0# --- Paso 2: Llamada recursiva en la parte superior --- s_high , r_high = SqrtRem ( hi , word_bits ) # raíz cuadrada aproximada de la mitad superior# --- Paso 3: Recombinar para aproximar la raíz cuadrada completa --- cuarto = desplazamiento // 2 numerador = ( r_alto << cuarto ) | ( lo >> cuarto ) # simular el paso DivRem de Zimmermann denominador = s_alto << 1 # término de 2*sq = numerador // denominador si denominador sino 0 # división entera s_candidato = ( s_alto << cuarto ) + q # recombinar alto y bajo# --- Paso 4: Verificación y corrección --- # Asegúrese de que el resto sea no negativo y s*s <= n < (s+1)*(s+1) s = s_candidate r = n - s * smientras r < 0 : # corrección de sobreestimación s -= 1 r = n - s * s mientras ( s + 1 ) * ( s + 1 ) <= n : # corrección de subestimación s += 1 r = n - s * sdevolver s , r

Ejemplo de uso

para n en [( 2 ** 32 ) + 5 , 12345678901234567890 , ( 1 << 1512 ) - 1 ]: s , r = SqrtRem ( n ) print ( f "SqrtRem( { n } ) = { s } , resto = { r } " )

En la entrada${\displaystyle c}$, un entero no negativo, el siguiente programa calcula la fracción continua simple de${\displaystyle {\sqrt {c}}}$La raíz cuadrada entera${\displaystyle \lfloor {\sqrt {c}}\rfloor }$se calcula una vez. ^{[ nota 5 ]} Solo se utiliza aritmética de enteros. El programa muestra${\displaystyle [a0,(a1,a2,...,am)]}$, donde el segundo elemento es la parte periódica.

def continue_fraction_sqrt ( c : int ) -> tuple [ int , tuple [ int , ... ]]: """  Calcula la fracción continua de sqrt(c) usando aritmética de enteros.  Devuelve [a0, (a1, a2, ..., am)] donde el segundo elemento es la parte periódica.  Para cuadrados perfectos, el período está vacío.  """ a0 = isqrt ( c ) # Cuadrado perfecto: devuelve el período vacío if a0 * a0 == c : return ( a0 , ())m , d , a = 0 , 1 , a0 periodo = [] visto = conjunto ()Mientras sea verdadero : m_next = d * a - m d_next = ( c - m_next * m_next ) // d a_next = ( a0 + m_next ) // d_nextSi ( m_next , d_next , a_next ) está en seen : breakvisto . agregar (( m_next , d_next , a_next )) período . agregar ( a_next ) m , d , a = m_next , d_next , a_nextdevolver ( a0 , tupla ( período ))

Ejemplo de uso

para c en lista ( rango ( 0 , 18 )) + [ 114 ] + [ 4097280036 ]: cf = continue_fraction_sqrt ( c ) print ( f "sqrt( { c } ): { cf } " )

Producción