in 2005 game company its software open
sourced the engine for their video game
quake 3 arena
in that source code fans of the game
discovered an algorithm that was so
ingenious it quickly became famous
and the only thing this algorithm does
is calculate the inverse of a square
root
if i had to write a piece of code that
would calculate the inverse of a square
root
this is how i would do it here i'm using
the c
programming language the same
programming language used for quake 3.
but to be fair i wouldn't actually write
the square root in there myself
people work more closely with the c
language or we see but design have
already figured out how to calculate the
square root
and the provided algorithm to us in a
math.h file
that we programmers can then just
include in our program
so what could possibly be so interesting
about the quake 3 algorithm
how does its software calculate inverse
square roots
at first glance it doesn't seem to make
any sense where does this number
0x5f3759df come from
what does this have to do with taking
square roots and why is there a
disgusting curse word in the second
comment
in this video i will show you how with
some cool bit manipulation
you can take in for square roots and the
algorithm that does this
bears the name the fast inverse square
root
first of all why would the game engine
want to calculate 1 divided by square
root of x
if you want to implement physics and
lighting or reflections in your game
engine
it helps if the vectors you're
calculating with are normalized to have
length 1
because otherwise your vectors might be
too short or too long and when you do
physics with them
things can go wrong
as all of you know the length of a
vector is square root of x squared plus
y squared plus z
squared if you don't happen to notice i
claim you've seen this for two
dimensions
it's just pythagoras theorem so if we
want to normalize the vector's length to
one
we have to scale everything down by the
length of the vector
i mean obviously because if we divide
the length of the vector by
the length of the vector we obviously
get one
so all that's left to do is to divide x
y and z by the length
or similarly multiply by one divided by
the length
you might already see where this is
going
calculating x squared plus y squared
plus z squared is easy
and more importantly really fast in code
you would implement it as x times
x plus y times y plus z times z and all
it is is just three multiplications
summed up additions and multiplications
are common operations
that have been designed to be very fast
the square root on the other hand
is a terribly slow operation and
division is not much better
this is not good if we have several
thousands of surfaces
where each has a vector that needs to be
normalized
but this also means that here we have an
opportunity for speed improvements
if we can find even just an
approximation of one divided by square
root of x
as long as it's fast we can save
precious time
the fast inverse square root is such an
approximation
with only an error of at most one
percent while being
three times as fast
looking at the code again we can see
that the beginning is pretty harmless
we are given a number called number as
input
the number we're supposed to take the
inverse square root of
first with the variable i we declare a
32-bit number
then we declare two 32-bit decimal
numbers
x2 and y and then we store 1.5 into the
variable
with the obvious name three halves
the next two lines simply copy half of
the input into x2
and the whole input into y but it's
after that where the magic happens
take a moment to look at it again
well the longer you look at it the less
it makes sense
and the comments on the right are not
really helpful either
but they do hint that there are three
steps to this algorithm
puzzling together these three steps will
show us the brilliancy of this algorithm
but before we start with these three
steps let's first take a look at binary
numbers
we said that in the first line we
declare a 32-bit integer
in a c programming language called a
long that means we're given 32 bits
and we can represent the number with it
but i think you all know how to do that
this is one two three four and so on
up to around two billion but in the next
line we declare two decimal numbers
in c called a float again we're given 32
bits
and we have to represent the decimal
number with it
how would you do that if you and i were
designing decimal numbers
this is probably one way we would do it
just put a decimal point in the middle
in front of the decimal point we count
in the usual way
1 2 3 4 and so on and after the decimal
point there are no surprises either
just remind yourself that this is binary
so instead of tenths hundredths
and thousands we have halves fourths
eights sixteens and any combination of
them
like a half and a fourth give you
three-fourths also known as 0.75
but this idea is actually terrible we've
decimated the range of numbers we can
represent
before we could represent numbers to
around 2 billion
now only to about 32 000.
luckily people much smarter than us have
found a better way to make use of those
32 bits they took inspiration from
scientific notation
the same way we can systematically
represent numbers like 23
000 as 2.3 times 10 to the 4 and .0034
as 3.4 times 10 to the minus 3
we can also represent them in a binary
system where here for example
1 1 0 0 0 could for example be 1.1 times
2 to the 4.
the standard they came up with takes the
name ieee
754
ieee standard 754 defines the following
we are as usual given 32 bits the first
bit is the sign bit
if it is 0 the number is positive if it
is 1 the number is negative
but the numbers quake 3 provides the
fast inverse square root
are always positive i mean obviously
they're positive
if we would have to calculate 1 divided
by square root of -5
something definitely has gone wrong so
for the rest of this video
we ignore the sign bit as it is always
0.
then the next 8 bits define the exponent
that means 2 to the 1 2 to the 2 2 to
the 3 2 to the 4 and so on
with 8 bits we can represent numbers
between 0 and 255
but that's not exactly what we need we
also want negative exponents
this is why everything is actually
shifted down by 127.
so instead of 2 to the 4 we actually
have 2 to the 4 minus 127
if we actually want the exponent to be 4
the bits need to be set to 131
because 131 minus 127
is 4. the last 23 bits
are the mantissa as usual in scientific
notation
we want to denote one digit followed by
the comma
followed by the decimal places but with
23 bits
we can represent numbers from 0 to but
not including
2 to the 23. again that's not exactly
what we need
for scientific notation we need the
mantissa to go from 1 to 10
or in binary scientific notation to go
from one to two
so we could do something that we've
already done before
put a comma after the first bit this
automatically gives us numbers from one
to two
but this naive approach is wasteful you
see the people that designed standard
754
realized that in binary something
happens that happens in no other base
look at the first digit in scientific
notation
the first digit is by definition always
non-zero
but in binary there is only one digit
that is not zero
one and if we know that the first digit
will always be a one
there is no need to store it thus we can
save one bit
by moving the comma one digit to the
left and fixing an extra one in the
number it represents
now our mantissa is between one and two
even though 23 bits gave us numbers
between 0 and 2 to the 23
we scaled them down to get numbers
between 0 and 1
and then we added an extra 1 to get
numbers between 1
and 2. and this already is the main part
of ieee standard 754
but just the so-called normalized
numbers
the informed viewer knows that the
standard also includes denormalized
numbers
not a number infinities and two zeros
but we won't go into those because in
quake 3 it just happens that these are
never inputs into our algorithm
otherwise something definitely has gone
wrong anyway
at no point should our game engine have
to normalize a vector
with infinite length for this algorithm
and for the rest of this video
it will be useful to think of the
mantissa and exponent as the binary
numbers they are
if we are given two numbers one being
the mantissa and 1 being the exponent
23 bits and 8 bits respectively we can
get a bit representation with
2 to the 23 times e plus m if you think
about it
because multiplying e by 2 to the 23
just
shifts e by 23 digits
so that's how one could write the bits
but we get the actual
number behind the bits with this formula
this should seem familiar to you here we
have the exponent with 127 subtracted
from it
and here we have the mantissa with the
extra one in front
but now something completely different
for no obvious reason at
all let's take the logarithm of that
expression
since we're doing computer science here
we take the logarithm base 2.
we simplify as much as we can take out
the exponent
but then we get stuck but not so the
creators of quake
developer gary tarouli knew a trick to
get rid of the logarithm
you see the trick is an approximation to
log of 1 plus x
for small values of x log of 1 plus x is
approximately equal to
x if you think about it this
approximation is actually correct
for x equals zero and x equals one
but we'll add an additional term mu this
correction term can be chosen freely
again with mu equal to zero this
approximation is correct at zero
and one but it turns out that setting mu
to this number
gives the smallest error on average for
numbers between zero and one
so going back to our formula we apply
our trick as
m divided by 2 to the 23 is indeed a
value between 0
and 1. we rearrange a little bit more
and we finally see why we did all those
calculations
m plus e times 2 to the 23 appears
that's our bit representation so let's
think
about what we just did we applied the
logarithm to our formula
and got the bit representation just
scaled and shifted by some constants
so in some sense the bit representation
of a number is its own logarithm
armed with this knowledge we can finally
start with the three steps of the fast
inverse square root
the first step is actually not
complicated it just looks complicated
because it's memory address trickery
so we stored our number into y and now
we want to do cool bit manipulation
tricks
floats unfortunately don't come with the
tools we need to do bit manipulation
the reason you cannot do bit
manipulation on floats is that they were
never designed to do so flows are
inherently tied to the ieee standard 754
longs on the other hand were designed to
do bit manipulation on them
for example here's one trick
bit shifting along to the left doubles
it
and bit shifting it to the right halves
it
and yes if your number is odd you do end
up rounding
but hey we're willing to accept such
inaccuracies as long as this means that
our algorithm is fast
c as pretty much all programming
languages do provides a way to convert
from a float to a long
this conversion does what most
programmers needed to do namely convert
from a decimal number to an ordinary
integer as best as it can do
so if we give it a float like 3.33 it
converts it to an integer
in this case 3 but this is not the
conversion we need here
first of all we don't care about the
resulting integer number
we want to somehow keep our float and
secondly
the bits that lie behind our number get
all messed up
we don't want this conversion to mess
with our bits the only thing we want to
do
is to put the bits one to one into a
long
the way you achieve this is to convert
the memory address
not the number first
we get the address of y this is the
address of a float
then you convert that address from a
floats address to a long's address
the address itself doesn't change but c
now thinks that the number living at
that address is now along
so then you read what's written at that
address because
c now thinks that this is an address of
a long it will read the number at that
address as if it were long
like this we tricked c by lifting the
conversion away from the number itself
to the address of that number and that's
how we get the bits of a number into i
i don't know what else to say that's
just how c works so let's just go to the
next step
the intuition behind the second step is
the following
remind yourself bit shifting a number to
the left doubles it
and shifting it to the right halves it
but what would happen if we did
something like this to an exponent
doubling an exponent squares the number
and halving the exponent gives us the
square root
but now also negating the exponent gives
us 1 divided by square root of x
that's exactly what we need so let's
remind ourselves what our goal is here
we have our numbers stored into y and
our goal was to calculate 1 divided by
square root of a y
as i've already said calculating this
directly is too hard and too expensive
but we've extracted the bits from y and
we've seen with the ieee standard 754
that the bits of a number are in some
sense its own logarithm
that means in i we have stored log of y
up to some scaling and shifting i claim
that our problem becomes
way easier if we work with logs instead
of trying so hard to calculate 1 divided
by square root of y
we instead calculate log of 1 divided by
square root of y
we rewrite this to log of y to the power
of minus a half
so we can take out the exponent
calculating this
is stupidly easy you might think oh no
we have a division in there
didn't you say in the beginning that
divisions are slow well yes but remember
that we can do bit shifts now
instead of dividing by 2 we just bit
shift once to the right
this already explains why we do minus i
bit shifted the ones to the right
but why is this number 0x5f37590f here
again
well because our logarithm is actually
scaled and shifted
so let's calculate and understand where
it comes from
let gamma be our solution then we know
that log of gamma equals to log of y to
the power minus a half
which equals to minus a half times log
of y
now we replace the logarithm with the
bit representation
and then we just solve for the bits of
gamma i'll spare us the details but this
is the result
the magic number turns out to be the
remnants of the error term mu
the scaling factor and the shifting
now we have the bits of resolution and
we can just reverse the steps from the
evil bit hack
to get back the actual solution from
those bits
well actually not the exact solution
just an approximation
this is why we need the third step
after the previous step we have a pretty
decent approximation
but we did pick up some error terms here
and there but thanks to newton's method
we can make a really good
approximation out of a decent one
newton's method is a technique that
finds a root for a given function
meaning it finds an x for which f x
equals zero
it does so by taking an approximation
and returning a better approximation
and usually you repeat this process
until you're close enough to the actual
solution
but it turns out that here we are
already close enough to the actual
solution
that one iteration suffices to get an
error within one percent
the only things newton's method needs is
the function and its derivative
and what newton's method does is that it
takes an x value
and tries to guess by how much it is off
from being a root
it does so by calculating f of x and its
derivative
we can write f of x as y and the
derivative as
d y over dx
we have the ratio between y and the x
offset
and y itself so to get the x offset
we just divide y by the ratio
so then we simply subtract this offset
to get our new
x the informed viewer can now verify
that the last line is one such newton
iteration
applied to the function f of y equals 1
divided by y
squared minus x notice that y being a
root of this function is
equivalent to y being the inverse square
root of x
i really encourage you to verify this
last line of code
since it's really surprising that even
though both the function
and newton's method have a division in
them the code does not
which means that our algorithm is and
stays fast
now we finally understand the fast
inverse square root
it only took us the knowledge of the
ieee standard 754
a trick to outsmart the c programming
language magic bit operations
and the calculus behind newton's method
you