c - Floating Point Square Root Reciprocal Method Correct Rounding -

i have implemented 32-bit ieee-754 floating point square root using newton-raphson method (in assembly) based upon finding reciprocal of square root. using round-to-nearest rounding method. square root method accepts normalized values , zeros, no denormalized values or special values (nan, inf, etc.)

i wondering how can achieve correct rounding (with assembly instructions) results correct (to ieee-754) inputs? basically, know how test if results correct, want adjust algorithm below obtain correctly rounded results. instructions should add algorithm?

see: determining floating point square root more information

thank you!

why not square result, , if it's not equal input, add or subtract (depending on sign of difference) least significant bit, square, , check whether have given better result?

better here mean less absolute difference. case tricky when "crossing" √2 mantissa, checked once , all.

edit

i realize above answer insufficient. squaring in 32-bit fp , comparing input doesn't give enough information. let's y = your_sqrt(x). compare y² x, find y²>x, subtract 1 lsb y obtaining z (y1 in comments), compare z² x , find not z²<x, but, within available bits, y²-x==x-z² - how choose between y , z? should either work bits (i guess looking for), or @ least more bits (which guess njuffa suggesting).

from comment of yours suspect on strictly 32-bit hardware, let me suppose have 32-bit 32-bit integer multiplication 64-bit result available (if not, can constructed). if take 23 bits of mantissa of y integer, put 1 in front, , multiply itself, have number that, except possible shift 1, can directly compare mantissa of x treated same way. way have 48 bits available comparison, , can decide without approximation whether abs(y²-x)≷abs(z²-x).

if not sure within 1 lsb final result (but sure not farther that), should repeat above until y²-x changes sign or hits 0. watch out edge cases, though, should cases when exponent adjusted because mantissa crosses power of 2.

it can helpful remember positive floating point numbers can correctly compared integers, @ least on machines 1.0f 0x3f800000.