Fast sqrt in C++ – algoritmi per radici quadrate veloci

Ho fatto giusto per curiosità alcune prove con complicati algoritmi per il calcolo (a dire degli autori) veloce delle radici quadrate sia di float che di double

Il programma è il seguente:

Penso di avere sfatato un mito

#include <cmath>
#include <iostream>
#include <cstdio>
#include <ctime>


#define SQRT_MAGIC_F 0x5f3759df
// 64  Bit float magic number
#define SQRT_MAGIC_D 0x5fe6ec85e7de30da

using namespace std;
static float fsqrt (float y)
{
    float x, z, tempf;
    size_t i=0;
    unsigned long *tfptr = ((unsigned long *)&tempf) + 1;

    tempf = y;
    *tfptr = (0xbfcdd90a – *tfptr) >> 1;
    x =  tempf;
    z = y*0.5;
    for (i=0; i<5; i++)
        x = x*1.5 – x*x*x*z;
   
    return x*y;
}

inline float fastSqrt_2(const float x) 
{
  union
  {
    int i;
    float x;
  } u;
  u.x = x;
  u.i = (1<<29) + (u.i >> 1) – (1<<22);
  return u.x;
}

inline double fastSqrt_2d(const double x) 
{
  union
  {
    long i;
    double x;
  } u;
  u.x = x;
  u.i = (((long)1)<<61) + (u.i >> 1) – (((long)1)<<51);
  return u.x;
}

inline float invSqrt(const float x)
{
  const float xhalf = 0.5f*x;
 
  union // get bits for floating value
  {
    float x;
    int i;
  } u;
  u.x = x;
  u.i = SQRT_MAGIC_F – (u.i >> 1);  // gives initial guess y0
  return u.x*(1.5f – xhalf*u.x*u.x);// Newton step, repeating increases accuracy
}

inline float fastSqrt_Q3(const float x)
{
  return x * invSqrt(x);
}

int main(void)
{
   
    clock_t start, end;
    float xf=14.42;
    float y1f,y2f;
    double y1d,y2d;   
    double xd=14.42;
   
    cerr << “– FLOAT VERSIONS –“<< endl;
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1f=fastSqrt_2(xf);
    }
    }   
    end=clock();
    printf(“SQRTF: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1f=fastSqrt_Q3(xf);
    }
    }   
    end=clock();
    printf(“SQRTF Q3: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1f=sqrt(xf);
    }
    }   
    end=clock();
    printf(“SQRT: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
   
    cerr << “– DOUBLE VERSIONS –“<< endl;
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1d=fastSqrt_2d(xd);
    }
    }   
    end=clock();
    printf(“SQRTF: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1d=(float)fastSqrt_Q3((float)xd);
    }
    }   
    end=clock();
    printf(“SQRTF Q3: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
    start=clock();
    for (int j=0; j<1E2; j++)
    {
    for (int i=0; i<1E5; i++)
    {
        y1d=sqrt(xd);
    }
    }   
    end=clock();
    printf(“SQRT: %.4lf seconds of processing\n”, (end-start)/(double)CLOCKS_PER_SEC);
   
       
    return 0;
}

 

Senza ottimizzazione in fase di compilazione c’è una certa differenza mentre compilato con g++ -O3 test.cpp -o test non mostra NESSUNA DIFFERENZA fra la sqrt(x) dell’header di cmath e le loro varie versioni più veloci. Suppongo che il livello di ottimizzazione -O3 sistemi già tutto.

(P.S g++ versione 4.1 compilato su Mac OsX)

Se pensate che ci sia qualcosa di metodologicamente sbagliato non esistate a contattarmi, sarò ben contento di scoprire che esiste qualcosa di meglio della semplice sqrt di #include <cmath>!

Riferimenti:
http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots

http://ilab.usc.edu/wiki/index.php/Fast_Square_Root

Link all'articolo originale: http://carlonicolini.altervista.org/index.php/Informatica-e-Web/Notizie-dal-web/Fast-sqrt-in-C++-algoritmi-per-radici-quadrate-veloci.html

Lascia un commento