fast reverse/inverse square root (Quake III style)

drsqrt.f

      pure function drsqrt(rev) result(res)
      use, intrinsic :: iso_fortran_env, only : int64
      implicit none
      real*8,intent(in) :: rev
      real*8            :: rsq
      integer(int64)    :: isq
      real*8            :: res
      real*8,parameter  :: onehalf    = 0.5
      real*8,parameter  :: threehalfs = 1.5
      integer(int64)    :: magic
      parameter(magic = Z'5fe6eb50c7b537a9')
      equivalence(rsq,isq)

      rsq = rev
      isq = magic - ishft(isq,-1)
      rsq = rsq * (threehalfs - onehalf * rev * rsq**2)
      res = rsq

      return
      end function drsqrt

rsqrt.f

      pure function rsqrt(rev) result(res)
      use, intrinsic :: iso_fortran_env, only : int32
      implicit none
      real,intent(in) :: rev
      real            :: rsq
      integer(int32)  :: isq
      real            :: res
      real,parameter  :: onehalf    = 0.5
      real,parameter  :: threehalfs = 1.5
      integer(int32)  :: magic
      parameter(magic = Z'5f375a86')
      equivalence(rsq,isq)

      rsq = rev
      isq = magic - ishft(isq,-1)
      rsq = rsq * (threehalfs - onehalf * rev * rsq**2)
      res = rsq

      return
      end function rsqrt