python - Generando dígitos de raíz cuadrada de 2

algorithm numerical-methods

Question

Quiero generar los dígitos de la raíz cuadrada de dos a 3 millones de dígitos.

Soy consciente de Newton-Raphson pero no tengo mucha idea de cómo implementarlo en C o C ++ debido a la falta de soporte de Biginteger. ¿Puede alguien apuntarme en la dirección correcta?

Además, si alguien sabe cómo hacerlo en python (soy un principiante), también lo apreciaría.

Answer 1

¿Para el trabajo? ¡Usa una biblioteca!

¿Por diversión? Bien por usted :)

Escribe un programa para imitar lo que harías con lápiz y papel. Comience con 1 dígito, luego 2 dígitos, luego 3, ..., ...

No te preocupes por Newton o por nadie más. Sólo hazlo a tu manera.

Answer 2

Aquí hay una función de raíz cuadrada entera más eficiente (en Python 3.x) que debe terminar en todos los casos. Comienza con un número mucho más cercano a la raíz cuadrada, por lo que toma menos pasos. Tenga en cuenta que int.bit_length requiere Python 3.1+. La comprobación de errores se dejó fuera por brevedad.

def isqrt(n):
    x = (n >> n.bit_length() // 2) + 1
    result = (x + n // x) // 2
    while abs(result - x) > 1:
        x = result
        result = (x + n // x) // 2
    while result * result > n:
        result -= 1
    return result

Answer 3

Aquí hay una versión corta para calcular la raíz cuadrada de un número entero a dígitos de precisión. Funciona encontrando la raíz cuadrada entera de a después de multiplicar por 10 elevados a los 2 dígitos .

def sqroot(a, digits):
    a = a * (10**(2*digits))
    x_prev = 0
    x_next = 1 * (10**digits)
    while x_prev != x_next:
        x_prev = x_next
        x_next = (x_prev + (a // x_prev)) >> 1
    return x_next

Sólo algunas advertencias.

Deberá convertir el resultado en una cadena y agregar el punto decimal en la ubicación correcta (si desea que se imprima el punto decimal).

Convertir un entero muy grande en una cadena no es muy rápido.

La división de enteros muy grandes tampoco es muy rápida (en Python).

Dependiendo del rendimiento de su sistema, puede tomar una hora o más para calcular la raíz cuadrada de 2 a 3 millones de lugares decimales.

No he probado que el bucle siempre termine. Puede oscilar entre dos valores que difieren en el último dígito. O puede que no.

Answer 4

Bueno, el siguiente es el código que escribí. Generó un millón de dígitos después del decimal para la raíz cuadrada de 2 en unos 60800 segundos para mí, pero mi computadora portátil estaba durmiendo cuando estaba ejecutando el programa, debería ser más rápido que eso. Puede intentar generar 3 millones de dígitos, pero puede tardar un par de días en obtenerlo.

def sqrt(number,digits_after_decimal=20):
import time
start=time.time()
original_number=number
number=str(number)
list=[]
for a in range(len(number)):
    if number[a]==''.'':
        decimal_point_locaiton=a
        break
    if a==len(number)-1:
        number+=''.''
        decimal_point_locaiton=a+1
if decimal_point_locaiton/2!=round(decimal_point_locaiton/2):
    number=''0''+number
    decimal_point_locaiton+=1
if len(number)/2!=round(len(number)/2):
    number+=''0''
number=number[:decimal_point_locaiton]+number[decimal_point_locaiton+1:]
decimal_point_ans=int((decimal_point_locaiton-2)/2)+1
for a in range(0,len(number),2):
    if number[a]!=''0'':
        list.append(eval(number[a:a+2]))
    else:
        try:
            list.append(eval(number[a+1]))
        except IndexError:
            pass
p=0
c=list[0]
x=0
ans=''''
for a in range(len(list)):
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    try:
        c=c*100+list[a+1]
    except IndexError:
        c=c*100
while c!=0:
    x=0
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    c=c*100
    if len(ans)-decimal_point_ans>=digits_after_decimal:
            break
ans=ans[:decimal_point_ans]+''.''+ans[decimal_point_ans:]
total=time.time()-start
return ans,total

Answer 5

En cuanto a números grandes arbitrarios, puede consultar la GMP (para C / C ++).

Answer 6

La mejor manera es probablemente usando la expansión de fracción continua [1; 2, 2, ...] [1; 2, 2, ...] la raíz cuadrada de dos.

def root_two_cf_expansion():
    yield 1
    while True:
        yield 2

def z(a,b,c,d, contfrac):
    for x in contfrac:
        while a > 0 and b > 0 and c > 0 and d > 0:
            t = a // c
            t2 = b // d
            if not t == t2:
                break
            yield t
            a = (10 * (a - c*t))
            b = (10 * (b - d*t))
            # continue with same fraction, don''t pull new x
        a, b = x*a+b, a
        c, d = x*c+d, c
    for digit in rdigits(a, c):
        yield digit

def rdigits(p, q):
    while p > 0:
        if p > q:
           d = p // q
           p = p - q * d
        else:
           d = (10 * p) // q
           p = 10 * p - q * d
        yield d

def decimal(contfrac):
    return z(1,0,0,1,contfrac)

decimal((root_two_cf_expansion()) devuelve un iterador de todos los dígitos decimales. t1 y t2 en el algoritmo son valores mínimos y máximos del siguiente dígito. Cuando son iguales, emitimos ese dígito.

Tenga en cuenta que esto no maneja ciertos casos excepcionales, como números negativos en la fracción continua.

(Este código es una adaptación del código de Haskell para manejar fracciones continuas que han estado flotando alrededor).

Answer 7

Podrías intentar usar el mapeo:

a/b -> (a+2b)/(a+b) comenzando con a= 1, b= 1 . Esto converge a sqrt (2) (de hecho, da las representaciones de fracción continua de él).

Ahora el punto clave: esto puede representarse como una multiplicación de matriz (similar a fibonacci)

Si a_n y b_n son los números n en los pasos, entonces

[1 2] [a_n b_n] ^T = [a_ (n + 1) b_ (n + 1)] ^T
[1 1]

que ahora nos da

[1 2] ⁿ [a_1 b_1] ^T = [a_ (n + 1) b_ (n + 1)] ^T
[1 1]

Por lo tanto, si la matriz de 2x2 es A, debemos calcular A ^{n, lo} que se puede hacer mediante la cuadratura repetida y solo usamos aritmética de enteros (para que no tenga que preocuparse por los problemas de precisión).

También tenga en cuenta que el a / b que obtenga siempre estará en forma reducida (como gcd (a, b) = gcd (a + 2b, a + b)), así que si está pensando en usar una clase de fracción para representar el intermedio resultados, no!

Dado que el nth denominadores es como (1 + sqrt (2)) ^ n, para obtener 3 millones de dígitos es probable que tenga que calcular hasta el 3671656 ^th término.

Tenga en cuenta que, aunque esté buscando el término de ~ 3.6 millones, la cuadratura repetida le permitirá calcular el enésimo término en las multiplicaciones y adiciones de O (Log n).

Además, esto puede hacerse fácilmente en paralelo, a diferencia de los iterativos como Newton-Raphson, etc.

Answer 8

Python ya soporta grandes enteros fuera de la caja, y si eso es lo único que lo retiene en C / C ++, siempre puede escribir una clase de contenedor rápido.

El único problema que mencionaste es la falta de grandes enteros. Si no quieres usar una biblioteca para eso, ¿estás buscando ayuda para escribir una clase así?

Answer 9

EDITAR : Me gusta esta versión mejor que la anterior. Es una solución general que acepta tanto números enteros como fracciones decimales; con n = 2 y precisión = 100000, toma alrededor de dos minutos. Gracias a Paul McGuire por sus sugerencias y otras sugerencias de bienvenida!

def sqrt_list(n, precision):
    ndigits = []        # break n into list of digits
    n_int = int(n)
    n_fraction = n - n_int

    while n_int:                            # generate list of digits of integral part
        ndigits.append(n_int % 10)
        n_int /= 10
    if len(ndigits) % 2: ndigits.append(0)  # ndigits will be processed in groups of 2

    decimal_point_index = len(ndigits) / 2  # remember decimal point position
    while n_fraction:                       # insert digits from fractional part
        n_fraction *= 10
        ndigits.insert(0, int(n_fraction))
        n_fraction -= int(n_fraction)
    if len(ndigits) % 2: ndigits.insert(0, 0)  # ndigits will be processed in groups of 2

    rootlist = []
    root = carry = 0                        # the algorithm
    while root == 0 or (len(rootlist) < precision and (ndigits or carry != 0)):
        carry = carry * 100
        if ndigits: carry += ndigits.pop() * 10 + ndigits.pop()
        x = 9
        while (20 * root + x) * x > carry:
                x -= 1
        carry -= (20 * root + x) * x
        root = root * 10 + x
        rootlist.append(x)
    return rootlist, decimal_point_index