¿Cómo puedo mejorar este método de raíz cuadrada?

Question

Sé que esto suena como una tarea de tarea, pero no lo es. Últimamente me interesaron los algoritmos utilizados para realizar ciertas operaciones matemáticas, como sinusoide, raíz cuadrada, etc. En este momento, estoy tratando de escribir el Babylonian method de cálculo de raíces cuadradas en C#.

Hasta ahora, tengo unas pocas cosas: método cada vez que

public static double SquareRoot(double x) { 
    if (x == 0) return 0; 

    double r = x/2; // this is inefficient, but I can't find a better way 
         // to get a close estimate for the starting value of r 
    double last = 0; 
    int maxIters = 100; 

    for (int i = 0; i < maxIters; i++) { 
     r = (r + x/r)/2; 
     if (r == last) 
      break; 
     last = r; 
    } 

    return r; 
} 

Funciona muy bien y produce exactamente la misma respuesta que Math.Sqrt de .NET Framework(). Sin embargo, como probablemente puedas adivinar, es más lento que el método nativo (alrededor de 800 tics). Sé que este método particular nunca será más rápido que el método nativo, pero me pregunto si hay optimizaciones que pueda hacer.

La única optimización que vi de inmediato fue el hecho de que el cálculo se ejecutaría 100 veces, incluso después de que la respuesta ya hubiera sido determinada (en ese punto, r siempre sería el mismo valor). Entonces, agregué una comprobación rápida para ver si el valor recién calculado es el mismo que el valor calculado previamente y salgo del ciclo. Desafortunadamente, no hizo mucha diferencia en la velocidad, pero parecía lo correcto.

Y antes dices "¿Por qué no utilizas Math.Sqrt() en su lugar?" ... Lo hago como ejercicio de aprendizaje y no pretendo utilizar este método en ningún código de producción.

Answer 1

En primer lugar, en vez de comprobar la igualdad (r == última), se debe comprobar la convergencia, en el que r es cercano a durar, donde cerca se define por una épsilon arbitraria:

eps = 1e-10 // pick any small number 
if (Math.Abs(r-last) < eps) break; 

A medida que la el artículo de la wikipedia que vinculó a las menciones (no calcula de manera eficiente las raíces cuadradas con el método de Newton), en cambio, usa logaritmos.

Answer 2

En lugar de romper el bucle y devolver r, simplemente puede devolver r. No puede proporcionar ningún aumento notable en el rendimiento.

Answer 3

Usted puede intentar r = x >> 1;

en lugar de/2 (también en otro lugar en el que el dispositivo por 2). Podría darle una ligera ventaja. También movería el 100 al lazo. Probablemente nada, pero estamos hablando de tics aquí.

acaba de comprobarlo ahora.

EDITAR: Se corrigió el> en >>, pero no funciona para los dobles, así que nunca importa. la entrada de 100 no me dio ningún aumento de velocidad.

Answer 4

Lo que está haciendo aquí es ejecutar Newton's method of finding a root. Entonces podrías usar un algoritmo de búsqueda de raíz más eficiente. Puede comenzar a buscarlo here.

Answer 5

Reemplazar la división por 2 con un cambio de bit es poco probable que haga una gran diferencia; dado que la división es constante, espero que el compilador sea lo suficientemente inteligente como para hacer eso por usted, pero también puede intentarlo.

Es mucho más probable que obtenga una mejora al salir temprano del ciclo, así que almacene el nuevo r en una variable y compárelo con r anterior, o almacene x/r en una variable y compárelo con r antes de hacer la adición y división.

Answer 6

Defina una tolerancia y regrese pronto cuando las iteraciones subsiguientes caigan dentro de esa tolerancia.

Answer 7

Bueno, la función nativa Sqrt() probablemente no está implementada en C#, lo más probable es que se haga en un lenguaje de bajo nivel, y ciertamente usará un algoritmo más eficiente. Así que tratar de igualar su velocidad es probablemente inútil.

Sin embargo, en lo que respecta a tratar de optimizar su función para el heckbook, la página de Wikipedia que enlazó recomienda la "conjetura inicial" de 2^piso (D/2), donde D representa el número de dígitos binarios en el número. Podría intentarlo, no veo mucho más que pueda ser optimizado significativamente en su código.

Answer 8

ya que dijo el código de abajo no era lo suficientemente rápido, intente esto:

static double guess(double n) 
    { 
     return Math.Pow(10, Math.Log10(n)/2); 
    } 

debe ser muy precisa y es de esperar rápido.

Aquí está el código para la estimación inicial descrita here. Parece ser bastante bueno. Use este código, y luego también debe iterar hasta que los valores converjan dentro de un épsilon de diferencia.

public static double digits(double x) 
    { 
     double n = Math.Floor(x); 
     double d; 

     if (d >= 1.0) 
     { 
      for (d = 1; n >= 1.0; ++d) 
      { 
       n = n/10; 
      } 
     } 
     else 
     { 
      for (d = 1; n < 1.0; ++d) 
      { 
       n = n * 10; 
      } 
     } 


     return d; 
    } 

    public static double guess(double x) 
    { 
     double output; 
     double d = Program.digits(x); 

     if (d % 2 == 0) 
     { 
      output = 6*Math.Pow(10, (d - 2)/2); 
     } 
     else 
     { 
      output = 2*Math.Pow(10, (d - 1)/2); 
     } 

     return output; 
    } 

Answer 9

Con su método, cada iteración dobla el número de bits correctos.

Usando una tabla para obtener los 4 bits iniciales (por ejemplo), tendrá 8 bits después de la 1ra iteración, luego 16 bits después del segundo y todos los bits que necesita después de la cuarta iteración (desde double tiendas 52 + 1 pedacitos de mantisa).

Para una búsqueda de tabla, puede extraer la mantisa en [0.5,1 [y exponente de la entrada (usando una función como frexp), luego normalizar la mantisa en [64,256 [usando la multiplicación por una potencia adecuada de 2.

mantissa *= 2^K 
exponent -= K 

Después de esto, su número de entrada es todavía mantissa*2^exponent. K debe ser 7 u 8 para obtener un exponente par. Puede obtener el valor inicial para las iteraciones a partir de una tabla que contiene todas las raíces cuadradas de la parte integral de la mantisa. Realice 4 iteraciones para obtener la raíz cuadrada r de la mantisa. El resultado es r*2^(exponent/2), construido con una función como ldexp.

EDITAR. Puse un código C++ a continuación para ilustrar esto. La función sr1 del OP con prueba mejorada requiere 2.78s para calcular 2^24 raíces cuadradas; mi función sr2 toma 1.42s, y el hardware sqrt toma 0.12s.

#include <math.h> 
#include <stdio.h> 

double sr1(double x) 
{ 
    double last = 0; 
    double r = x * 0.5; 
    int maxIters = 100; 
    for (int i = 0; i < maxIters; i++) { 
    r = (r + x/r)/2; 
    if (fabs(r - last) < 1.0e-10) 
     break; 
    last = r; 
    } 
    return r; 
} 

double sr2(double x) 
{ 
    // Square roots of values in 0..256 (rounded to nearest integer) 
    static const int ROOTS256[] = { 
    0,1,1,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6, 
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9, 
    9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,11,11,11,11,11, 
    11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,12,12,12, 
    12,12,12,12,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13, 
    13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14, 
    14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15, 
    15,15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16 }; 

    // Normalize input 
    int exponent; 
    double mantissa = frexp(x,&exponent); // MANTISSA in [0.5,1[ unless X is 0 
    if (mantissa == 0) return 0; // X is 0 
    if (exponent & 1) { mantissa *= 128; exponent -= 7; } // odd exponent 
    else { mantissa *= 256; exponent -= 8; } // even exponent 
    // Here MANTISSA is in [64,256[ 

    // Initial value on 4 bits 
    double root = ROOTS256[(int)floor(mantissa)]; 

    // Iterate 
    for (int it=0;it<4;it++) 
    { 
     root = 0.5 * (root + mantissa/root); 
    } 

    // Restore exponent in result 
    return ldexp(root,exponent>>1); 
} 

int main() 
{ 
    // Used to generate the table 
    // for (int i=0;i<=256;i++) printf(",%.0f",sqrt(i)); 

    double s = 0; 
    int mx = 1<<24; 
    // for (int i=0;i<mx;i++) s += sqrt(i); // 0.120s 
    // for (int i=0;i<mx;i++) s += sr1(i); // 2.780s 
    for (int i=0;i<mx;i++) s += sr2(i); // 1.420s 
} 

Answer 10

He estado mirando esto también para propósitos de aprendizaje. Puede que le interesen dos modificaciones que probé.

La primera consistía en utilizar una primera orden de aproximación de Taylor series en x0:

Func<double, double> fNewton = (b) => 
    { 
     // Use first order taylor expansion for initial guess 
     // http://www27.wolframalpha.com/input/?i=series+expansion+x^.5 
     double x0 = 1 + (b - 1)/2; 
     double xn = x0; 
     do 
     { 
      x0 = xn; 
      xn = (x0 + b/x0)/2; 
     } while (Math.Abs(xn - x0) > Double.Epsilon); 
     return xn; 
    }; 

El segundo era intentar un tercer orden (más caro), iterar

Func<double, double> fNewtonThird = (b) => 
    { 
     double x0 = b/2; 
     double xn = x0; 
     do 
     { 
      x0 = xn; 
      xn = (x0*(x0*x0+3*b))/(3*x0*x0+b); 
     } while (Math.Abs(xn - x0) > Double.Epsilon); 
     return xn; 
    }; 

creé un ayudante método para sincronizar las funciones

public static class Helper 
{ 
    public static long Time(
     this Func<double, double> f, 
     double testValue) 
    { 
     int imax = 120000; 
     double avg = 0.0; 
     Stopwatch st = new Stopwatch(); 
     for (int i = 0; i < imax; i++) 
     { 
      // note the timing is strictly on the function 
      st.Start(); 
      var t = f(testValue); 
      st.Stop(); 
      avg = (avg * i + t)/(i + 1); 
     } 
     Console.WriteLine("Average Val: {0}",avg); 
     return st.ElapsedTicks/imax; 
    } 
} 

El método original era más rápido, pero de nuevo, podría b e interesante :)

Answer 11

float InvSqrt (float x){ 
    float xhalf = 0.5f*x; 
    int i = *(int*)&x; 
    i = 0x5f3759df - (i>>1); 
    x = *(float*)&i; 
    x = x*(1.5f - xhalf*x*x); 
    return x;} 

Esta es mi raíz cuadrada rápida favorita. En realidad es el inverso de la raíz cuadrada, pero puedes invertirlo si quieres ... No puedo decir si es más rápido si quieres la raíz cuadrada y no la raíz cuadrada inversa, pero es igual de genial. .
http://www.beyond3d.com/content/articles/8/

Answer 12

Sustitución "/ 2" por "x 0,5" hace que este ~ 1,5 veces más rápido en mi máquina, pero por supuesto no es tan rápido como la implementación nativa.