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.

• ardnew: can you clarify what is meant by You will get the floating point number as first argument after program name and you should implement it. ?
• beary605: What is 0x5f3759df supposed to be in decimal?
• ardnew: @beary605: the value is discussed pretty extensively on the wikipedia article
• null: @ardnew: Clarified that by examples of calling the program.
• z0rberg's: I remember this number. First time I saw this, I was amazed. Wasn't it Michael Abrash's discovery?

Answer 1

Jelly, 41 bytes

l×ばつ+®μ23 ×ばつ"9®ġ’_Hμ‘_×ばつ33_×ばつ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:

l×ばつ+®
l2 Logarithm of the input, base 2
 Ḟ Round down to the integer below (produces the exponent)
 © Store the exponent in a variable
 .* Take (1⁄2 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 1⁄2 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 
 .* 1⁄2 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 231⁄2).

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 ×ばつ y2)), where x2 is half the original input; this is shorter to express as
• y ÷ 2 ×ばつ (3 - (x ×ばつ y2)), where x is the original input.

And here's how it looks in Jelly:

×ばつ33_×ばつH
μ set this point as the new default for missing arguments
 2 square
 ×ばつ3 multiply (×ばつ) by the original input (3)
 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.

• The Empty String Photographer: Uhm, actually, it is defined for negative numbers. It just returns a complex number. (e.g. 1/sqrt(-1) = -i)

Answer 2

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.

• DJMcMayhem: Welcome to the site!
• rafa11111: You can suppress the first line, since, as there is no module being declared in the same source code, to define a 'program' is irrelevant. In the line before 'end' you can change '.5*' for '/2'. It saves 11 bytes ;)
• jvriesem: The people upvoting the Fortran answer are the ones who came looking not for code golf, but for a fast inverse square root routine. ;-)

Answer 3

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 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

• ceased to turn counterclockwis: +1 for the ridiculosity of using a BASIC dialect to golf a performance-related problem.
• Toothbrush: Smaller solution: Type x:q&: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
• SeanC: @toothbrush, symbol for single is !. If I try either type declaration change in Excel, I get an error: Compile Error: Statement invalid inside Type block

Answer 5

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 6

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 7

X86-64 Machine Code, Microsoft Calling Convention - (削除) 65 (削除ここまで) 61 bytes

48 8B 4A 08 mov rcx,qword ptr [rdx+8] 
48 83 EC 28 sub rsp,28h 
E8 C5 FD FF FF call strtof (07FF62D2A1642h)
48 83 C4 28 add rsp,28h 
66 0F 7E C0 movd eax,xmm0 
D9 E8 fld1 
D9 E8 fld1 
D8 C0 fadd st,st(0) 
D8 F9 fdivr st,st(1) 
DC C1 fadd st(1),st 
50 push rax 
D8 0C 24 fmul dword ptr [rsp] 
5A pop rdx 
D1 F8 sar eax,1 
B9 DF 59 37 5F mov ecx,5F3759DFh 
2B C8 sub ecx,eax 
51 push rcx 
D9 04 24 fld dword ptr [rsp] 
D9 C0 fld st(0) 
DE C9 fmulp st(1),st 
DE C9 fmulp st(1),st 
DE E9 fsubp st(1),st 
D8 0C 24 fmul dword ptr [rsp] 
D9 14 24 fst dword ptr [rsp] 
58 pop rax 
C3 ret 

Might be a little shorter in 32 bit. About 20 bytes used in string conversion.

-4 bytes by using strtof

Answer 8

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 9

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 10

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 11

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!

(削除ここまで)

• sergiol: Failed try to outgolf myself
• sergiol: Same 110 byte count: tio.run/…