Based on John Carmark password detailed _c language

Someone in the source code of Quake III found such a paragraph to be used to square root codes:

/*================squarerootfloat================*/

float Squarerootfloat (float number) {
Long I;
float x, y;
const float F = 1.5F;
x = number * 0.5F;
y = number;
i = * (long *) &y;
i = 0x5f3759df-(i >> 1); Pay attention to this line
y = * (float *) &i;
y = y * (f-(x * y * y));
y = y * (f-(x * y * y));
return number * y;
}

X5f3759df? What is this thing? Studied the numerical analysis to know that the algorithm inside the square root is generally used
is a method of infinite approximation, such as Newton's iterative method, I'm sorry I was too bad in numerical analysis, I can't tell
。 In simple terms, for example, ask for the square root of 5, choose a guess value such as 2, then we can calculate

/2 = 2.5; 2.5+2/2 = 2.25; 5/2.25 = xxx; 2.25+XXX/2 = xxxx ...
This repeated iteration, the result must converge to sqrt (5), yes, the general square root is so calculated
。 And the difference between carmack is that he chose a mysterious guessing value 0x5f3759df as the starting point, making
The convergence rate of the whole approximation process is soaring, and for the 10 minus three times of the precision required by Quake III, only one
Iterations to get results.

Well, if that's not bull B, then look.

Chris Lomont, a mathematician at Purdue University, found it interesting and decided to study the Carmack.
What is the mystery of this conjecture. Lomont is also a cow, after careful study of the theory also deduced a
Best guesses, and Carmack numbers are very close, 0x5f37642f. Carmack, is he an alien?

The legend does not end here. Lomont to calculate the result is very satisfied, so take their own calculation of the starting
The value and Carmack of the mysterious number of games to see whose numbers can be faster and more accurate to obtain the square root. The result is
Carmack won ... No one knows how Carmack found the number.

Finally Lomont anger, the use of violent methods a number of a number try to come, finally found a more than Carmack number
Word is better than a little bit of numbers, although in fact these two numbers produce a very approximate result, this storm
The number of forces is 0x5f375a86.

Lomont wrote a paper for this, "Fast inverse Square Root".

I use C # for this function to rewrite:

Copy Code code as follows:

Using System;
Using System.Collections.Generic;
Using System.Text;

Namespace ConsoleApplication1
{
Class Program
{
static void Main (string[] args)
{
Console.WriteLine ("Carmark ' s Method:");
Console.WriteLine (Squarerootfloat (3.0f). ToString ());
Console.WriteLine ("Use Math.sqrt () method:");
Console.WriteLine (((float) math.sqrt (3.0)). ToString ());
Console.read ();
}

private static float Squarerootfloat (float number)
{

Long I;
float x, y;
const float F = 1.5F;
x = number * 0.5F;
y = number;
Unsafe
{
i = * (long *) &y;
i = 0x5f3759df-(i >> 1); Pay attention to this line
y = * (float *) &i;
}
y = y * (f-(x * y * y));
y = y * (f-(x * y * y));
return number * y;
}
}
}

32nd, 33 lines used two Newton iterative method to achieve a certain degree of precision, of course, you can also control the accuracy of the square root of Y to find the reciprocal, so finally returned to Number*y.

The most critical sentence of the Squarerootfloat function is i=0x5f3759df-(i>>1);
The following is a partial explanation of it:

The most important part of Newton's iterative method is to estimate the first approximate root. If the approximate root is close enough to the true root, then only a few iterations are needed to get a satisfactory solution.

Then we try to estimate the first approximate root. This is also the most magical place on the function above. It calculates an approximate root that is very close to the real root by some means, so it only needs to use an iterative process to obtain a satisfactory solution. How does it do that? All that is in this line:

i = 0x5f3759df-(i >> 1); Compute the first approximate root

Super inexplicable statement, isn't it? But if you think about it, it's understandable: float-type data is represented on a 32-bit system.

Bits:31 30 ... 031: Sign bit 30-23: Total 8 digits, save index (E) 22-0: Total 23 digits, save Mantissa (M)

So, the 32-bit floating-point number in decimal notation is: m*2^e. Open the root and then the countdown is: m^ ( -1/2) *2^ (-E/2). It's very clear now. Statement i>>1 its job is to divide the exponent by 2 to achieve the 2^ (E/2) part. And the front is subtracting it with a constant to get m^ (1/2) to invert all indices at the same time.