Fast Square Roots

Native Hardware

The first category gives the right answer using hardware.

We can use the SSE hardware that is present basically everywhere, which boils down to the sqrtss instruction. Running std::sqrt(⋯) will get you this, but with a major performance overhead unless you pass the optimization flags specified above^[6]. You can also call intrinsic functions to access the instruction directly. Although this has disadvantages, it will work without the flags (advantageous especially for MSVC).

/*
Note clarity, auto-vectorization, and bitwise-perfect answer for all inputs.
Change `inline` -> `constexpr` in C++26.

To avoid overhead, should specify:
	"-fno-math-errno -fno-trapping-math" (Clang, GCC)
	"/fp:fast" (MSVC)
*/
#include <cmath>
[[nodiscard]] inline float sqrtf_std( float val ) noexcept
{
	return std::sqrt(val);
}

//x86(-64) only, doesn't require special flags
#include <immintrin.h>
[[nodiscard]] inline float sqrtf_sse( float val ) noexcept
{
	return _mm_cvtss_f32( _mm_sqrt_ss(_mm_set_ps1(val)) );
}

There is also the x87 FPU, which does the operation in 80-bit precision. This is overkill for 32-bit, and the use of the x87 FPU is discouraged anyway. Since compilers (finally, quite sensibly) default to using the SSE unit instead, we must access the x87 through inline assembly or extra compiler flags^[7].

//x87 only (x86(-64) discouraged hardware)
[[nodiscard]] inline float sqrtf_x87( float val ) noexcept
{
	float result;
	#if   defined __clang__ //Clang
		asm(
			"fsqrt\n"
			: "=t"(result)
			: "t"(val)
			: "cc"
		);
	#elif defined __GNUC__ //GCC, etc.
		asm(
			"fld dword ptr %1\n"
			"fsqrt\n"
			: "=&t"(result)
			: "m"(val)
			: "cc"
		);
	#elif defined _MSC_VER && defined _M_IX86 //32-bit MSVC
		__asm
		{
			fld dword ptr val
			fsqrt
			fstp dword ptr result
		}
	#else
		#error
	#endif
	return result;
}

Approximation

The second category gives an approximate answer using a result that we compute ourselves, possibly aided by hardware.

The hardware provides an inverse-square-root function rsqrtss, which by itself is faster. We can therefore compute by \(\sqrt{s}~\)\(= (s)(1/\!\sqrt{s})\):

//Note incorrect result for 0 (and subnorms)
#include <immintrin.h>
[[nodiscard]] inline float rsqrtf_sse( float val ) noexcept
{
	return _mm_cvtss_f32( _mm_rsqrt_ss(_mm_set_ps1(val)) );
}
[[nodiscard]] inline float sqrtf_ssersqrt( float val ) noexcept
{
	return val * rsqrtf_sse(val);
}

There is also an ingenious bithack way to get a square-root. This comes from Wikipedia, perhaps as a simplification of a method first derived by Kahan and Ng in the same unpublished work that led to inverse-square-root.

The explanation (in both locations) is lacking, but it's not notable enough to deserve corrections and the article it needs for a proper treatment, in either academia or on Wikipedia. Roughly speaking, the algorithm works by computing \(\sqrt{s}~\)\(= \exp_2\!\left(\log_2(s)/2\right)\), where we can compute the \(\log_2\) (and so also \(\exp_2\)) efficiently with \(\log_2(s)~\)\(\approx s_\text{int}/2^{23} - 127\). However, I have written an article with the complete explanation, which has been moved to here.

#include <bit>
#include <cstdint>
[[nodiscard]] constexpr float sqrtf_bithack( float val ) noexcept
{
	uint32_t tmp = std::bit_cast<uint32_t>(val);

	constexpr int TWEAK = 0;
	tmp = ((1<<29)-(1<<22)+TWEAK) + (tmp>>1);

	return std::bit_cast<float>(tmp);
}

Note the std::bit_cast<⋯>(⋯) reinterprets the bits to and from floating-point and integers. Commonplace pointer-reinterpreting is dangerously broken in both C and C++, and union hacks are problematic in C++. Before bit-cast, the recommended way was memcpy(⋯). This works surprisingly well, but it does break constexpr-ness.

Newton–Raphson Improvements

Since for this second category the answer is approximate, there is also the option to improve the result using Newton–Raphson iteration (NR, N–R), which for square-root is:

—where \(x_k\) is the current guess and \(s=x^2\) is the number we're trying to take the square root of.

The basic case for NR is of course to do zero iterations:

[[nodiscard]] constexpr float nr0( float /*val*/, float xk ) noexcept { return xk; }

For one iteration, we can just do the formula blithely. The coefficient is nominally 0.5, but small perturbations can result in more accuracy, so we pull it out for later tweaking (see below!):

