Fast InvSqrt in Max

After reading this article:

http://www.beyond3d.com/articles/fastinvsqrt/

And having some time to spare while waiting on the Sint i made a BM version of the original algorithm.

Original c version:

float math_InvSqrt (float x){
   float xhalf = 0.5f*x;
   int i = *(int*)&x;
   i = 0x5f3759df - (i >> 1);
   x = *(float*)&i;
   x = x*(1.5f - xhalf*x*x);
   return x;
}

float math_Sqrt (float x)
{
   float xhalf = 0.5f*x;
   int i = *(int*)&x;
   i = (0xbe6f0000 - i)>>1;
   x = *(float*)&i;
   x = x*(1.5f - xhalf*x*x);
   return x;
}

and my Blitzmax version:

Import "InvSqrt.c"

Extern 
	
	Function InvSqrt:Float(x:Float) = "math_InvSqrt"
	Function Sqrt:Float(x:Float) 	= "math_Sqrt"
	
End Extern

Function InvSqrtKnotz:Float(x:Float)

	Local xhalf:Float = 0.5*x
	Local i:Int = (Int Ptr(Varptr(x)))[0]
	i = $5f3759df - (i Shr 1)
	x = (Float Ptr(Varptr(i)))[0]
	x = x*(1.5 - xhalf*x*x)
	Return x

End Function

Function InvSqrtAntMan:Float(x:Float)

	Local i = $5f3759df - (Int Ptr(Varptr(x))[0] Shr 1)
	Local ox:Float = (Float Ptr(Varptr(i))[0])
	Return  ox *(1.5 - ((0.5*x)*ox*ox))
	
End Function

Function InvSqrtFredborg:Float(x:Float)

	Local i:Int = $5f3759df - (Int Ptr(Varptr(x))[0] Shr 1)
	Return  (Float Ptr(Varptr(i))[0]) *(1.5 - ((0.5*x)*(Float Ptr(Varptr(i))[0])*(Float Ptr(Varptr(i))[0])))
	
End Function

Function InvSqrtSOS:Float(x:Float)

    Local xhalf:Float = 0.5*x
	Local i:Int = $5f3759df - (Int Ptr(Varptr(x))[0] Shr 1)
	Local ox:Float = (Float Ptr(Varptr(i))[0])
	Return ox *(1.5 - (xhalf*ox*ox))
	
End Function

Function InvSqrtAza:Float(x:Float)

	Local xhalf:Float = 0.5*x
	Local ip:Int Ptr=Int Ptr(Varptr(x))
	ip[0]=$5f3759df-(ip[0] Shr 1)
	Return x *(1.5 - (xhalf*x*x))

End Function

Function InvSqrtPlain:Float(x:Float)

	Return 1.0/Sqr(x)
	
End Function

Function SqrtKnotz:Float(x:Float)

	Local xhalf:Float = 0.5*x
	Local i:Int = (Int Ptr(Varptr(x)))[0]
	i = ($be6f0000 - i) Shr 1
	x = (Float Ptr(Varptr(i)))[0]
	x = x*(1.5 - xhalf*x*x)
	Return x

End Function

Function SqrtAntMan:Float(x:Float)

	Local i = ($be6f0000 - (Int Ptr(Varptr(x))[0])) Shr 1
	Local ox:Float = (Float Ptr(Varptr(i))[0])
	Return  ox *(1.5 - ((0.5*x)*ox*ox))
	
End Function

Function SqrtFredborg:Float(x:Float)

	Local i:Int = ($be6f0000 - (Int Ptr(Varptr(x))[0])) Shr 1
	Return  (Float Ptr(Varptr(i))[0]) *(1.5 - ((0.5*x)*(Float Ptr(Varptr(i))[0])*(Float Ptr(Varptr(i))[0])))
	
End Function


Function SqrtSOS:Float(x:Float)

    Local xhalf:Float = 0.5*x
	Local i:Int = ($be6f0000 - (Int Ptr(Varptr(x))[0])) Shr 1
	Local ox:Float = (Float Ptr(Varptr(i))[0])
	Return ox *(1.5 - (xhalf*ox*ox))
	
End Function

Function SqrtAza:Float(x:Float)

	Local xhalf:Float = 0.5*x
	Local ip:Int Ptr=Int Ptr(Varptr(x))
	ip[0]=($be6f0000 - ip[0]) Shr 1
	Return x *(1.5 - (xhalf*x*x))

End Function

Function SqrtPlain:Float(x:Float)

	Return Sqr(x)
	
End Function

Function Benchmark(label:String, func:Float(x:Float))

	Local startTime:Long = MilliSecs()
	For Local ii:Int=1 To 100000000
		func(Float(ii))
	Next
	Local endTime:Long = MilliSecs() - startTime
	Print label + ":~t~t" + endTime + "ms"

End Function

Print "~nInvSqrt:"
Benchmark("Raw C",InvSqrt)
Benchmark("Knotz",InvSqrtKnotz)
Benchmark("Artman",InvSqrtAntman)
Benchmark("Fredborg",InvSqrtFredborg)
Benchmark("SculptureOfSoul",InvSqrtSOS)
Benchmark("Azathoth",InvSqrtAza)
Benchmark("1/Sqr()",InvSqrtPlain)

Print "~nSqrt:"
Benchmark("Raw C",Sqrt)
Benchmark("Knotz",SqrtKnotz)
Benchmark("Artman",SqrtAntman)
Benchmark("Fredborg",SqrtFredborg)
Benchmark("SculptureOfSoul",SqrtSOS)
Benchmark("Azathoth",SqrtAza)
Benchmark("Sqr()",SqrtPlain)

Some more information can be found here:

http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf

I've included the entries from the posts in the code. I've also added a fast implementation of the regular Sqr().

Have a nice weekend.