Cody の SQRT

program sqrt2
    implicit none
    real(16) :: x, y, z, y2, y3, f, e
    integer :: i, j, k
    REAL(16), PARAMETER :: sqrt05 = SQRT(0.5q0)
   DO j = 1, 1000
    CALL RANDOM_NUMBER(x)
    x = x * 1000.0q0 
!    print *, 'x=', x
    f = fraction(x)
    i = exponent(x)
!    print *, f, i
 !   print '(b64.64)', x
    y = 0.41731q0 + 0.59016q0 * f
    z = y + f / y 
    y = set_exponent(z, exponent(z) - 2) + f / z ! 0.25 * z + f / z
    DO k = 1, 2 ! 0 single, 1 double, 2 quadruple
     y = y + f / y
     y = set_exponent(y, exponent(y) - 1) ! 0.5 * y
    END DO
    IF (MOD(i, 2) == 0) THEN
     e = set_exponent(y, i / 2)
    ELSE
     e = set_exponent(y * sqrt05, (i + 1) / 2)
    END IF 
    print *, e - sqrt(x)
  end do
end program sqrt2