Now, if you want a super fast sqrt on x86, you can use Greg Walsh's (in)famous approximation
float Sqrt(float x){
int i = *(int*)&x;
i = 0x5f3759df - (i>>1);
float r = *(float*)&i;
r = r*(1.5f - 0.5f*x*r*r);
return r * x;
}
Just don't ask me how it works...
We have no evidence it's really Greg Walsh the author if i'm not wrong.
I can pretty much guarantee that sqrt is not the source of your performance problems...
Micro-optimizations like "fast SQRT" rarely get you anything these days. Modern processors are so fast that you're going to get much more bang for your buck by looking at memory usage patterns and lowering cache misses.
I was writing some SIMD code the other day and replacing a single instruction with an approximation gave me a 25% performance boost. "Fast SQRT" could be really damn useful if you are calculating 4 billion of them every second.
Just don't ask me how it works...float Sqrt(float x){ int i = *(int*)&x; i = 0x5f3759df - (i>>1); float r = *(float*)&i; r = r*(1.5f - 0.5f*x*r*r); return r * x; }
That code computes a fast approximation to 1/sqrt(x) and multiplies it by x.
In modern CPUs casting from float to int and then back will cause moving data among different registers and potentially memory, which can cause stalls.
I wouldn't be surprised if it performs poorly today.
SSE2 and AXV2 give you integer SIMD instructions that operate on the same register file as the floating point instructions. This means you get float -> int punning for free. See clb's post below, apparently on Intel anyway there is a penalty for switching between float and integer SIMD instructions. It seems this has to do with the way their pipeline is structured and how forwarding of results from one operation to the next is handled.
if (std::sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)) < radius) ... // BAD!
if ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) < radius * radius) ... // GOOD!If you really need a square root fast and you aren't to concerned about the precision you can use this. May actually not be that much faster depending on your architecture.
//use 0x1FBD1DF5 for (0.0, 1.0) RME: 1.42%
//use 0x1FBD22DF for (0.0, 1000.0) RME: 1.62%
inline float sqrt(float f)
{
int i = 0x1FBD1DF5 + (*(int*)&f >> 1);
return *(float*)&i;
}
Edit: If you want to do the same for doubles you might have luck with this magic number 0x1FF776E75B46DF3D.
I can pretty much guarantee that sqrt is not the source of your performance problems...
Micro-optimizations like "fast SQRT" rarely get you anything these days. Modern processors are so fast that you're going to get much more bang for your buck by looking at memory usage patterns and lowering cache misses.
I was writing some SIMD code the other day and replacing a single instruction with an approximation gave me a 25% performance boost. "Fast SQRT" could be really damn useful if you are calculating 4 billion of them every second.
SSE2 and AXV2 give you integer SIMD instructions that operate on the same register file as the floating point instructions. This means you get float -> int punning for free.
There is a crossover penalty, or a "bypass delay" at least on Intel hardware for crossing computation of data between float and int modes, see https://software.intel.com/en-us/forums/topic/279382 , so float -> int punning comes at the cost of one cycle.
If you really need a square root fast and you aren't to concerned about the precision you can use this. May actually not be that much faster depending on your architecture.
//use 0x1FBD1DF5 for (0.0, 1.0) RME: 1.42% //use 0x1FBD22DF for (0.0, 1000.0) RME: 1.62% inline float sqrt(float f) { int i = 0x1FBD1DF5 + (*(int*)&f >> 1); return *(float*)&i; }
Edit: If you want to do the same for doubles you might have luck with this magic number 0x1FF776E75B46DF3D.
Testing this on a Haswell Core i7 5960X, at https://github.com/juj/MathGeoLib/commit/e6f667e89a6b51f7c1e9bf2c35c074f3ad3509ac , with VS 2013 Community update 4 targeting AVX, yields (reposting the above numbers since different computer and compiler):
Carmack's: 6.16 clocks, precision: 1.75e-3
mulss(rsqrtss(x)): 2.08 clocks, precision: 3e-4
Looks like using the SSE1 instructions is better at least on the Haswell.
Carmack's: 6.16 clocks, precision: 1.75e-3
sqrt_ss(x): 5.92 clocks, precision: 8.4e-8
mulss(rsqrtss(x)): 2.08 clocks, precision: 3e-4
0x1FBD1DF5: 2.32 clocks, precision: 4.3e-2
Though I think we all agree that the sqrt itself is rather irrelevant when the number of actual CPU cycles is this close to 1, it is certainly interesting, and I did some tests as well.
For one I think you're measuring loads and stores.. as I managed to reproduce almost exactly those numbers (also Haswell) and the ~2 clocks for rsqrt were identical to mem0 = mem1 (unless your test compensates for the load/store time?)
I glanced at the code and it looked as if it was included in the timing.
I guess in order to make anything like this meaningful at all we have to assume that we're at a point in a program where we have our float values in an XMM register and want to measure only the time it takes to get the sqrt of that value into another XMM register, which means we measure latency.
However if we have the sqrt in a loop and the compiler is smart, it might order operations so that the throughput matters. The loop in your test does independent sqrts, so the second sqrt doesn't wait for the first, and therefore measures throughput, but hides latency.
I tried it in ASM and for thoughput got something like sqrtss ~= 7, mulss rsqrtss ~= 2.5, psrld paddd ~= 1.5.
When measuring latency however, sqrtss ~= 10, mulss rsqrtss ~= 10, psrld paddd ~= 2. Looking at Intels specifications this matches pretty closely, as mulss for example has a throughput of 1 cycle, but a latency of 4 cycles. sqrtss is actually listed as 13 cycles latency (if I managed to look at the right column), which probably means my test wasn't well enough constructed to hide such a long latency.
psrld have 1 latency and 1 throughput, whereas paddd has 1 latency but 0.5 throughput, so it can actually do 2 adds per cycle in throughput.
If you have a long sequence of floating point calculations and at one point want to at minimal latency turn once of your floats into the sqrt of itself and don't care for very much precision, psrld paddd is probably the way to go.
Also note that the doc linked by Chris_F actually talked about GPUs and doing sqrts in a pixel shader for example, and we're discussing using the same methods on completely different architectures. I would guess these things matter a lot more on GPUs nowadays, as we're much more often at a time where a few instructions can actually make a difference in performance, and latency is already hidden as good as possible by the GPU so once it starts to matter parallelism has already been maximized and is no longer improvable.
On CPUs I would guess this case is more rare and much harder to predict in general code... and depending on the code around the sqrt different methods can give best results.
I think a script that brute-force tested every sqrt in a program individually with all these methods would find that performance would benefit from one of them in some cases and from another in another case.
Very good points. My benchmark does include loads and stores (identical in each tested variant) to hot cache, and it stresses more throughput, rather than latency. It is difficult to say which one is more important, since it depends if the computation of the overall algorithm has something else to do to hide the latency. Sometimes it's possible to hide latency that it doesn't matter, and other times it is not. Also indeed, none of the timing applies to GPUs or ARM or anything like that.
