Good seed for the square-root algorithm - Dahlquist's book on Numerical methods

Quasar Chunawala

Hi all,

I want to implement Heron's(Baybylonian) method of computing square root of a floating point number $$a$$ in C++. The number of iterations required to get the desired precision depends on the initial seed $$x_{0}$$. The Heron's method to find the square root of $$a$$ is:

$$x_{n+1}=\frac{1}{2}\large(x_{n}+\frac{a}{x_{n}}\right)$$

Dahlquist's book on numerical methods suggests the following:

1. If we first shift the mantissa so that the exponent becomes even, $$a=c2^{2e}$$ and $$1/2{\le}c{\lt}2$$, then

$$\sqrt{a}=\sqrt{c}\cdot{2^{e}}$$

2. We need only consider the reduced range $$1/2{\le}c{\le}1$$, since for $$1{\lt}c{\le}2$$, we can compute $$\sqrt{1/c}$$ and invert.

3. To find an initial approximation $$x_{0}$$ to start the Newton iterations when $$1/2{\le}c{\lt}1$$, we can use the linear interpolation of $$x=\sqrt{c}$$ between the endpoints $$1/2$$ and $$1$$ giving,

\begin{aligned} x_{0}(c)&=\frac{(1/{\sqrt{2}}(1-c)+1(c-1/2)}{1-1/2}\\ &=\sqrt{2}(1-c)+2(c-1/2) \end{aligned}

In the above steps, I am not clear with respect to a few things.

First, how do I get my floating point number into the form $$a=c2^{2e}$$? I know to find the exponent of a float as:

C++:
float a = 8.5;
int i = *(int*)&a;
int exp = i >> 23;
cout << "Exponent = " << exp;

Second, what is it that we do if $$1{\lt}c{\le}2$$?

And third, what is the intuition behind the interpolation to get an initial estimate of $$x_{0}$$?

Aside: I found this beautiful post explaining how the legendary game programmer John Carmack came up with this evil floating point bit level hack to calculate inverse square roots. The below code snippet is from Quake 3's source code. Here's the explanation - 0x5f3759df

C++:
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}

Daniel Duffy

C++ author, trainer
If the mapping is a contraction, any seed will do!

Fixed-point iteration - Wikipedia

Speed of convergence can be increased by Aitken acceleration.

Your solution looks a wee bit over-engineered.

Another cool solution is least square minimization of

min (x^2 - a)^2 in a given interval. Use e.g. Brent's Method (in Boost C++ libs).

Last edited:

Quasar Chunawala

If the mapping is a contraction, any seed will do!

Fixed-point iteration - Wikipedia

Speed of convergence can be increased by Aitken acceleration.

Your solution looks a wee bit over-engineered.

Another cool solution is least square minimization of

min (x^2 - a)^2 in a given interval. Use e.g. Brent's Method (in Boost C++ libs).

I shall try to write a code using the Brent's method using C++. I registered for the numerical methods and ODE courses on NPTEL(IIT's MOOCs) . The written exams would be conducted on 24th October. The best part is, I also get a certificate. It's kinda nice to work towards an end objective(deadline), it keeps me disciplined.

Replies
2
Views
899
Replies
0
Views
1K
Replies
0
Views
1K
Replies
15
Views
20K
Replies
7
Views
36K