Skip to content

Instantly share code, notes, and snippets.

@certik
Last active Oct 13, 2021
Embed
What would you like to do?
Testing 1/sqrt(x) implementations
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment