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
Vai articolo originale: http://carlonicolini.altervista.org/index.php/Informatica-e-Web/Notizie-dal-web/Fast-sqrt-in-C++-algoritmi-per-radici-quadrate-veloci.html