Come posso migliorare questo metodo di radice quadrata?

Question

So che sembra un compito a casa, ma non lo è. Ultimamente mi sono interessato agli algoritmi utilizzati per eseguire determinate operazioni matematiche, come seno, radice quadrata, ecc. Al momento, sto cercando di scrivere lo Babylonian method di calcolare le radici quadrate in C#.

Finora, ho questo: metodo ogni volta

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; 
} 

funziona bene e produce la stessa risposta esatta come Math.Sqrt del .NET Framework(). Come puoi probabilmente intuire, però, è più lento del metodo nativo (di circa 800 tick). So che questo particolare metodo non sarà mai più veloce del metodo nativo, ma mi chiedo solo se ci sono delle ottimizzazioni che posso fare.

L'unica ottimizzazione che ho visto immediatamente era il fatto che il calcolo sarebbe stato eseguito 100 volte, anche dopo che la risposta era già stata determinata (a quel punto, r avrebbe sempre lo stesso valore). Così, ho aggiunto un controllo rapido per vedere se il nuovo valore calcolato è lo stesso del valore calcolato in precedenza e uscire dal ciclo. Sfortunatamente, non ha fatto molta differenza in termini di velocità, ma mi è sembrata la cosa giusta da fare.

E prima di dire "Perché non usare solo Math.Sqrt() invece?" ... Lo sto facendo come esercizio di apprendimento e non intendo utilizzare effettivamente questo metodo in alcun codice di produzione.

Answer 1

Prima, invece di controllare l'uguaglianza (r == ultimo), si dovrebbe verificherà la convergenza, in cui r è vicino a durare, dove vicino è definito da un epsilon arbitraria:

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

Come l'articolo di wikipedia che hai collegato a menzioni - non si calcolano in modo efficiente le radici quadrate con il metodo di Newton - invece, si usano i logaritmi.

Answer 2

Invece di interrompere il ciclo e quindi restituire r, è possibile solo restituire r. Potrebbe non fornire alcun aumento notevole delle prestazioni.

Answer 3

Si può provare r = x >> 1;

invece di/2 (anche in altro luogo voi dispositivo per 2). Potrebbe darti un leggero vantaggio. Vorrei anche spostare il 100 nel ciclo. Probabilmente niente, ma stiamo parlando di zecche qui.

basta controllarlo ora.

MODIFICA: Corretto il> in >>, ma non funziona per il doppio, quindi non importa. l'inlining del 100 non mi ha dato alcun aumento di velocità.

Answer 4

Quello che stai facendo qui è eseguire Newton's method of finding a root. Quindi potresti semplicemente usare un algoritmo di individuazione delle radici più efficiente. È possibile iniziare a cercarlo here.

Answer 5

Sostituire la divisione di 2 con uno spostamento di bit è improbabile che faccia una grande differenza; dato che la divisione è costante, spero che il compilatore sia abbastanza intelligente da farlo per te, ma puoi anche provare a vederlo.

È molto più probabile ottenere un miglioramento uscendo dal ciclo in anticipo, quindi è possibile memorizzare nuovo valore r in una variabile e confrontarlo con il vecchio r, o memorizzare x/r in una variabile e confrontarlo con r prima di fare l'aggiunta e la divisione.

Answer 6

Definire una tolleranza e tornare in anticipo quando le iterazioni successive rientrano in tale tolleranza.

Answer 7

Bene, la funzione nativa Sqrt() probabilmente non è implementata in C#, molto probabilmente verrà eseguita in un linguaggio di basso livello e utilizzerà sicuramente un algoritmo più efficiente. Quindi cercare di eguagliare la sua velocità è probabilmente inutile.

Tuttavia, per quanto riguarda il solo tentativo di ottimizzare la funzione per l'heckuvit, la pagina di Wikipedia che hai collegato suggerisce "l'ipotesi iniziale" per essere 2^piano (D/2), dove D rappresenta il numero di cifre binarie in il numero. Potresti dare un tentativo, non vedo molto altro che potrebbe essere ottimizzato in modo significativo nel tuo codice.

Answer 8

Dal momento che lei ha detto il codice qui sotto non era abbastanza veloce, provate questo:

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

dovrebbe essere molto accurata e, si spera veloce.

Ecco il codice per la stima iniziale descritta here. Sembra essere abbastanza buono. Utilizzare questo codice e quindi si dovrebbe anche iterare fino a quando i valori convergono all'interno di un epsilon di differenza.

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 il metodo, ogni iterazione raddoppia il numero di bit corretti.

Utilizzando una tabella per ottenere i primi 4 bit (ad esempio), si avranno 8 bit dopo la 1a iterazione, quindi 16 bit dopo il secondo e tutti i bit necessari dopo la quarta iterazione (poiché un archivio double 52 + 1 bit di mantissa).

Per una ricerca tabella, è possibile estrarre la mantissa in [0.5,1 [ed esponente dall'input (utilizzando una funzione come frexp), quindi normalizzare la mantissa in [64,256 [utilizzando la moltiplicazione con una potenza adeguata di 2.

mantissa *= 2^K 
exponent -= K 

Dopo questo, il vostro numero d'ingresso è ancora mantissa*2^exponent. K deve essere 7 o 8, per ottenere un esponente pari. È possibile ottenere il valore iniziale per le iterazioni da una tabella contenente tutte le radici quadrate della parte integrale della mantissa. Eseguire 4 iterazioni per ottenere la radice quadrata r della mantissa. Il risultato è r*2^(exponent/2), costruito utilizzando una funzione come ldexp.

MODIFICA. Ho messo qualche codice C++ sotto per illustrare questo. La funzione sr1 dell'OP con test migliorato impiega 2,78 per calcolare 2^24 radici quadrate; la mia funzione sr2 richiede 1.42 s, e l'hardware sqrt richiede 0.12 s.

#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

Sono stato a guardare questo anche per scopi di apprendimento. Potresti essere interessato a due modifiche che ho provato.

La prima era di utilizzare un primo ordine taylor serie approssimazione in 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; 
    }; 

La seconda era di provare un terzo ordine (più costoso), iterate

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; 
    }; 

ho creato un helper metodo per cronometrare le funzioni

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; 
    } 
} 

Il metodo originale era più veloce, ma ancora, potrebbe b e interessante :)

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;} 

Questa è la mia radice quadrata veloce preferita. In realtà è l'inverso della radice quadrata, ma puoi invertirlo dopo se vuoi .... Non posso dire se è più veloce se vuoi la radice quadrata e non la radice quadrata inversa, ma è strana come la stessa .
http://www.beyond3d.com/content/articles/8/

Answer 12

Sostituzione "/ 2" con "* 0.5" che rende questo ~ 1,5 volte più veloce sulla mia macchina, ma ovviamente non così veloce come l'implementazione nativa.