> ...
> one at a time subtract a**2 (an argument squared) in descending
> order of magnitude
> ...

But that doesn't really help unless the sum of squares was computed without care to begin with. Can do as well by skipping that but instead computing the original sum of squares in increasing order of magnitude.

Cheapest way I know of that "seemingly always" reproduces the Decimal result (when that's rounded back to float) combines fsum(), Veltkamp splitting, and the correction trick from the paper:

    def split(x, T27=ldexp(1.0, 27)+1.0):
        t = x * T27
        hi = t - (t - x)
        lo = x - hi
        assert hi + lo == x
        return hi, lo

    # result := hypot(*xs)
    parts = []
    for x in xs:
        a, b = split(x)
        parts.append(a*a)
        parts.append(2.0*a*b)
        parts.append(b*b)
    result = sqrt(fsum(parts))
    a, b = split(result)
    parts.append(-a*a)
    parts.append(-2.0*a*b)
    parts.append(-b*b)
    x = fsum(parts) # best float approx to sum(x_i ** 2) - result**2
    result += x / (2.0 * result)