#include <bit>
[[nodiscard]] constexpr float nr1( float val, float xk ) noexcept
{
	constexpr float coeff = std::bit_cast<float>( 1'056'964'608 ); //`0.5f`
	return coeff*( xk + val/xk );
}

Floating-point division is infamously slow, so we can do the same thing but using the SSE rcpss instruction:

#include <bit>
#include <immintrin.h>
[[nodiscard]] inline float rcp_sse( float val ) noexcept
{
	return _mm_cvtss_f32( _mm_rcp_ss( _mm_set_ps1(val) ) );
}
[[nodiscard]] inline float nr1_ssercp( float val, float xk ) noexcept
{
	constexpr float coeff = std::bit_cast<float>( 1'056'964'608 ); //`0.5f`
	return coeff*( xk + val*rcp_sse(xk) );
}

This is a good use-case for a fused multiply-add (FMA) instruction, which is supported on some AVX and all AVX2 processors. To access it, pass "-mfma"^[8] (Clang, GCC, etc.) or "/arch:AVX2" (MSVC). You could also use the _mm_fmadd_ss(⋯) intrinsic, which might be useful when you have FMA but not AVX2.

//Pass "-mfma" (Clang, GCC) or "/arch:AVX2" (MSVC)
#include <cmath>
#include <bit>
[[nodiscard]] inline float rcp_sse( float val ) noexcept;
[[nodiscard]] inline float nr1_ssercp_fmaA( float val, float xk ) noexcept
{
	constexpr float coeff = std::bit_cast<float>( 1'056'964'608 ); //`0.5f`
	return coeff*std::fma( val,rcp_sse(xk), xk );
}

The final multiplication by the coefficient has to wait for the FMA to finish. To get better pipelining, then, we can try distributing the coefficient. The idea is that even though there are now two multiplies, they can be executing at the same time as the reciprocal:

//Pass "-mfma" (Clang, GCC) or "/arch:AVX2" (MSVC)
#include <cmath>
#include <bit>
[[nodiscard]] inline float rcp_sse( float val ) noexcept;
[[nodiscard]] inline float nr1_ssercp_fmaB( float val, float xk ) noexcept
{
	constexpr float coeff = std::bit_cast<float>( 1'056'964'608 ); //`0.5f`
	return std::fma( coeff*val,rcp_sse(xk), coeff*xk );
}

Finally, here's a version that does two NR iterations. The two 0.5 factors can be combined into a single factor of 0.25:

[[nodiscard]] constexpr float nr2( float val, float xk ) noexcept
{
	constexpr float coeff = std::bit_cast<float>( 1'048'576'000 ); //`0.25f`
	//xk = 0.5f*( xk + val/xk );
	//xk = 0.5f*( xk + val/xk );
	xk =       xk + val/xk;
	xk = coeff*xk + val/xk;
	return xk;
}

Summary Table

Here is a summary table comparing the variants:

I am providing the raw data for the general statistics and latency: statistics.csv, cycles_100iters.csv. These include more statistics than are presented above.

Pipelined Cycles

The performance cost is presented as "Pipelined Cycles", computed using "llvm-mca" 16.0.6, a cycle-accurate hardware simulator. This allows us to get accurate results on many different architectures without physically possessing them or dealing with timing issues^[9]. The results above are for the "znver4" (i.e. AMD Zen 4) architecture, but results for several others can be downloaded.

An end-to-end latency of the operation would be deeply misleading since other instructions interleave around it, so a standard practice is to repeat the operation many times and take the average (I used 100 times). Additionally, to reduce the influence of register dependencies, I unrolled the operation 4 times and disabled auto-vectorization to prevent optimization between iterations^[10].

Relative Error

The relative error is fairly straightforward. The result of std::sqrt(⋯) is bitwise-perfect, and so is taken as the reference. The relative error of a value is the absolute difference of the value from the reference, divided by the reference. All error calculations were done in double-precision.

There are a few edge cases to consider. A correct result of 0 or ∞ would result in a relative error of NaN, but is defined instead to be 0. A finite result where ∞ is expected is an error of ∞. The relative error of a NaN results in (and is left as) a NaN (though this doesn't show up in the table above; see the raw data).

The approximate methods have different parameters used to tweak the result to be more accurate. It is sensible to tweak the result either to minimize the average relative error or the maximum relative error^[11]. The results above are for coefficient values that minimize the average relative error^[12] (i.e. the middle column below).

Analysis

The built-in (sqrtf_std(⋯)) function is a good place to start. It just calls the sqrtss instruction (as does sqrtf_sse(⋯), which is identical), producing the exact answer for each value. On the Zen 4 architecture, sqrtss has a latency of 15 cycles, but it can issue a new square-root every 5 cycles^[13], meaning that, if we have a lot of square-roots, on average we'll be running three at once and the runtime will be a third as long^[14]. If you don't have a lot of square-roots, then it will be other instructions interleaving. Thus, the pipelined cycles is reported as 5.

The sqrtf_x87(⋯) also gives the exact answers. It is 80-bit, whereas only 32-bit is required, and so it takes commensurately longer. The x87 on the Zen 4 is not pipelined, so it doesn't get the benefit of interleaving, either—however, even with pipelining (e.g. on Skylake), the result is still slower, as we'd expect.

The sqrtf_ssersqrt(⋯) variants launch us into approximate algorithms. The only variant here worth considering is the NR nr0(⋯) one (i.e. no NR at all). Every variant incorrectly gives NaN for 0 and ∞ for subnormals, but at least it is faster (5 → 2 cycles on Zen and 3 → 2 on most Intel architectures). Significantly, it also would auto-vectorize just as well. Still, the fact that it can't handle 0 or subnormals makes it hard to recommend. The NR implementations are not worth it, either; they don't improve the error by much, and they cut into the performance benefit.

The star of the show here is sqrtf_bithack(⋯), which can crank out a square root every 1 cycle on average, on every architecture. With nr0(⋯), the largest relative error for normals is 1.5%, which is pretty good. For 0 and subnormals it doesn't produce a great result, but at least it does produce a small value which is unlikely to cause trouble. The other variant worth considering is nr1_ssercp_fmaB(⋯): a definite slowdown, but the error is much less, improving normals by two orders of magnitude.