Fast inverse square root

Jelly, 41 bytes

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.

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.

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