itchyny/fastinvsqrt

Fast inverse square root in programming languages

Fast inverse square root in programming languages

This is a repository for my challenge of writing Fast inverse square root algorithm in many languages.

Writing one algorithm in many languages is fun. I wrote some codes in languages I have never experienced. I learned the differences and similarities between the languages, how languages are influenced by others.

This challenge is just a passing point on my way. Now I want to continue learning some languages I've never thought of writing before this challenge. For example, I have never experienced Rust, OCaml and Erlang but now I think I should continue learning these languages.

Why did I choose this algorithm?

Here is an implementation in C.

float fastInvSqrt(float x) {
  int i = *(int*)&x;
  i = 0x5f3759df - (i >> 1);
  float y = *(float*)&i;
  return y * (1.5F - 0.5F * x * y * y);
}

The pointer casting magic is the most important part of this algorithm. I do not explain why the code works here, please refer to the other references link.

First of all, the algorithm to choose should not be too easy and should not be too difficult. If the problem was easy, I could write the codes after some glance on the language tutorials. On the other hand, I could give up the challenge if it was too difficult.

Secondary, pointer casting is not such an easy task in some languages. Some languages do not allow us to touch the pointer addresses of the variables. For scripting languages, I have to look through the document whether I can get the byte pattern of floating point numbers.

Thirdly, I like the algorithm. The algorithm is largely based on how the floating point number is represented on our computers. When I came across this algorithm, I got very surprised that it effectively uses the casting pointers between different types of number.

Some examples in languages

All the codes are included in the repository but I want to introduce some interesting implementations here.

Groovy

Groovy is like an un-fancy Java.

Scanner scanner = new Scanner(System.in)
while (scanner.hasNext()) {
  try {
    println fastInvSqrt(Float.parseFloat(scanner.nextLine()))
  } catch (NumberFormatException e) {}
}

def fastInvSqrt(float x) {
  int i = Float.floatToRawIntBits(x)
  float y = Float.intBitsToFloat(0x5f3759df - (i >> 1))
  y * (1.5F - 0.5F * x * y * y)
}

Clojure

Clojure is also a JVM language so the Float/floatToRawIntBits function routes in the function of Java.

(defn fast-inv-sqrt [x]
  (let [i (Float/floatToRawIntBits x)
        y (Float/intBitsToFloat (- 0x5f3759df (bit-shift-right i 1)))]
    (* y (- 1.5 (* 0.5 x y y)))))

(doseq [line (line-seq (java.io.BufferedReader. *in*))]
  (try (println (fast-inv-sqrt (read-string line)))
      (catch java.lang.RuntimeException e ())))

Crystal

This language looks like Ruby, but is one of the compiled languages. The code of pointer casting looks pretty.

def fastinvsqrt(x : Float32) : Float32
  i = pointerof(x).as(Int32*).value
  i = 0x5f3759df - (i >> 1)
  y = pointerof(i).as(Float32*).value
  y * (1.5 - 0.5 * x * y * y)
end

while line = gets
  puts fastinvsqrt(line.to_f32)
end

Rust

There are f32::to_bits, f32::from_bits to convert the number types keeping the bits.

use std::io::BufRead;

fn fast_inv_sqrt(x: f32) -> f32 {
    let i = x.to_bits();
    let j = 0x5f3759df - (i >> 1);
    let y = f32::from_bits(j);
    y * (1.5 - 0.5 * x * y * y)
}

fn main() {
    let stdin = std::io::stdin();
    for line in stdin.lock().lines().filter_map(|x| x.ok()) {
        if let Ok(x) = line.parse() {
            println!("{}", fast_inv_sqrt(x));
        }
    }
}

Perl

The pack and unpack is useful in this challenge for scripting languages (PHP, Ruby, Python etc). Somehow this code runs very fast as compiled languages.

use strict;
use warnings;
use feature qw(say);

while (<>) {
  say fastInvSqrt($_);
}

sub fastInvSqrt {
  my ($x) = @_;
  my $i = unpack("l", pack("f", $x));
  $i = 0x5f3759df - ($i >> 1);
  my $y = unpack("f", pack("l", $i));
  $y * (1.5 - 0.5 * $x * $y * $y);
}

JavaScript

I thought it is impossible to implement in this language. There's no difference between float and integer, 32-bit and 64-bit. However, JavaScript has TypedArray, which we can read and write the numbers in specific bit widths.

require('readline').createInterface({
    input: process.stdin,
}).on('line', (line) => {
    console.log(fastInvSqrt(parseFloat(line)));
});

function fastInvSqrt(x) {
    const i = floatToUInt32(x);
    const j = 0x5f3759df - (i >> 1);
    const y = uint32ToFloat(j);
    return y * (1.5 - 0.5 * x * y * y);
}

function floatToUInt32(x) {
    const buf = Float32Array.of(x).buffer;
    return new Uint32Array(buf)[0];
}

function uint32ToFloat(i) {
    const buf = Uint32Array.of(i).buffer;
    return new Float32Array(buf)[0];
}

J

Just difficult to write, difficult to read.

fastInvSqrt =: 3 : 0
  i =. _2 ic 1 fc y
  i =. (dfh'5f3759df') - _1 (33 b.) i
  z =. _1 fc 2 ic i
  z * (1.5 - 0.5 * y * z * z)
)

