Implementación de la función de cálculo de raíz

Question

La implementación de funciones matemáticas para varias cosas es bastante simple. int mul(int,int);, int pow(int,int);, incluso double div(float,float); son fáciles de hacer y se pueden implementar con bucles o recursión. (Estos son los mismos métodos que se utilizan para realizar estas funciones a mano o en la cabeza). Para multiplicar, simplemente agregue el número repetidamente. Para dividir, restarlo repetidamente. Para obtener el poder, multiplícate repetidamente. Y así.

Sin embargo, una función matemática que siempre me he preguntado son las raíces. Por ejemplo, ¿cómo escribirías una función para calcular la raíz cuadrada (o cubo, etc.) de un número (es decir, double root(float num, float root);)? Traté de mirar alrededor y no pude encontrar un algoritmo o método para hacerlo.

Cuando trato de calcular una raíz a mano, generalmente utilizo el método de adivinar (comience con un número aproximado, agregue una fracción, multiplique, vea qué tan lejos está, agregue una fracción más pequeña, multiplique, vuelva a verificar, y repite hasta que esté satisfecho). Supongo que podría funcionar, pero seguramente hay un método mejor y más rápido (independientemente de la velocidad con la que una computadora pueda hacerlo).

Obviamente, las LUT no son relevantes ya que tendrían que ser lo suficientemente genéricas como para tomar cualquier operando (a menos que esté escribiendo un juego con un conjunto finito de datos). El Wikipedia article menciona el método de suposición y enumera algunos antiguos (desde mucho antes de que se inventaran las computadoras), así como algunos métodos matemáticos puros e incluso de cálculo (incluidos algunos que tienen como componente el "infinito"). Los únicos que parecen tener algo que ver con la electrónica usan trucos o logiritmos. (Y eso es solo para square-roots, y mucho menos cube-roots, etc.)

¿No hay un método de cálculo de raíz fácil? ¿Cómo lo hacen las calculadoras? ¿Cómo lo hacen las computadoras? (No, simplemente haciendo double pow(a,0.5); no funcionará porque entonces ¿cómo se implementaría double pow(float,float)?)

¿Acabo de agrupar incorrectamente funciones de raíz con funciones más simples? ¿Son más complejos de lo que parecen?

Answer 1

Hay un N-th root algorithm simple de implementar vinculado directamente desde el artículo que indique. Se deriva de Newton's method.

Answer 2

Hay un par de posibilidades. Hay un par de enfoques iterativos diferentes, como la bisección o Newton. En lo que respecta al uso de pow, algunas computadoras (x86, por ejemplo) tienen instrucciones para hacer (al menos parte de) elevar un número a una potencia, por lo que se trata simplemente de escribir un poco de marco sobre eso.

Aquí hay una implementación en lenguaje ensamblador del método de Newton para una raíz cuadrada, en este caso trabajando solo con enteros de 16 bits, pero la misma idea básica se aplica a otros tipos. Escribí esto hace unos 20 años, por lo que era para CPU de 16 bits sin hardware de coma flotante.

isqrt proc uses di, number:word 
; 
; uses bx, cx, dx 
; 
    mov di,number 
    mov ax,255 
start_loop: 
    mov bx,ax 
    xor dx,dx 
    mov ax,di 
    div bx 
    add ax,bx 
    shr ax,1 
    mov cx,ax 
    sub cx,bx 
    cmp cx,2 
    ja start_loop 
    ret 
isqrt endp 

Aquí hay un código para un x87 para calcular poderes arbitrarios:

pow proc 
    ; x^y = 2^(log2(x) * y) 
    fyl2x  
    fld st(0) 
    frndint 
    fld1  
    fscale 
    fxch st(2) 
    fsubrp 
    f2xm1  
    fld1  
    faddp  
    fmulp  
    ret 
endp 

Tenga en cuenta, sin embargo, que no es común que desea implementar la multiplicación con sólo añadir varias veces, o división con sólo resta repetida . Por el contrario, desea desplazar y sumar/restar potencias sucesivas de dos para obtener el resultado mucho más rápido.

Aquí hay algo de código que muestra la idea general:

mult proc 
; multiplies eax by ebx and places result in edx:ecx 
    xor ecx, ecx 
    xor edx, edx 
mul1: 
    test ebx, 1 
    jz mul2 
    add ecx, eax 
    adc edx, 0 
mul2: 
    shr ebx, 1 
    shl eax, 1 
    test ebx, ebx 
    jnz mul1 
done: 
    ret 
mult endp 

Esto es bastante inútil para x86, ya que tiene instrucciones de multiplicación incorporados, pero en los procesadores de mayor edad (PDP-11, 8080, 6502, etc.) código como este era bastante común.

Answer 3

Depende de qué tan general desee ser. Si desea calcular, por ejemplo, (-4.2) ^0.23, va a necesitar una aritmética compleja. Como señala Mat, existen algoritmos rápidos para calcular x ^1/n para enteros ny positivos x. Si quiere x ^y para x positivo y cualquier y, entonces los registros y exponenciales funcionarán.