I should point out the "denormalized" numbers completely blow out the "SQRT hack", regardless of the FP format. The function asinh(x) acts like exp(x)/2 for large x, and like x for small x -- i.e., like a "normalized" FP number for large x (albeit with a radix of e), and like a "denormalized" number for small x. So asinh(x) is the "smoothest" floating point format that includes denormalized numbers. So I tested Newton given a first approximation sinh(asinh(N)/2) to sqrt(N). In order to compute a full precision (in GNU Common Lisp) IEEE sqrt for N in the range 2^1 up to 2^512, exactly _five_ Newton iterations sufficed. However, when N started falling in the "denomalized" range -- i.e., N<1, the number of Newton iterations skyrocketed, basically requiring N Iters reqd 2^0 6 2^-10 11 2^-20 16 2^-30 21 2^-40 26 2^-50 31 etc. So, if you care about the square roots of denormalized numbers, you will have to treat them as a special case.