The iterative procedure y(n+1)=(1+y(n)^2)/(2y(n)+x) converges to sqrt(1 + (x/2)^2) - x/2 (NB note the minus sign as opposed to the original question). Easy to prove by solving the equation for the limit value y: y=(1+y^2)/(2y+x) which transforms into a quadratic equation. It can be derived by working out the Newton-Raphson algorithm for this case.

With a starting value y(0)=1/x it converges very quickly for any x>10 and because it doesn't need the square of x there are no overflow issues for large x. Simply terminate the iterations when y(n+1)=y(n) (i.e. machine precision).

I suppose this is what is under the bonnet for most hypot-implementations