- Joined
- 11/5/14
- Messages
- 295
- Points
- 53
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:
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:
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
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:
- 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}}\)
- 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.
- 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:
Code:
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
Code:
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;
}
