Q_rsqrt.py

import ctypes

# quakes "fast" inverse square root
def Q_rsqrt(value):

    # base values
    c_value = ctypes.c_float(value)
    c_value2 = ctypes.c_float(value * 0.5)

    # cast to long so that we can bit manipulate the number
    c_pLongValue = ctypes.cast(ctypes.pointer(c_value), ctypes.POINTER(ctypes.c_long))
    c_longValue = ctypes.c_long(c_pLongValue[0])

    # the magic happens here
    magic_number = 0x5f3759df

    longShift = (c_longValue.value >> 1)

    #print("{:0032b} orignal number in bit form".format(c_longValue.value))
    #print("{:0032b} right shifted once".format(longShift))
    #print("{:0032b} magic number".format(magic_number))
    
    # this computes a very good starting point that is usually 3.4% wrong
    c_newLongValue = ctypes.c_long(magic_number - longShift)

    # cast back to float
    c_pNewFloatValue = ctypes.cast(ctypes.pointer(c_newLongValue), ctypes.POINTER(ctypes.c_float))
    c_value = ctypes.c_float(c_pNewFloatValue[0])

    #print("{:0032b} first guess ({})".format(c_newLongValue.value, c_value.value))

    # runs 1 round of newtons method to lower inaccuracy to 0.17%
    c_value = ctypes.c_float( c_value.value * (1.5 - (c_value2.value * c_value.value * c_value.value )))

    return c_value.value

def test():
    # accuracy as a percentage
    testAccuracy = [0] * 10000
    overall_accuracy = 0

    x = 1
    idx = 0
    while(x <= 100):
        real = x**(-1/2)
        fast = Q_rsqrt(x)
        accuracy = (fast/real)*100
        testAccuracy[idx] = accuracy
        print("index: {} x: {} has accuracy of {}%".format(idx, x, accuracy))
        idx = idx + 1
        x = x + 0.01
    
    totalAccuracy = 0
    total = len(testAccuracy)
    for acc in testAccuracy:
        totalAccuracy += acc
    overall_accuracy = totalAccuracy / total

    print("overall accuracy is {}%".format(overall_accuracy))

test()

# real inverse square root
#x = 5.5**(-1/2)
#print("real inverse square root of 5.5 is {}".format(x))

#y = Q_rsqrt(5.5)
#print("fast inverse square root of 5.5 is {}".format(y))

#print("difference is {} or {}% accurate and {}% inaccurate".format(x - y, (y/x)*100, 100 - ((y/x)*100)))