MikeDSutton's picture

Fast inverse square root fix

Here's a fixed version of the double overload for the fast inverse square root method from OpenTK.MathHelper

public static double InverseSqrtFastB(double x) {
    unsafe {
        double xhalf = 0.5 * x;
        long i = *(long*)&x;               // Read bits as integer.
        i = 0x5fe6ec85e7de30da - (i >> 1); // Make an initial guess for Newton-Raphson approximation
        x = *(double*)&i;                  // Convert bits back to float
        x = x * (1.5 - xhalf * x * x);     // Perform left single Newton-Raphson step.
        return x;
    }
}

It's using a 64-bit 'magic' number from "Fast inverse square root" by Chris Lomont (the same paper the float version's magic number comes from.)
Although this version works without casting the inputs and outputs, it's slower than the equivalent using the standard System.Math library:

public static double InverseSqrt(double x) {
    return 1.0 / Math.Sqrt(x);
}

Slightly better performance and accuracy can be attained using a second Newton-Raphson step with the float version:

public static float InverseSqrtFast(float x) {
    unsafe {
        float xhalf = 0.5f * x;
        int i = *(int*)&x;              // Read bits as integer.
        i = 0x5f375a86 - (i >> 1);      // Make an initial guess for Newton-Raphson approximation
        x = *(float*)&i;                // Convert bits back to float
        x = x * (1.5f - xhalf * x * x); // Perform left single Newton-Raphson step.
        x = x * (1.5f - xhalf * x * x); // Perform left single Newton-Raphson step.
        return x;
    }
}