Testing 1/sqrt(x) implementations

test_rsqrt.f90

program test_rsqrt
! Prints:
! $ gfortran test_rsqrt.f90 && ./a.out
!   0.81649658092772603     
!   0.81649658092772626     
!   0.81649658092772603     
!   0.81649658092772603     

implicit none

integer, parameter :: sp = kind(0.0)
integer, parameter :: dp = kind(0.d0)
real(dp) :: x, r

x = 1.5_dp
r = 1/sqrt(real(x,sp))
r = r * (1.5_dp - 0.5_dp * x * r**2)
r = r * (1.5_dp - 0.5_dp * x * r**2)
print *, r

r = 1/sqrt(real(x,sp))
r = r * (0.0625_dp*x**4*r**8 - 0.5625_dp*x**3*r**6 + 1.6875_dp*x**2*r**4 - 2.4375_dp*x*r**2 + 2.25_dp)
print *, r

r = 1/sqrt(real(x,sp))
r = r * (0.375_dp*x**2*r**4 - 1.25_dp*x*r**2 + 1.875_dp)
print *, r

r = rsqrt(x)
print *, r

contains

    real(dp) function rsqrt(s) result(x)
    real(dp), intent(in) :: s
    real(dp) :: y, a, b
    real(sp) :: q
    q = real(s, sp)
    q = 1/sqrt(q)
    x = real(q, dp)
    y = s * x * x
    a = y * 0.375_dp
    a = a * y
    b = y * 1.25_dp
    b = b - 1.875_dp
    y = a - b
    x = x * y
    end function

end program