Fast inverse square root

Question

Inspired by http://meta.codegolf.stackexchange.com/questions/547/write-a-code-golf-problem-in-which-script-languages-are-at-a-major-disadvantage question on meta, I've decided to make a question which could be problematic for scripting languages.

The goal is to calculate fast inverse square root, just like it was done in Quake III Arena. You will get the floating point number as first argument after program name and you should implement it. Simply doing ** -0.5 is disallowed as it doesn't implement the algorithm.

Your program will be called like this. The 12.34 could be other value.

$ interpreter program 12.34 # for interpreted languages
$ ./a.out 12.34 # for compiled languages

For comparison, this is original Quake III Arena implementation.

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 heck?
        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;
}

You have to do just one iteration because the second was commented out.

Winning condition: Shortest code.

Answer 1

Perl, 89 85 chars

includes the necessary symbols for declaring and implementing a function

edit: standalone program. receives input parameter as command line argument

$_=unpack'f',pack'i',0x5f3759df-unpack('i',pack'f',$n=shift)/2;
print$_*1.5-$n/2*$_**3

Answer 2

Jelly, 41 bytes, language postdates challenge

l2Ḟ©.*×+®µ23 .*ד9®ġ’_Hµ‘_Ḟ×Ḟ2*$$µ²×³3_×H

Try it online!

I know this challenge was designed in the hope that golfing languages would have a hard time. In a way, they do; Jelly doesn't have any primitives for looking at the representation of a floating point number in memory. However, it's still possible to solve the challenge via working out "manually" what the representation would look like using the basic definitions of floating point arithmetic, and in fact there's some amount of mathematical interest in doing things "the hard way". Jelly's so much terser than (say) C, that the fact it has to do more work doesn't prevent the program being considerably shorter.

Explanation

The input is a floating-point number, as a number. In order to run the fast inverse square root algorithm, we need to find how it would be represented in memory. Jelly doesn't have a way to do that by looking at the bit pattern, but we can do it arithmetically.

First, note that the input must be positive (or its inverse square root would be undefined). As such, it's laid out in memory as follows, from most to least significant:

• A zero bit (the sign bit, zero means nonnegative so this will always be zero);
• The exponent (dividing the number by 2 to the power of the exponent scales it into the range 1 to 2), represented in 8 bits as its value minus 127;
• The mantissa (the result of the above division), represented in 23 bits via subtracting 1, then multiplying by 2²³.

The results of each of these calculations can be represented directly in Jelly. As such, we could generate the same output as C's convert-float-to-int memory hack does like this:

• Calculate the exponent
• Calculate the mantissa
• Take ((exponent + 127) × 2²³) + (mantissa × 2²³ - 1)

However, the last expression shown here simplifies to the following rather simpler form:

• 2²³ × (exponent + mantissa + 126)
• alternatively, 8388608 × (exponent + mantissa) + 1056964608 (the result of multiplying the above out)

Let's call the exponent + mantissa of the input floating point number em for short. (The exponent + mantissa of a floating point number uniquely defines it.) In other words, after the // evil floating point bit level hacking comment, the C program is currently working with i = 8388608 × em + 1056964608.

The next steps in the C program are to halve i, and subtract it from 0x5f3759df (which, in decimal, is 1597463007). i is (8388608 × em + 1056964608); halving it gives (4194304 × em + 528482304); the subtraction gives (1068980703 - 4194304 × em).

Then, C convert this number back to a floating point number y. Let's call the exponent + mantissa of y em'. What the C program is therefore effectively doing in the floating point representational hacking is solving the following equation:

• (1068980703 - 4194304 × em) = (8388608 × em' + 1056964608), which simplifies to:
• 12016095 = 4194304 × em + 8388608 × em'
• em' = (12016095 - 4194304 × em) ÷ 8388608 = (12016095 × 2^-23) - em / 2

Now, to convert this into Jelly. We have a nice arithmetic definition of em and of em', so we can translate it directly. First, here's how to calculate em:

l2Ḟ©.*×+®
l2                 Logarithm of the input, base 2
  Ḟ                Round down to the integer below (produces the exponent)
   ©               Store the exponent in a variable
    .*             Take (½ to the power of the exponent)
      ×            Multiply by {the original value, by default}
       +®          Add the value of the variable (i.e. the exponent)

The original value multiplied by ½ to the power of the exponent is equal to the original value divided by 2 to the power of the exponent, i.e. the mantissa, so this is em.

Getting em' is very straightforward from here, if we want it. However, we're going to need to convert the exponent+mantissa format back into a floating point number. This can be done unambiguously (the exponent's always an integer, the mantissa runs from 1 inclusive to 2 exclusive). However, to extract the exponent, we're going to want to floor and subtract 1. As such, it's actually shorter to generate em' -1.

• em' = (12016095 × 2^-23) - em / 2, so
• em'-1 = (3627487 × 2^-23) - em / 2

We can encode this formula in Jelly directly:

µ23 .*ד9®ġ’_H
µ                  set this point as the new default for missing arguments
 23                restart with 23  
    .*             ½ to the power of that (i.e. 2 to the power -23)
      ×            times
       “9®ġ’       3627487 (compressed representation)
            _H     minus half {the value as of the last µ command}

Incidentally, we need the space because 23. is a single token in Jelly (which evaluates to 23½).

The next step is to convert from em'-1 to an actual floating point number y. We can extract the exponent using Ḟ; then the mantissa is em'-1 + 1 - the exponent. To produce the floating point number, we want to calculate mantissa × 2^exponent:

µ‘_Ḟ×Ḟ2*$$
µ                  set this point as the new default for missing arguments
 ‘                 increment (producing em'-1 + 1)
  _Ḟ               subtract the exponent (producing the mantissa)
    ×              multiply by
     Ḟ2*           2 to the power of the exponent
        $$         parse the previous three links as a group

Finally, we just need to handle the line marked // 1st iteration. This is just regular arithmetic, so encodes into Jelly really easily. The formula is:

• y × (1.5 - (x2 × y²)), where x2 is half the original input; this is shorter to express as
• y ÷ 2 × (3 - (x × y²)), where x is the original input.

And here's how it looks in Jelly:

µ²×³3_×H
µ                  set this point as the new default for missing arguments
 ²                 square
  ׳               multiply (×) by the original input (³)
    3_             subtract (_) from 3
      ×H           multiply by half {the value as of the last µ command}

Verification

Running this program on the input 2 produces the output 0.7069300386983334. This is the same value (allowing for differences in the float-to-string conversion) as produced by this VBA answer, and not equal to the mathematically correct value for 2^-0.5.

Answer 3

FORTRAN, 124 113 Chars

saved 11 Bytes thanks to rafa11111

I know that this question is rather old, but it deserves a FORTRAN answer. It cannot outgolf Perl or Jelly, but it is at least not the worst.

character*9 c
equivalence(r,i)
call getarg(1,c)
read(c,*)y
r=y
i=Z'5f3759df'-ishft(i,-1)
print*,r/2*(3-y*r*r)
end

The hardest thing is to get command line arguments into FORTAN. Fortunately there is getarg() as alias for get_command_argument(). Note the use of an equivalence statement to avoid type casting. r and i are simply the same bytes but can be addressed as different datatypes. Since FORTRAN pointers cannot be casted into another datatypes like in C this is the way to go in FORTRAN.

Answer 4

vba, 171

Type x
q As Long
End Type
Type z
w As Single
End Type
Function q(n)
Dim i As x,y As z
x=n/2
y.w=n
LSet i=y
i.q=&H5F3759DF-i.q/2
LSet y=i
q=y.w*(1.5-x*y.w*y.w)
End Function

finding a way to do the pointer cast was the hardest part. unfortunately, as there is no real direct way, I had to define my own type, which added to the length of the program

call from the immediate window with ?q(number)

result:

?2^-0.5
0.707106781186548
?q(2)
0.706930038698333

Answer 5

C#, 157 chars

Horrible for golfing, but why not...

class C{unsafe static void Main(string[]a){var n=float.Parse(a[0]);var i=*(int*)&n/-2+0x5f3759df;var y=*(float*)&i;System.Console.Write(y*1.5f-y*n/2f*y*y);}}

Readable version:

class C
{
    unsafe static void Main( string[] a )
    {
        var n = float.Parse( a[0] );
        var i = *(int*) &n / -2 + 0x5f3759df;
        var y = *(float*) &i;
        System.Console.Write( y * 1.5f - y * n / 2f * y * y );
    }
}

Answer 6

Julia, 96 characters

F=Float32;\ =reinterpret;n=parse(F,ARGS[]);y=F\(0x5f3759df-UInt\n>>1);print(y*(1.5f0-(n/2*y*y)))

Julia (function), 73 characters

n->(\ =reinterpret;y=Float32\(0x5f3759df-UInt\n>>1);y*=(1.5f0-(n/2*y*y)))

Needs to be run on a 32-bit system for UInt to mean UInt32. Unfortunately, there isn't a similar alias for Float32, and if I tried to just parse(ARGS[]), it would be Float64.

Answer 7

Portable C99, 114 bytes

float g(f,t,e)float f,t;{t=1.93243-frexp(f,&e);t-=.5*(e&1);t=ldexp(t+(t<1),-(e>>1)-(t<1));return(1.5-.5*f*t*t)*t;}

Try it online!

A solution in a more portable subset of C99. Uses frexp and ldexp from math.h instead of bit twiddling. Has very slight rounding errors. Here's a documented version that should be exact:

float my_rsqrt(float f) {
    int e;
    float r = frexp(f, &e);

    // Subtract from magic mantissa.
    int prev_mode = fegetround();
    fesetround(FE_UPWARD);
    r = (0x0.3759dfp1f + 1.5f) - r;
    fesetround(prev_mode);

    // Handle LSB of exponent that spills into mantissa.
    if (e & 1) r -= 0.5f;

    // Handle overflow.
    if (r < 1.0f) {
        r += 1.0f;
        e += 2;
    }

    // Scale with adjusted exponent.
    r = ldexp(r, -(e>>1));

    return r * (1.5f - 0.5f * f * r * r);
}

Answer 8

Java, 175 bytes

^{_{I decided to give this question a Java answer since it didn't have one, and then it got tweeted and bumped while I was golfing... weird.}}

As a lambda, 101 bytes:

f=n->{float y=Float.intBitsToFloat(0x5f3759df-(Float.floatToIntBits(n)>>1));return y*(1.5f-n/2*y*y);}

As a full program, 175 bytes:

class Q{public static void main(String[]a){float n=Float.parseFloat(a[0]),y=Float.intBitsToFloat(0x5f3759df-(Float.floatToIntBits(n)>>1));System.out.print(y*(1.5f-n/2*y*y));}}

_{^{Darn it, two bytes more than portable C. It's the darn Float.floatToIntBits()&Float.intBitsToFloat(), for Java doesn't like it when you try to say "You know this float that I just defined? Make it an int. You know that int I just defined? Make it a float."}}

Answer 9

C, 173 chars

Might as well put code golfed version of example above. Could be less if I would ignore portability.

#include<stdio.h>
#include<stdint.h>
main(int _,char**a){float n,g;sscanf(a[1],"%f",&n);int32_t o=0x5f3759df-(*(long*)&n>>1);g=*(float*)&o;printf("%g\n",g*(1.5-(n/2*g*g)));}

Answer 10

Tcl, 110 bytes

rename binary B
proc R n {B s [B f f $n] i i
B s [B f i [expr 0x5f3759df-($i>>1)]] f y
expr $y*1.5-$n/2*$y**3}

Try it online!

---

Tcl, 112 bytes

rename binary B
proc R n {B s [B f f $n] i i
B s [B f i [expr 0x5f3759df-($i>>1)]] f y
expr $y*(1.5-$n/2*$y*$y)}

Try it online!

---

Tcl, 113 bytes

rename binary B
proc R n {B s [B f f $n] i i
B s [B f i [expr 0x5f3759df-($i>>1)]] f y
expr $y*(1.5-$n*.5*$y*$y)}

Try it online!

---

Tcl, 129 bytes

rename binary B
proc R n {B scan [B format f $n] i i
B scan [B format i [expr 0x5f3759df-($i>>1)]] f y
expr $y*(1.5-$n*.5*$y*$y)}

Try it online!

Tcl, 133 bytes

proc R n {binary scan [binary format f $n] i i
binary scan [binary format i [expr 0x5f3759df-($i>>1)]] f y
expr $y*(1.5-$n*.5*$y*$y)}

Try it online!

Tcl, 142 bytes

proc R n {binary scan [binary format f $n] i i
set i [expr 0x5f3759df-($i>>1)]
binary scan [binary format i $i] f y
expr $y*(1.5-$n*.5*$y*$y)}

Try it online!