One of our numerical analysis(NA) professors had worked at NASA, at IBM and had enormous real-world experience in NA. Several weeks into the class he challenged us to write a correct square root routine in FORTRAN. Over the next two weeks submissions were presented in class, where our professsor showed how each had multiple critically fatal errors. He let it stew awhile: two weeks before semester's end he presented and explicated to us a correct implementation of the square root, which turned out to be far more lengthy and involved than we had imagined. He said there was always the possibility of more bugs being present. That was really an eye-opening class.
I saw that done years ago in assembler for UNIVAC mainframes. Tear apart the exponent and mantissa. Shift mantissa if exponent is odd. Shift exponent right one bit to get square root of exponent. Look up starting guess from table using a few high bits of mantissa. Run two iterations of Newton's method. Reassemble mantissa and exponent. Normalize. Done.
musl sqrt is here, about hundred lines of code with nice explanation on top. Considering that it has 1993 copyright stamp, I'd imagine it is fairly canonical implementation
