So I just discovered the very interesting Quake III inverse square root hack. After learning how it works and all, I decided to test it. I found that the hack only outperformed math.h 1/sqrt(X) when compiled with optimizations enabled.
The hack's implementation:
float q_sqrt(float x) {
float x2 = x * 0.5F;
int i = *( int* )&x; // evil floating point bit hack
i = 0x5f3759df - (i >> 1); // what the fuck?
x = *( float* )&i;
x = x * ( 1.5F - ( (x2 * x * x) ) ); //1st iteration
//y = y * ( 1.5F - ( (x2 * y * y) ) ); //2nd iteration, can be removed
return x;
}
To test how fast 1/sqrt(x) runs compared to q_sqrt(x):
//qtest.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
/*
Implementation of 1/sqrt(x) used in tue quake III game
*/
float q_sqrt(float x) {
float x2 = x * 0.5F;
int i = *( int* )&x; // evil floating point bit hack
i = 0x5f3759df - (i >> 1); // what the fuck?
x = *( float* )&i;
x = x * ( 1.5F - ( (x2 * x * x) ) ); //1st iteration
//y = y * ( 1.5F - ( (x2 * y * y) ) ); //2nd iteration, can be removed
return x;
}
int main(int argc, char *argv[]) {
struct timespec start, stop;
//Will work on floats in the range [0,100]
float maxn = 100;
//Work on 10000 random floats or as many as user provides
size_t num = 10000;
//Bogus
float ans = 0;
//Measure nanoseconds
size_t ns = 0;
if (argc > 1)
num = atoll(argv[1]);
if (num <= 0) return -1;
//Compute "num" random floats
float *vecs = malloc(num * sizeof(float));
if (!vecs) return -1;
for (int i = 0; i < num; i++)
vecs[i] = maxn * ( (float)rand() / (float)RAND_MAX );
fprintf(stderr, "Measuring 1/sqrt(x)\n");
clock_gettime( CLOCK_REALTIME, &start);
for (size_t i = 0; i < num; i++)
ans += 1 / sqrt(vecs[i]);
clock_gettime( CLOCK_REALTIME, &stop);
ns = ( stop.tv_sec - start.tv_sec ) * 1E9 + ( stop.tv_nsec - start.tv_nsec );
fprintf(stderr, "1/sqrt(x) took %.6f nanosecods\n", (double)ns/num );
fprintf(stderr, "Measuring q_sqrt(x)\n");
clock_gettime( CLOCK_REALTIME, &start);
for (size_t i = 0; i < num; i++)
ans += q_sqrt(vecs[i]);
clock_gettime( CLOCK_REALTIME, &stop);
ns = ( stop.tv_sec - start.tv_sec ) * 1E9 + ( stop.tv_nsec - start.tv_nsec );
fprintf(stderr, "q_sqrt(x) took %.6f nanosecods\n", (double)ns/num );
//Side by side
//for (size_t i = 0; i < num; i++)
// fprintf(stdout, "%.6f\t%.6f\n", 1/sqrt(vecs[i]),q_sqrt(vecs[i]));
free(vecs);
}
On my system (Ryzen 3700X) I get:
gcc -Wall -pedantic -o qtest qtest.c -lm
./qtest
Measuring 1/sqrt(x)
1/sqrt(x) took 4.470000 nanosecods
Measuring q_sqrt(x)
q_sqrt(x) took 4.859000 nanosecods
gcc -Wall -pedantic -O1 -o qtest qtest.c -lm
./qtest
Measuring 1/sqrt(x)
1/sqrt(x) took 0.378000 nanosecods
Measuring q_sqrt(x)
q_sqrt(x) took 0.497000 nanosecods
gcc -Wall -pedantic -O2 -o qtest qtest.c -lm
qtest.c: In function ‘q_sqrt’:
qtest.c:11:14: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
11 | int i = *( int* )&x; // evil floating point bit hack
|
qtest.c:13:10: warning: dereferencing type-punned pointer will break strict-aliasing rules [-Wstrict-aliasing]
13 | x = *( float* )&i;
|
./qtest
Measuring 1/sqrt(x)
1/sqrt(x) took 0.500000 nanosecods
Measuring q_sqrt(x)
q_sqrt(x) took 0.002000 nanosecods
My expectation was that q_sqrt(x) was going to work better than 1/sqrt(X) out of the box. After reading some more I now know that either libm is way better optimized or my CPU is equipped with a hardware solution for sqrt(X). After all, CPUs have changes by leaps and bounds since the development of the quick inverse root hack.
What I don't understand is what type of optimizations would the compiler be applying to make it so much faster. Of course maybe my benchmark is ill conceived?
Thanks for any help!!
As you said, most of modern CPUs include a Floating Point Unit that usually provides a hardware instruction to compute square root. FPUs also provide division instructions so I would expect your processor (although I don't know it) to be able to compute an inverse sqrt in only a few assembly instructions. Your results are a bit surprising: you should check whether the FPU is really used. I don't know Ryzen but on ARM processors you can compile your software to use either hardware floating point instructions or software libraries.
Now to answer your questions: GCC optimizations are a complex story and it is usually impossible to predict precisely the effect of a given level on performance. So run some tests as you did, or have a look here for theory.