Transferred from: http://www.guokr.com/post/91389/
After looking at the mathematical principle of 0x5f3759df after spending two days to read the @sheldon to the paper, thanks to @uroboros's explanation, as a person with no computer background, can read (well, actually not all see) This article I am very fulfilling, So write down this post, basically all translation, understand the computer can hold the wrong mentality to see.
Basic Knowledge 1:i>>1
The action i>>1 indicates that the binary number I moves one bit to the right, that is, the last one is deleted and added 0 in the first place
Notice that we delete the last digit of an n-bit decimal number M and then it becomes the n-1 bit, which can be seen as dividing the decimal number by 10 and then rounding down the floor (M/10)
Similarly, when a binary number is deleted, the last one is equivalent to dividing the number by 2 and down to the floor (M/2)
This looks i>>1 like floor (I/2), because the first derivative of the function f (x) =1/sqrt (x) is -1/2* (x) ^-3/2, just a -1/2 in front, can not help feeling 0x5f3759df-(i >> 1) is the first-order Taylor expansion of a function 1/sqrt (x) at a certain point. In fast inverse Square Root There is this paragraph:
Eberly ' s explanation was, it produced linear approximation, but is incorrect; we'll see the guess is piecewise Linear, and the function being approximated is not the same in all cases.
Eberly's explanation is that this is equivalent to a linear approximation, but the explanation is wrong. We will see that this estimate is piecewise linear, and the approximate function is not the same in all cases. ”
Why do you say that? Here's the basics. 2: How to store floating point numbers
Various types of floating-point numbers can be stored by looking at IEEE745 (completely unaware of what it is) to understand
Here is a 32-bit single-precision floating-point number and is always a positive number, expressed in the following way:
0| e| M
where 0 represents a positive number
E is the exponent, e=0 is equivalent to 2^-127
M represents a number with an absolute value less than 1, but notice that one is omitted here. That is, when the conversion bit decimal should be expressed as (1+m)
Then the number after conversion should be: (1+m) 2^ (E-127)
In this way we find that i>>1 is not entirely floor (I/2) but rather a number is changed (Floor (M/2) +1) *2^floor ((E-127)/2)
And depending on the parity of E, the mantissa may also need to be added 1/2
Let e=e-127, notice that 1/sqrt (x) is to make the exponent of the index into-E/2, so Carmack may just want to produce a E/2 and use the displacement, next is to find a subtraction after generating-e/2 and let the mantissa as close as (1+m) ^-1/2
Since this number is necessarily positive, the value is assumed to be:
0| T1: R2
The next step is to discuss the situation:
Assuming that e is even, when the right shift of the exponent does not affect the value of the mantissa:
This time E is an odd number, so that the e=2d+1
is subtracted after the exponential portion becomes:
Note that the subtraction here is actually to convert the floating-point number to Integer (that is, regularization) and then subtract, rather than the normal floating-point number plus minus.
because the initial value must be positive, so we need the formula must be positive, because 0=<e<=256, so r1>=128
because we are talking about the E is even, that is, the last number must be 0, so do not consider the effect of his right shift to M, So the result of subtracting the two numbers is:
Note that this is used M/2 instead of floor (M/2) because the error of this point is too small for the other error
of course, there is a situation, that is r2< M/2, in fact, binary subtraction and decimal almost, if the mantissa is small, then it will be like a higher number of "borrow", that is, the result of subtraction after this situation is:
We define:
Then we can combine these two cases into a function:
This function is our approximation of the function 1/sqrt (x), so our goal is to make the relative error of the function | ( Y-Y0)/y| as small as possible:
so we get an error equation:
:
Note that the 1/8 here is actually made up, the specific method can be assumed to be less than a positive number a, because of the 0<r2<1,0<m<1, we can be expanded by the R1 about a range. Make a as small as possible, making the range less than one. According to the R1 must be an integer character, we can determine the R1 that make the error minimum. Here r1=190, bring it into hex and move right (notice that there is a 0 of the symbol at the beginning) is 0x5f, just the first few of the Black Magic constants.
What happens when E is odd? In fact the same way, if E is odd, then m/2 need to add 1/2 (the first bit of the mantissa is equivalent to 1/2), according to the same method, we can get another relative error function:
which
Interested can calculate how much R1 should be in this case, the author is very lazy to say that because of the need to let the constant in both cases, so let r1=190.
After that is to determine the value of the R2, various sections of discussion, too tangled we will not look at (anyway, finally did not calculate the right, bay stretched out), determine down about 0.45, and then through the software to calculate the optimal solution. It's worth noting that Chris Lomont is using Mathematica 4, and is proud to be a LZ with Mathematica 8, perhaps with more advanced software to get better values?
Well, finally finished, if you find someone else has written similar ... Please don't tell me.
Turn: About black Magic constants