Detailed explanation of password Based on John Carmark

Someone found such a piece of code in the source code of Quake III to calculate the square root:

/* ======================== 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 * (f-(x * y ));
Return number * y;
}

X5f3759df? What is this? After learning the numerical analysis, we know that the square root of the algorithm generally uses
It's an infinitely approximate method, such as the Newton Iteration Method. Sorry, I was so bad at numerical analysis that I could not tell clearly.
. Simply put, for example, finding the square root of 5 and selecting a guess value such as 2, we can calculate this.

/2 = 2.5; 2.5 + 2/2 = 2.25; 5/2. 25 = xxx; 2.25 + xxx/2 = xxxx...
In this way, the result will surely converge to sqrt (5). That's right. The square root is calculated in this way.
. The difference is that he chose a mysterious guess value 0x5f3759df as the starting point, making
The convergence speed of the entire approximation process is soaring. For Quake III's requirement of accuracy 10, the negative cubic power only needs one
The results can be obtained through the next iteration.

Well, if this is not a good B, let's take a look.

Chris Lomont, a mathematician at the Pudu University, thought it was interesting after reading it. He decided to study what camark got.
What are the mysteries of this guess. Lomont is also a cool man. After careful research, it is also theoretically derived from
The best guess value, which is very close to the number of the card mark, 0x5f37642f. Kamakok, is he an alien?

The legend is not over here. Lomont is very satisfied with the calculation result, so he calculates the start point by himself.
The value matches the mysterious number of camark to see who can obtain the square root faster and more accurately. The result is
Mark won... no one knows how Mark found this number.

Finally, Lomont got angry. He tried a number and a number using the brute force method and finally found a card mark.
Although the results produced by these two numbers are very similar, this violent
The force output is 0x5f375a86.

For this reason, Lomont wrote the next paper "Fast Inverse Square Root ".

I used C # To rewrite this function:

Copy codeThe Code is 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 * (f-(x * y ));
Return number * y;
}
}
}

Lines 32nd and 33 use two Newton Iteration Methods to achieve certain accuracy. Of course, you can also control the accuracy by yourself and find the reciprocal of the square root of y, therefore, number * y is returned.

The most critical sentence of the SquareRootFloat function is I = 0x5f3759df-(I> 1 );
Some of the explanations are as follows:

The most important part of the Newton iteration method is to estimate the first approximate root. If the approximate root is close enough to the actual root, a satisfactory solution can be obtained only after a few iterations.

Next, we will try to estimate the first approximate root. This is also the most amazing part of the above functions. It calculates an approximate root which is very close to the real root through some method. Therefore, it only needs to use an iteration process to obtain a satisfactory solution. How does it do it? All the mysteries of this line are:

I = 0x5f3759df-(I> 1); // calculate the first approximate Root

Super inexplicable statement, isn't it? However, if you think about it, you can understand that float data is represented in a 32-bit system.

Bits: 31 30... 031: Symbol bit 30-23: a total of 8 bits, storage index (E) 22-0: a total of 23 bits, storage ending number (M)

Therefore, the 32-bit floating point number is represented by a decimal real number: M * 2 ^ E. Start the root and then the reciprocal is: M ^ (-1/2) * 2 ^ (-E/2 ). Now it is very clear. Statement I> 1: divide the index by 2 to implement the 2 ^ (E/2) Part. Instead, subtract it with a constant to get the M ^ (1/2) and reverse all the exponent symbols at the same time.