One more implication: since the quality of the initial square root doesn't really much matter, instead of

    result = sqrt(to_float(parts))
    a, b = split(result)
    parts = add_on(-a*a, parts)
    parts = add_on(-2.0*a*b, parts)
    parts = add_on(-b*b, parts)
    x = to_float(parts)
    result += x / (2.0 * result)

at the end, it should work just as well (in fact, probably a little better, due to getting more beneficial cancellation) to do:

    a = parts[0] - 1.0
    result = sqrt(a)
    x, y = split(result)
    result += (a - x*x - 2.0*x*y - y*y + parts[1]) / (2.0 * result)

at the end. Although this crucially relies on the doing the last-line chain of subtracts and adds "left to right".