Skip to content

Instantly share code, notes, and snippets.

@dorchard
Last active February 7, 2024 17:08
Show Gist options
  • Save dorchard/c9c0d2307f5bf875b0c664a0f3c34a44 to your computer and use it in GitHub Desktop.
Save dorchard/c9c0d2307f5bf875b0c664a0f3c34a44 to your computer and use it in GitHub Desktop.
Fast inverse square root in Haskell - Two ways
import Data.Bits
import Data.Word
import Foreign ( Ptr, new, castPtr, Storable(peek, poke) )
import System.IO.Unsafe
import Unsafe.Coerce
qsqrt :: Float -> Float
qsqrt x = 1 / qRsqrt x
qsqrt' :: Float -> Float
qsqrt' x = 1 / qRsqrt' x
-- Quake 3 inverse square root, using C-style pointer casting
qRsqrt :: Float -> Float
qRsqrt number = unsafePerformIO $ do
let x2 = number * 0.5
y <- new number
let i = castPtr y :: (Ptr Word32)
magic <- peek i
let magic' = 0x5f3759df - (magic `shiftR` 1)
poke i magic'
y' <- peek y
return $ y' * (1.5 - (x2 * y' * y'))
-- Quake 3 inverse square root, using unsafe coerce
qRsqrt' :: Float -> Float
qRsqrt' number =
let x2 = number * 0.5
magic = unsafeCoerce number :: Word32
magic' = 0x5f3759df - (magic `shiftR` 1)
number' = unsafeCoerce magic' :: Float
in (number' * (1.5 - (x2 * number' * number')))
-- Compare sqrt 1 and sqrt 2 for both methods
test () = map (\x -> (qsqrt x, qsqrt' x, qsqrt x == qsqrt' x)) (take 2 [1.0..])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment