Square root approximation:

Some times ago, I asked myself how the square root function might be implemented and approximated. After a quick look on Wikipedia I went on and implemented the Babylonian method:

• Start from a rough estimate y = 0.5f*(1.0f + x)
• Iterate the following: y = 0.5f*(y + t/y)

I will not give details on why it converges, but assuming it does, the limit verifies 2.0*y – y = t/y which would mean y^2 = t if we worked with reals instead of float.

However, the number of iterations needed can be quiet important for small and big numbers.

Another minor drawback is that the convergence point is not always equal to sqrtf(x). There may be a difference of +/-1 on the binary representations of res and sqrtf(x) because of rounding:

For example for x = 1.17549477e-38f we have:

res = 1.08420243e-19f
0.5f*(res + x / res) = res

and

res2 = 1.08420230e-19f
0.5f*(res2 + x / res2) = res

but

x-static_cast<double>(res)*res = -1.4012991325159946e-45
x-static_cast<double>(res2)*res2 = 1.4012982972770227e-45

so res2 is nearer the correct value (which can’t be represented exactly as float) even if res is the fixed point.

To take this into account, I added a final step that does this correction:

double tmp = static_cast<double>(res);

//when x - res*res == 0.0f, the direction to search is unclear if done with float instead
double diff = x - tmp*tmp; 
uint32_t bits = bit_cast<uint32_t>(res);
if (diff < 0)
{
   --bits;
}
else
{
   ++bits;
}

float res2 = bit_cast<float>(bits);
double tmp2 = static_cast<double>(res2);
double diff2 = x - tmp2*tmp2;

if (abs(diff2) < abs(diff))
{
   res = res2;
}

On my machine in the default rounding mode, this gives the exact same result as sqrtf but it is a lot slower!

One trick I found on stackoverflow (by user mvp) is to use frexp and ldexp to move the computation to a smaller interval where it converges faster.

This gives a pretty good approximation for normalized floats and converges in at most three iterations of y = 0.5f*(y + t/y).

However, the following technique found on Wikipedia gives faster implementations (between 1.4 to 2.8 faster):

uint32_t bits = (bit_cast<uint32_t>(x) >> 1) + (1 << 29) - (1 << 22);
res = bit_cast<float>(bits);

The bit trick above followed by 0-4 iterations of the Babylonian method gives the following:

float testSqrt1(float x) noexcept
{
   float res = 0.0f;

   uint32_t bits = (bit_cast<uint32_t>(x) >> 1) + (1 << 29) - (1 << 22);

   res = bit_cast<float>(bits);

   return res;
}

float testSqrt2(float x) noexcept
{
   float res = 0.0f;

   uint32_t bits = (bit_cast<uint32_t>(x) >> 1) + (1 << 29) - (1 << 22);

   res = bit_cast<float>(bits);
   res = (res + x / res) * 0.5f;

   return res;
}

float testSqrt3(float x) noexcept
{
   ...
   
   res = (res + x / res) * 0.5f;
   res = (res + x / res) * 0.5f;

   return res;
}

Where testSqrt4 an testSqrt5 do respectively 3 and 4 steps of res = (res + x / res) * 0.5f.

After some tests, the following gives better results:

uint32_t bits = (bit_cast<uint32_t>(x) >> 1) + 532362861;

Where 532362861 is (1 << 29) - (1 << 22) - 313747 and was found by searching values in the neighbourhood of (1 << 29) - (1 << 22): it is the value that minimize the maximal relative error on all normalized numbers.

This leads to the following corresponding functions:

float testSqrt6(float x) noexcept
{
   float res = 0.0f;
	
   uint32_t bits = (bit_cast<uint32_t>(x) >> 1) + 532362861;

   res = bit_cast<float>(bits);

   return res;
}

float testSqrt7(float x) noexcept {...}
float testSqrt8(float x) noexcept {...}
float testSqrt9(float x) noexcept {...}
float testSqrt10(float x) noexcept {...}

that are similar to testSqrt1, …, testSqrt10 but with the modified constant.

The a simplified version of the micro benchmark I ran was:

float sum = 0;
for (uint32_t i = 8388608; i < 2139095040; ++i)
{
   float val = bit_cast<float>(i);
   float res = testSqrt8(val);
   sum += res;
}
	
cout << "sum:" << sum << endl;

It computes the sum of the square roots of all the normalized float.

This is done because compilers are able to optimize away loops without side effects so you can’t just time the following:

for (uint32_t i = 8388608; i < 2139095040; ++i)
{
   float val = bit_cast<float>(i);
   float res = testSqrt8(val);
}

In general it is a good idea to check the assembly code when doing micro-benchmarks as compilers are quiet smart. The exact code I used is available in the repository.

With clang, I got the following timings:

As you can see, having three iterations of more is too expensive compared to regular sqrtf, while having lower precision as described below.

The following table shows the maximal relative error of the different algorithms on normalized numbers. Given this error bound, the algorithms can be extended to bigger intervals than normalized float, the lower bounds (as bit_cast) are also showed below:

Relative errors become absurdly big on denormalized for all 10 versions (testSqrt10 has 17523.9 for example…), however absolute error stays pretty low.

On denormalized and positive zero: