need some help understanding the Fast Inverse Sqrt

This is the place for queries that don't fit in any of the other categories.

need some help understanding the Fast Inverse Sqrt

Postby Tcll » Fri Feb 07, 2014 4:30 am

Code: Select all
# constants:

global infinity; infinity = float('inf')
global minus_infinity; minus_infinity = float('-inf')
global min_float; min_float = 1.1754943508222875e-38 # 0x00800000
global max_float; max_float = 3.4028234663852886e+38 # 0x7F7FFFFF

from math import sqrt
def InverseSqrt(x): # C4 game engine:
    if (x < min_float): return infinity

    r = sqrt(x) # __frsqrte(x)
    r = 0.5 * r * (3.0 - x * r * r)
    r = 0.5 * r * (3.0 - x * r * r)
    return r

from ctypes import cast, pointer, c_int, POINTER, c_float
def InverseSqrt2(x): # http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/
    xhalf = 0.5*x

    i = cast(pointer(c_float(x)), POINTER(c_int)).contents.value
    i = 0x5f3759d5 - (i >> 1)
    x = cast(pointer(c_int(i)), POINTER(c_float)).contents.value
    x = x*(1.5 - xhalf*x*x)
    return x

after the reference in the code, I found this:
http://stackoverflow.com/questions/1778 ... n-modern-c

I just need some comparisons to help me understand...

here's the C++ representation for the C4 engine:
Code: Select all
const float K::infinity =         *(float *) &hex_7F800000;
const float K::min_float =         *(float *) &hex_00800000;

   float C4::InverseSqrt(float x)
   {
      if (x < K::min_float)
      {
         return (K::infinity);
      }

      float r = __frsqrte(x);
      r = 0.5F * r * (3.0F - x * r * r);
      r = 0.5F * r * (3.0F - x * r * r);
      return (r);
   }


I just want to know if my translation is working right, and so far, both return different results :/
I'm not sure if __frsqrte() is defined locally or globally... >_>
(I know 50x more about Python than I do C++ sooo...)
User avatar
Tcll
 
Posts: 100
Joined: Wed Jan 01, 2014 6:36 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby Mekire » Fri Feb 07, 2014 6:09 am

Haven't had time to look at your problem as I'm at work, but just a note; the way you are using the global keyword is doing absolutely nothing in your code. Anything defined in the global scope is global. The global keyword is used for assigning objects to the global scope while inside another scope. In python one indicates that variables are global constants by writing them in ALL_CAPS (and imports should be at the very top of the module). Other python programmers will expect you to adhere to these conventions, and you slow code comprehension by not doing so.

-Mek
User avatar
Mekire
 
Posts: 1011
Joined: Thu Feb 07, 2013 11:33 pm
Location: Amakusa, Japan

Re: need some help understanding the Fast Inverse Sqrt

Postby Tcll » Fri Feb 07, 2014 6:14 am

wow I'm finding nothing on the src or python examples for frsqrte...

I'm hardly finding any C++ examples...
just documentation and mail...

meanwhile I'm stuck with no idea how this function should work >.>

as for testing, I've built a while loop, and function2 returns results almost exact to (1/sqrt(x))...
(where the primary function is way off)

I'm tempted to use that, but then I get hit with the fact that C4 may have a different use for it...
I need to know what the results would be in C++

if they're the same as function2, then I'll use function2

Mekire wrote:Haven't had time to look at your problem as I'm at work, but just a note; the way you are using the global keyword is doing absolutely nothing in your code. Anything defined in the global scope is global. The global keyword is used for assigning objects to the global scope while inside another scope. In python one indicates that variables are global constants by writing them in ALL_CAPS (and imports should be at the very top of the module). Other python programmers will expect you to adhere to these conventions, and you slow code comprehension by not doing so.

-Mek


the globals are copied (modified values) from the C4 src and are meant for an even larger interface

as for imports... I'll do what's needed when I have the src I need
User avatar
Tcll
 
Posts: 100
Joined: Wed Jan 01, 2014 6:36 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby stranac » Fri Feb 07, 2014 8:42 am

That's some interesting code.
It seems to be based on this piece of code I found on wikipedia(with original comments, which should help you understand):
Code: Select all
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
 
   return y;
}


I especially like the wtf comment :P

There are probably some further explanations on the page:
http://en.wikipedia.org/wiki/Fast_inverse_square_root
Friendship is magic!

R.I.P. Tracy M. You will be missed.
User avatar
stranac
 
Posts: 1209
Joined: Thu Feb 07, 2013 3:42 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby Mekire » Fri Feb 07, 2014 9:14 am

Yeah I saw that too. Been trying to implement it as a c extension and failing miserably. I think that is about the only way this will be worth while though. If you don't write it in c, you might as well just use 1/math.sqrt(x).

-Mek
User avatar
Mekire
 
Posts: 1011
Joined: Thu Feb 07, 2013 11:33 pm
Location: Amakusa, Japan

Re: need some help understanding the Fast Inverse Sqrt

Postby Tcll » Fri Feb 07, 2014 12:45 pm

stranac wrote:I especially like the wtf comment :P

lol IKR
I saw that article as well :P

but I need to know if the results of (1/sqrt(x)) and the results of the C4 code match or not.
or if the C4 values I'm getting in my python translation are correct.

so basically, it's a matter of using thing1 or thing2... heh

could I get the printed results from the C++ code that I can test with my python codes plox. :)


EDIT: @Mekire:
use this:
Code: Select all
from ctypes import cast, pointer, c_int, POINTER, c_float
def InverseSqrt2(x): # http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/
    xhalf = 0.5*x

    i = cast(pointer(c_float(x)), POINTER(c_int)).contents.value
    i = 0x5f3759d5 - (i >> 1)
    x = cast(pointer(c_int(i)), POINTER(c_float)).contents.value
    x = x*(1.5 - xhalf*x*x)
    return x

this is the fast version

and the results (almost) match. (as it's supposed to)
User avatar
Tcll
 
Posts: 100
Joined: Wed Jan 01, 2014 6:36 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby Mekire » Fri Feb 07, 2014 1:30 pm

That is not "implement[ing] it as a c extension". If you are concerned with the speed of this operation you need to actually implement it in c. Obfuscating python to do something as simple as 1/sqrt(x) is just ridiculous. If you really need the speed, do it in a way that will give it to you.

Demonstrating:
Code: Select all
import math

from timeit import timeit
from ctypes import cast, pointer, c_int, POINTER, c_float


def InverseSqrt2(x): # http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/
    xhalf = 0.5*x

    i = cast(pointer(c_float(x)), POINTER(c_int)).contents.value
    i = 0x5f3759d5 - (i >> 1)
    x = cast(pointer(c_int(i)), POINTER(c_float)).contents.value
    x = x*(1.5 - xhalf*x*x)
    return x

def invsqrt(x):
    return 1./math.sqrt(x)
Code: Select all
>>> InverseSqrt2(5)
0.44714101924152144
>>> invsqrt(5)
0.4472135954999579
>>> timeit(lambda : InverseSqrt2(5))
8.292665217028755
>>> timeit(lambda : invsqrt(5))
0.5192085575675165
>>>
The ctypes implementation is massively slower than doing it the idiomatic way.

-Mek
User avatar
Mekire
 
Posts: 1011
Joined: Thu Feb 07, 2013 11:33 pm
Location: Amakusa, Japan

Re: need some help understanding the Fast Inverse Sqrt

Postby Tcll » Fri Feb 07, 2014 2:54 pm

oh I see...
alright thanks :)

EDIT:
so what about a few tests from the C4 code??
is it the same as the tests for (1/sqrt(x))??
User avatar
Tcll
 
Posts: 100
Joined: Wed Jan 01, 2014 6:36 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby Mekire » Fri Feb 07, 2014 3:42 pm

Okay, if this interests you, I got my c extension working.

Here is the c extension, cmathext.c (keep in mind I'm garbage at c):
Code: Select all
#include <Python.h>


/*
 * Function to be called from Python
 */
static PyObject* py_invsqrt(PyObject* self, PyObject* args)
{
    const float number;
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    PyArg_ParseTuple(args, "f", &number);

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );
    return Py_BuildValue("f", y);
}


/*
 * Bind Python function names to our C functions
 */
static PyMethodDef module_methods[] = {
    {"invsqrt", py_invsqrt, METH_VARARGS},
    {NULL, NULL}
};


/*
 * Python calls this to let us initialize our module
 */
void initcmathext()
{
    (void) Py_InitModule("cmathext", module_methods);
};

Here is my setup.py file to compile it:
Code: Select all
from distutils.core import setup, Extension

module1 = Extension('cmathext',
                    sources = ['cmathext.c'])

setup (name = 'cmathext',
       version = '1.0',
       description = 'Math functions in C.',
       ext_modules = [module1])

And finally testing it against the two previous methods:
Code: Select all
import math
import cmathext

from timeit import timeit
from ctypes import cast, pointer, c_int, POINTER, c_float


def InverseSqrt2(x):
    xhalf = 0.5*x
    i = cast(pointer(c_float(x)), POINTER(c_int)).contents.value
    i = 0x5f3759d5 - (i >> 1)
    x = cast(pointer(c_int(i)), POINTER(c_float)).contents.value
    x = x*(1.5 - xhalf*x*x)
    return x


def invsqrt(x):
    return 1./math.sqrt(x)


if __name__ == "__main__":
    print("Ctypes: {}".format(timeit(lambda : InverseSqrt2(5))))
    print("Math module: {}".format(timeit(lambda : invsqrt(5))))
    print("Custom C extension: {}".format(timeit(lambda : cmathext.invsqrt(5))))
Code: Select all
Ctypes: 7.85409705254
Math module: 0.496447274005
Custom C extension: 0.325943794404

The custom c extension edges out the standard python method by a little bit. If you are making millions of calls to this function per second it would make a difference.

You will have to test it to see if it gives you the values you want.

Here is the compiled file in case you can't get it to compile:
http://www.mediafire.com/download/e8x3s822m032kil/cmathext.pyd

-Mek
User avatar
Mekire
 
Posts: 1011
Joined: Thu Feb 07, 2013 11:33 pm
Location: Amakusa, Japan

Re: need some help understanding the Fast Inverse Sqrt

Postby Tcll » Fri Feb 07, 2014 4:01 pm

indeed nice :3

now I just need to figure out if the C4 code will output similar results
and I'm even worse than you are at C++, heh

EDIT:
Tcll wrote:
here's the C++ representation for the C4 engine:
Code: Select all
const float K::infinity =         *(float *) &hex_7F800000;
const float K::min_float =         *(float *) &hex_00800000;

   float C4::InverseSqrt(float x)
   {
      if (x < K::min_float)
      {
         return (K::infinity);
      }

      float r = __frsqrte(x);
      r = 0.5F * r * (3.0F - x * r * r);
      r = 0.5F * r * (3.0F - x * r * r);
      return (r);
   }

User avatar
Tcll
 
Posts: 100
Joined: Wed Jan 01, 2014 6:36 pm

Re: need some help understanding the Fast Inverse Sqrt

Postby stranac » Fri Feb 07, 2014 4:03 pm

Interesting...the ctypes code just raises a MemoryError for me...

@C API:
Why u so ugly?
Friendship is magic!

R.I.P. Tracy M. You will be missed.
User avatar
stranac
 
Posts: 1209
Joined: Thu Feb 07, 2013 3:42 pm


Return to General Coding Help

Who is online

Users browsing this forum: liveact, pam, Yoriz and 4 guests