echo & fastInvSqrt & (0 & ".) & (-. & CRLF) ;. 2 &. stdin ''
exit ''

F#

F# is a .NET language so looking through the MSDN document gave me the way to convert to/from the byte arrays.

let fastinvsqrt(x: float32): float32 =
    let i = System.BitConverter.ToInt32(System.BitConverter.GetBytes(x), 0)
    let j = 0x5f3759df - (i >>> 1)
    let y = System.BitConverter.ToSingle(System.BitConverter.GetBytes(j), 0)
    y * (1.5F - 0.5F * x * y * y)

seq { while true do yield stdin.ReadLine() }
    |> Seq.takeWhile ((<>) null)
    |> Seq.iter (fun x -> try float32 x |> fastinvsqrt |> printfn "%f"
                          with :? System.FormatException -> ())

Haskell

Haskell is a functional programming language with advanced features of type system mainly for the research of this field. But it also provides an interface to read and write pointers.

import Control.Monad ((<=<))
import Data.Bits (shiftR)
import Data.Traversable (mapM)
import Data.Word (Word32)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Ptr (castPtr)
import Foreign.Storable (peek, poke)
import Prelude hiding (mapM)
import Text.Read (readMaybe)

main :: IO ()
main = mapM_ (mapM (print <=< fastInvSqrt) . readMaybe) . lines =<< getContents

fastInvSqrt :: Float -> IO Float
fastInvSqrt x =
  alloca $ \ptr -> do
    poke ptr x
    i <- peek (castPtr ptr)
    poke (castPtr ptr) $ 0x5f3759df - (i :: Word32) `shiftR` 1
    y <- peek ptr
    return $ y * (1.5 - 0.5 * x * y * y)

Erlang

How cool the binary pattern match is. I have never experienced this language but I think I want to learn more. The test results show that the performance of this code is bad. If you have any idea of improving the performance, please let me know from the issue tracker.

-module(fastinvsqrt).
-export([main/0, fast_inv_sqrt/1]).

main() ->
  case io:get_line("") of
    eof -> ok;
    {error, _} -> ok;
    Line -> catch io:format("~f~n", [fast_inv_sqrt(parse_float(Line))]), main()
  end.

fast_inv_sqrt(X) ->
  <<I:32/integer>> = <<X:32/float>>,
  J = 16#5f3759df - (I bsr 1),
  <<Y:32/float>> = <<J:32/integer>>,
  Y * (1.5 - 0.5 * X * Y * Y).

parse_float(Str) ->
  case string:to_float(Str) of
    {error, _} ->
      case string:to_integer(Str) of
        {error, _} -> error;
        {Int, _} -> float(Int)
      end;
    {Float, _} -> Float
  end.

Author

itchyny (https://github.com/itchyny)

Some codes were written by other people. Thank you so much.

• The credit of the Smalltalk version goes to sumim #2.

fastinvsqrt's Issues

Language Addition: Solidity

Current findings

Currently there's a few implementations that work using bytes32 or by splitting the sign, exponent & mantissa for bitwise operations however no full implementation appears to be complete.

https://github.com/BANKEX/solidity-float-point-calculation/blob/master/contracts/FloatMath.sol
https://github.com/cag/solidity-float-point-calculation

Current suggestions is to use the Babylonian method of finding a root:

function sqrt(uint x) returns (uint y) {
        if (x == 0) return 0;
        else if (x <= 3) return 1;
        uint z = (x + 1) / 2;
        y = x;
        while (z < y)
        /// @why3 invariant { to_int !_z = div ((div (to_int arg_x) (to_int !_y)) + (to_int !_y)) 2 }
        /// @why3 invariant { to_int arg_x < (to_int !_y + 1) * (to_int !_y + 1) }
        /// @why3 invariant { to_int arg_x < (to_int !_z + 1) * (to_int !_z + 1) }
        /// @why3 variant { to_int !_y }
        {
            y = z;
            z = (x / z + z) / 2;
        }
    }```

Implementation for Delphi

手元にgit入れてなく、hatenaアカウントももっていないため、Issueで失礼します。

試しにDelphiのコード書いてみました・・・
というかほぼCそのままの写しで済んでしまったわけですが・・・。

program Project1;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  System.SysUtils;

function FastInvSqrt(x: single): single;
var
  i: integer;
  y: single;
begin
  i := PInteger(@x)^;
  i := $5f3759df - (i shr 1);
  y := PSingle(@i)^;

  Result := y * (1.5 - 0.5 * x * y * y);
end;

var
  line: string;
  f: single;
begin
  try
    { TODO -oUser -cConsole メイン : ここにコードを記述してください }
    repeat
      Readln(line);
      if line = '' then Break;
      if Single.TryParse(line, f) then begin
        Writeln(FastInvSqrt(f));
      end;
    until false;
  except
    on E: Exception do
      Writeln(E.ClassName, ': ', E.Message);
  end;
end.

Use `f32::to_bits` and `f32::from_bits` for Rust

Hi there, just saw this repo. To avoid unsafe in Rust you can use the mentioned methods! See this StackOverflow Q&A for more information.

Matlab implementation

Hi,

do you plan in implementing it also in Matlab?

Edit: Already managed it on my own!

Thanks.

Roman