Add fast inverse square root algorithm (#647)

* Add fast inverse square root algorithm

* PR improvements see comments in #647

* Update math/binary/fast_inverse_sqrt.go file docstring

Co-authored-by: Taj <[email protected]>

* Update math/binary/fast_inverse_sqrt.go function documentation

Co-authored-by: Taj <[email protected]>

* Fix function documentation grammar

Co-authored-by: Taj <[email protected]>

---------

Co-authored-by: Taj <[email protected]>

math/binary/fast_inverse_sqrt.go

@@ -0,0 +1,29 @@
+// Calculating the inverse square root
+// [See more](https://en.wikipedia.org/wiki/Fast_inverse_square_root)
+
+package binary
+
+import (
+	"math"
+)
+
+// FastInverseSqrt assumes that argument is always positive,
+// and it does not deal with negative numbers.
+// The "magic" number 0x5f3759df is hex for 1597463007 in decimals.
+// The math.Float32bits is alias to *(*uint32)(unsafe.Pointer(&f))
+// and math.Float32frombits is to *(*float32)(unsafe.Pointer(&b)).
+func FastInverseSqrt(number float32) float32 {
+	var i uint32
+	var y, x2 float32
+	const threehalfs float32 = 1.5
+
+	x2 = number * float32(0.5)
+	y = number
+	i = math.Float32bits(y)   // evil floating point bit level hacking
+	i = 0x5f3759df - (i >> 1) // magic number and bitshift hacking
+	y = math.Float32frombits(i)
+
+	y = y * (threehalfs - (x2 * y * y)) // 1st iteration of Newton's method
+	y = y * (threehalfs - (x2 * y * y)) // 2nd iteration, this can be removed
+	return y
+}

math/binary/sqrt.go

@@ -7,18 +7,4 @@
 
 package binary
 
-import (
-	"math"
-)
-
-const threeHalves = 1.5
-
-func Sqrt(n float32) float32 {
-	var half, y float32
-	half = n * 0.5
-	z := math.Float32bits(n)
-	z = 0x5f3759df - (z >> 1) // floating point bit level hacking
-	y = math.Float32frombits(z)
-	y = y * (threeHalves - (half * y * y)) // Newton's approximation
-	return 1 / y
-}
+func Sqrt(n float32) float32 { return 1 / FastInverseSqrt(n) }

math/binary/sqrt_test.go

@@ -10,7 +10,7 @@ import (
 	"testing"
 )
 
-const epsilon = 0.2
+const epsilon = 0.001
 
 func TestSquareRootCalculation(t *testing.T) {
 	tests := []struct {